On Sat, Jan 31, 2015 at 5:41 PM, Patrice Dumas <[email protected]> wrote: > Hello, > > grass 6.4.4 on Centos 6 (from EPEL). > > Before I try to turn this into a bug I'd like to be sure that I am not > doing something stupid with projections (happened in the past...). > > I have a Lambert equal area location and a lon-lat location, the lon-lat > location covers the whole earth, the Lambert does not. I would like to > project some points from the lon-lat location to the Lambert one, > selecting only those points in a specific region on the Lambert map > which is much smaller than the whole Lambert map. > > To that aim, I first do a "grid" vector over the region of interest in > the Lambert location with "v.mkgrid map=region_countries_whole_grid > grid=40,40" (a rather loose grid such that it is easier to project on > different locations). I import that grid in the lon-lat location, use > it to select points and then, back in the Lambert location use v.proj to > import the points. Except that it gives me > > > v.proj --verbose input=lon_lat_demand_relocated > location=whole_region_demand_map_asiamed_cities
Is lon_lat_demand_relocated really the point selection? > Input Projection Parameters: +proj=longlat +a=6378137 +rf=298.257223563 > +no_defs +towgs84=0.000,0.000,0.000 > Input Unit Factor: 1 > Output Projection Parameters: +proj=laea +lat_0=55 +lon_0=20 +x_0=0 +y_0=0 > +no_defs +a=6370997 +b=6370997 > Output Unit Factor: 1 > Reprojecting primitives: WARNING: pj_transform() failed: latitude or > longitude exceeded limits > ERROR: Error in pj_do_transform > > > Visually it seems that all the points are on the "grid" vector (which is > obviously not a regular grid once projected on the lon-lat location). > What is very strange is that projecting the "grid" vector back to the > Lambert location works, but not projecting the points back?? >From your description, the workflow is # create grid in LAEA v.mkgrid map=region_countries_whole_grid grid=40,40" # reproject grid to latlong v.proj in=region_countries_whole_grid location= mapset= # select points with grid in latlong v.select ainput=points binput=region_countries_whole_grid atype=point btype=area operator=overlap output=points_selected # reproject selected points to LAEA v.proj in=points_selected location=whole_region_demand_map_asiamed_cities mapset= # test: project grid back to LAEA v.proj in=region_countries_whole_grid location=whole_region_demand_map_asiamed_cities mapset= > > Am I doing something silly, or does it looks like a bug? > If you are sure that v.proj fails with the selected points, it is a bug. v.proj will likely fail with all points if the points are distributed all over the world, that would be correct. > > The Lambert location map comes from HYDRO1K (the eu basins): > http://edcftp.cr.usgs.gov/pub/data/gtopo30hydro/ > https://lta.cr.usgs.gov/HYDRO1K > https://lta.cr.usgs.gov/HYDRO1KReadMe > I used the grass helper to setup the region and it didn't set a datum > (but it set an ellipsoid). In the .prj file that comes with vectors, > there is > > PROJCS["Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area",GEOGCS["GCS_Sphere_ARC_INFO",DATUM["D_Sphere_ARC_INFO",SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",20.0],PARAMETER["Latitude_Of_Origin",55.0],UNIT["Meter",1.0]] > but I couldn't find the "D_Sphere_ARC_INFO" datum in the datum proposed > by grass. It's a simple sphere with radius 6370997.0, without datum. Markus M > The other information from the .prj file were correctly > extracted by the grass helper, as far as I can tell. > > -- > Pat > _______________________________________________ > grass-user mailing list > [email protected] > http://lists.osgeo.org/mailman/listinfo/grass-user _______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
