Opened 5 years ago

Last modified 5 years ago

#5415 new defect

MorphToESRI creates an ambiguous WKT for EPSG:2903 (NM Central HARN , NAD83)

Reported by: roomthily Owned by: warmerdam
Priority: lowest Milestone:
Component: OGR_SRS Version: 1.9.1
Severity: normal Keywords:
Cc:

Description

We've been building tile indexes for MapServer? and, while testing, noticed an issue when generating a file for EPSG:2903 (NAD83(HARN) / New Mexico Central (ftUS)).

The errors include failing to render the geometries in QGIS and MapServer? dumping the memory map. When looking at the properties in QGIS, the CRS is marked as 2258 instead of 2903. The shapefile is created using the 2903 spatial reference with the morphtoesri output in the PRJ file.

The WKTs:

>>> from osgeo import osr
>>> s2258 = osr.SpatialReference()
>>> s2258.ImportFromEPSG(2258)
0
>>> s2903 = osr.SpatialReference()
>>> s2903.ImportFromEPSG(2903)
0
>>> s2258.ExportToWkt()
'PROJCS["NAD83 / New Mexico Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",31],PARAMETER["central_meridian",-106.25],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",1640416.667],PARAMETER["false_northing",0],UNIT["US survey foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2258"]]'
>>> s2903.ExportToWkt()
'PROJCS["NAD83(HARN) / New Mexico Central (ftUS)",GEOGCS["NAD83(HARN)",DATUM["NAD83_High_Accuracy_Reference_Network",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6152"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4152"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",31],PARAMETER["central_meridian",-106.25],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",1640416.667],PARAMETER["false_northing",0],UNIT["US survey foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2903"]]'
>>> mte2903 = osr.SpatialReference()
>>> mte2903.ImportFromEPSG(2903)
0
>>> mte2903.MorphToESRI()
0
>>> mte2903.ExportToWkt()
'PROJCS["NAD83_HARN_New_Mexico_Central_ftUS",GEOGCS["GCS_NAD83(HARN)",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",31],PARAMETER["central_meridian",-106.25],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",1640416.667],PARAMETER["false_northing",0],UNIT["Foot_US",0.30480060960121924]]'
>>> 

Example of the MapServer? error:

*** glibc detected *** XXXX           : free(): invalid next size (fast): 0x00007f242e494610 ***
======= Backtrace: =========/lib/x86_64-linux-gnu/libc.so.6(+0x7e626)[0x7f243aa2a626]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msShapefileClose+0x4b)[0x7f2422816e6b]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msSHPLayerClose+0x15)[0x7f2422816e95]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msDrawRasterLayerLow+0x24b)[0x7f242284bbdb]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msDrawRasterLayer+0xd5)[0x7f24228610a5]
/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msDrawLayer+0x465)[0x7f24228679b5]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(msDrawMap+0x434)[0x7f2422868644]/opt/modwsgi/XXXX_env/lib/python2.7/site-packages/MapScript-6.0.3-py2.7-linux-x86_64.egg/_mapscript.so(+0x49ebc)[0x7f24227c5ebc]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x53a5)[0x7f2437649845]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x614b)[0x7f243764a5eb]/usr/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x855)[0x7f2437614605]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x5420)[0x7f24376498c0]
/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x614b)[0x7f243764a5eb]/usr/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x855)[0x7f2437614605]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x5420)[0x7f24376498c0]
/usr/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x855)[0x7f2437614605]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x5420)[0x7f24376498c0]/usr/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x855)[0x7f2437614605]/usr/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x5420)[0x7f24376498c0]/usr/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x855)[0x7f2437614605]

We're using osgeo 1.9.1 (python 2.7), MapServer? 6.0.3, and QGIS 1.7.2. Of course, ArcGIS 10 has no problem rendering the shapefile. And there's obviously no issue with creating and using a shapefile in EPSG:2258.

The simplest solution is just not to use MorphToESRI with 2903.

Attachments (1)

tkt5415_tile2903shp.zip (777.4 KB) - added by roomthily 5 years ago.

Download all attachments as: .zip

Change History (8)

comment:1 Changed 5 years ago by Kyle Shannon

Could you attach a small dataset that fails to render? It'd be helpful.

comment:2 Changed 5 years ago by Jukka Rahkonen

Hi,

QGIS 1.7.2 is probably using pretty old GDAL. Could you try how the current QGIS 2.2 behaves?

comment:3 Changed 5 years ago by Even Rouault

I think there are different defects, in MapServer? and in QGIS. For MapServer? I doubt the SRS is the cause of the crash you see... As far as QGIS is concerned, I think it bases its recognition of SRS on the proj.4 string and EPSG:2258 and EPSG:2903 have the same proj.4 string :

  • before morphing to ESRI : '+proj=tmerc +lat_0=31 +lon_0=-106.25 +k=0.9999 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs '
  • after morphing to ESRI : '+proj=tmerc +lat_0=31 +lon_0=-106.25 +k=0.9999 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs '

I'm not sure if the difference between both EPSG code is that significant : NAD83 vs NAD83(HARN) ?

Well, looking at https://trac.osgeo.org/proj/wiki/HarnGrids, it seems that we could possibly turn the NAD83(HARN) into "+datum=NAD83 +nadgrids=azhpgn.gsb", but as that grid is not in the base proj.4 grids, I'm not sure this is a good idea.

Changed 5 years ago by roomthily

Attachment: tkt5415_tile2903shp.zip added

comment:4 Changed 5 years ago by roomthily

The current QGIS 2.2 can render the file. The CRS is empty for it, though.

As far as the MapServer? issue being a different bug, it's possible but the only difference that we made between the WMS request against the original, non-functional shapefile and a functional shapefile was removing MorphToESRI from the build process.

comment:5 Changed 5 years ago by Even Rouault

The MapServer? issue should be reported on MapServer? bug tracker, with mapfile, data and request that triggers the crash.

comment:6 Changed 5 years ago by roomthily

Well, sure. And you can police your forum as you see fit.

But in this case, it might be a bit unrealistic to assume that a platform higher up the toolchain will correctly handle an issue like this. Especially if it's not a known issue and clearly not one bothering a whole heck of a lot of people. So in the interest of finding absolutely no information anywhere and feeling a bit community minded, I thought perhaps having some information about the cause and the symptoms at the source would be beneficial. My bad.

comment:7 Changed 5 years ago by Even Rouault

roomthily, you've misunderstood me. I'm also a MapServer? developer. My feeling is that the crash in MapServer? is due to a bug in MapServer?, not in OGR. I can be wrong, but without any direct mean of replicating, it is hard to tell. And the easier it is to replicate, the more likely it is that someone will spend a bit of time to analyze to see what happens.

I will also note that you don't use the latest stable versions of software. Perhaps upgrading to MapServer? 6.4 would solve the crash. And QGIS 2.2 is now out (1.7.3 is pretty old). I'm not sure that this will solve the confusion between the 2 SRS, but better be up-to-date.

As far as the global issue, I'm not sure how and if that can be fixed. SRS issues over the whole software stack tend to be complex ones. I already shared some thoughts above. And on my volunteer time... You can find a list of service providers for OSGeo software at http://www.osgeo.org/search_profile

Note: See TracTickets for help on using tickets.