ogr2ogr segfault with LCC target projection in ESRI style

[neteler]

Frank,

reprojecting a SHAPE file from LatLong?/WGS84 to North Carolina LCC leads to a segfault when the target .prj file is in ESRI style, otherwise it works:

When .prj looks like this, it fails:

g.proj -wef
PROJCS["Lambert Conformal Conic",GEOGCS["grs80",DATUM["D_North_American_1983",SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",34.33333333333334],PARAMETER["standard_parallel_2",36.16666666666666],PARAMETER["latitude_of_origin",33.75],PARAMETER["central_meridian",-79],PARAMETER["false_easting",609601.22],PARAMETER["false_northing",0],UNIT["Meter",1]]

ogr2ogr -s_srs epsg:4326 -t_srs "`g.proj -wef`" tiger_landmark_points_NCSPM.shp tgr37183lpt.shp
OGR: Auto register /usr/local/lib/gdalplugins/ogr_GRASS.so using RegisterOGRGRASS.
OGR: OGROpen(tgr37183lpt.shp/0x8052b08) succeeded as ESRI Shapefile.
OGRCT: Source: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
ERROR 6: Failed to initialize PROJ.4 with `+ellps=GRS80 +units=m +no_defs '.
OGRCT: Target: +ellps=GRS80 +units=m +no_defs
Failed to create coordinate transformation between the
following coordinate systems.  This may be because they
are not transformable, or because projection services
(PROJ.4 DLL/.so) could not be loaded.
Segmentation fault

while it works like this:

# no ESRI style
g.proj -wf
PROJCS["Lambert Conformal Conic",GEOGCS["grs80",DATUM["North_American_Datum_1983",SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",34.33333333333334],PARAMETER["standard_parallel_2",36.16666666666666],PARAMETER["latitude_of_origin",33.75],PARAMETER["central_meridian",-79],PARAMETER["false_easting",609601.22],PARAMETER["false_northing",0],UNIT["Meter",1]]

ogr2ogr -s_srs epsg:4326 -t_srs "`g.proj -wf`" tiger_landmark_points_NCSPM.shp tgr37183lpt.shp
OGR: Auto register /usr/local/lib/gdalplugins/ogr_GRASS.so using RegisterOGRGRASS.
OGR: OGROpen(tgr37183lpt.shp/0x8052b58) succeeded as ESRI Shapefile.
OGRCT: Source: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
OGRCT: Target: +proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
Shape: 131 features read on layer 'tgr37183lpt'.

Apparently the ESRI morph fails in this case.

The source SHAPE is attached.

Best, Markus

[neteler]

TIGER extracted points in LatLong?

[neteler]

Note: Hobu mentioned

http://trac.osgeo.org/gdal/changeset/10845, at , by warmerdam: added LCC 2SP morph test per bug 187

in an IRC chat.

[warmerdam]

Mateusz,

Could you please try and reproduce this and if so fix it?

[mloskot]

The segmentation fault was caused by dereferencing of null pointer in ogr2ogr.cpp and is not directly related to this problem with ESRI style of projection def.

I applied the fix in r11636.

[mloskot]

I applied one more improvement for handling ESRI style of WKT from file (r11638).

The OGRSpatialReference::SetFromUserInput?() function also accepts name of file in which projection definition is stored. There was missing test of ESRI:: prefix in WKT string read from given file.

NOTE: According to the  manual, the ESRI:: prefix is required in order to trigger morphToESRI() on WKT string given in ESRI style. Otherwise, the projection definition may be imported incorrectly, incompletely or an error can be thrown.

Markus,

As I see, g.proj -wef command returns ESRI WKT string without the prefix mentioned above, and that's why the transformation failes (note, the segfault was caused by yet another issue, already fixed). The WKT in ESRI style should be passed to OGR as follows:

ESRI::PROJCS["Lambert Conformal Conic",GEOGCS["grs80",DATUM["D_North_American_1983",SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",34.33333333333334],PARAMETER["standard_parallel_2",36.16666666666666],PARAMETER["latitude_of_origin",33.75],PARAMETER["central_meridian",-79],PARAMETER["false_easting",609601.22],PARAMETER["false_northing",0],UNIT["Meter",1]]

and your command should read as:

ogr2ogr -s_srs epsg:4326 -t_srs "ESRI::`g.proj -wef`" tiger_landmark_points_NCSPM.shp tgr37183lpt.shp

Could you check if this command works for you?

I tested this approach (with the prefix) and it works well for me. I'm going to attach a few files with my tests.

[mloskot]

One additional note to my analysis above. If I create OGRSpatialReference instance and import the ESRI-style WKT without the ESRI:: prefix (using SetFromUserInput?() function):

PROJCS["Lambert Conformal Conic",GEOGCS["grs80",DATUM["D_North_American_1983",SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",34.33333333333334],PARAMETER["standard_parallel_2",36.16666666666666],PARAMETER["latitude_of_origin",33.75],PARAMETER["central_meridian",-79],PARAMETER["false_easting",609601.22],PARAMETER["false_northing",0],UNIT["Meter",1]]

and next I will export the SRS definition to PROJ.4 (using exportToProj4() function), the output I get is incomplete:

+ellps=GRS80 +units=m +no_defs 

instead of

+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs

This incomplete form of PROJ.4 definition is the main reason of the transformation failure Markus is experiencing. The transformation fails during PROJ.4 initialization using this incomplete definition, in ogrct.cpp:517

psPJTarget = pfn_pj_init_plus( pszProj4Defn );

where the returned psPJTarget is NULL.

Summarizing, the problem was caused by missing ESRI:: prefix in WKT of ESRI-style.

[mloskot]

Markus' WKT definition in ESRI style without the prefix

[mloskot]

Markus' WKT definition in ESRI style with the mandatory prefix

[mloskot]

Simple test of processing WKT in ESRI-style from string with the prefix to PROJ.4 definition

[mloskot]

Simple test of processing WKT in ESRI-style read directly from file, with the ESRI:: prefix, to PROJ.4 definition

[warmerdam]

Mateusz,

Your analysis looks good. I think r11636 should go into 1.4 branch, but r11638 (as a new feature) should likely only be in trunk.

[mloskot]

Frank,

I backported r11636 changeset to the 1.4 branch (r11646).

[mloskot]

I added two test cases validating the ticket fix (r11647).

[mloskot]

I'm closing this ticket. I believe the fix has been tested well.

Markus, please reopen it if it doesn't work for you.

[neteler]

Hi,

thanks for the fixes (still have to test).

For the record, I have made a ticket for GRASS to add the ESRI:: prefix in case the -e flag is used:

 http://wald.intevation.org/tracker/index.php?func=detail&aid=434&group_id=21&atid=204

Markus