[GRASS-user] v.proj pj_transform() failed for points but not for a vector "grid"

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
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??

Am I doing something silly, or does it looks like a bug?

Additional info (similar report and lack of datum):

I found a possibly similar issue reported:
  http://lists.osgeo.org/pipermail/grass-user/2012-June/065049.html
but no definitive answer.
In some answer to another thread that started with the same issue
erroneously,
  http://lists.osgeo.org/pipermail/grass-user/2012-October/065805.html
Markus asked to do a g.region -b, but in fact this does not work in the
Lambert location for me, as there is no datum information. Lack of
datum may be a cause of the v.proj issue above, although I doubt it.

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. The other information from the .prj file were correctly
extracted by the grass helper, as far as I can tell.

--
Pat

On Sat, Jan 31, 2015 at 5:41 PM, Patrice Dumas <pertusus@free.fr> 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
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

On Mon, Feb 02, 2015 at 09:29:18AM +0100, Markus Metz wrote:

On Sat, Jan 31, 2015 at 5:41 PM, Patrice Dumas <pertusus@free.fr> wrote:
>
> v.proj --verbose input=lon_lat_demand_relocated location=whole_region_demand_map_asiamed_cities

Is lon_lat_demand_relocated really the point selection?

It was not. I couldn't see that visually, but in fact I had relocated
some points after the v.select and the visual inspection was not enough
to check that they were out of the "grid".

>From your description, the workflow is

It was almost that, as said above, with a relocation of points that
moved some of them outside of the grid. I readded a last v.select and now
it works correctly.

Thanks a lot, and sorry for the noise.

--
Pat