Opened 18 months ago

Closed 18 months ago

Last modified 18 months ago

#6523 closed defect (fixed)

Incorrect parameters read for Mercator projection in adf-style prj file

Reported by: Mike Toews Owned by: warmerdam
Priority: normal Milestone: 1.11.5
Component: default Version: unspecified
Severity: normal Keywords:
Cc:

Description

Reading "NZ 250m ESRI ASCII grid" from this website (563 MB) with QGIS shows the raster in the wrong location. The .prj file appears to be a copy of the prj.adf file from the "NZ 250m ESRI binary grid" download option, and looks like this:

Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters    
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

Further information for the dataset describes the projection as EPSG:3994. However, with GDAL 2.0.2:

$ gdalinfo nzbathymetry_2016_ascii-grid.txt
Driver: AAIGrid/Arc/Info ASCII Grid
Files: nzbathymetry_2016_ascii-grid.txt
       nzbathymetry_2016_ascii-grid.prj
Size is 12116, 15391
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",100],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",-41],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]
Origin = (4795705.599999999600000,-2067835.799999999800000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Corner Coordinates:
Upper Left  ( 4795705.600,-2067835.800)
Lower Left  ( 4795705.600,-5915585.800)
Upper Right ( 7824705.600,-2067835.800)
Lower Right ( 7824705.600,-5915585.800)
Center      ( 6310205.600,-3991710.800)
Band 1 Block=12116x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

Note the discrepancies for latitude_of_origin (latitude 100?), central_meridian, and false_easting. This raster displays correctly in ArcGIS 10.2.

Change History (9)

comment:1 Changed 18 months ago by Mike Toews

To confirm, the same erroneous happens with the ESRI binary grid file, which uses prj.adf. From ArcCatalog, here is the Coordinate System info:

WGS_1984_Mercator
Authority: Custom

Projection: Mercator
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 100.0
Standard_Parallel_1: -41.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563

Looking at gdal/ogr/ogr_srs_esri.cpp ~L818, the order of PARAM_0 to PARAM_3 for SetMercator looks correct (as documented by Esri for Mercator > Parameters > Workstation), and SetMercator is behaving normally with the OSR Python bindings.

Here is a reduced version of the issue with OSR through the Python bindings:

from osgeo import gdal, osr
print(gdal.VersionInfo())

sr1 = osr.SpatialReference()
sr1.ImportFromEPSG(3994)
print('sr1 PROJ.4: ' + sr1.ExportToProj4())
print('sr1 WTK: ' + sr1.ExportToPrettyWkt())

sr2 = osr.SpatialReference()
sr2.SetMercator(-41, 100, 1, 0, 0)
print('sr2 PROJ.4: ' + sr2.ExportToProj4())
print('sr2 WTK: ' + sr2.ExportToPrettyWkt())

sr3 = osr.SpatialReference()
sr3.ImportFromESRI('''\
Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters    
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)'''.split('\n'))
print('sr3 PROJ.4: ' + sr3.ExportToProj4())
print('sr3 WTK: ' + sr3.ExportToPrettyWkt())

shows

2000200
sr1 PROJ.4: +proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
sr1 WTK: PROJCS["WGS 84 / Mercator 41",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",-41],
    PARAMETER["central_meridian",100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","3994"]]
sr2 PROJ.4: +proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
sr2 WTK: PROJCS["unnamed",
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",-41],
    PARAMETER["central_meridian",100],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
sr3 PROJ.4: +proj=merc +lon_0=0 +lat_ts=100 +x_0=-41 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sr3 WTK: PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",100],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",-41],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]

comment:2 Changed 18 months ago by Jukka Rahkonen

Could you prepare a small test data with prj.adf and some rows of ascii grid data and attach them to the ticket? I was about to have a look but downloading the 563 MB zip from New Zealand to Finland would have taken 4 hours.

comment:3 in reply to:  2 Changed 18 months ago by Mike Toews

Replying to jratike80:

Could you prepare a small test data with prj.adf and some rows of ascii grid data and attach them to the ticket? I was about to have a look but downloading the 563 MB zip from New Zealand to Finland would have taken 4 hours.

The reduced version of this issue was demonstrated in 1, so you don't need to download anything.

comment:4 Changed 18 months ago by Even Rouault

Resolution: fixed
Status: newclosed

In 34265:

importFromESRI(): fix import of multi line MERCATOR SRS (closes #6523)

comment:5 Changed 18 months ago by Even Rouault

In 34266:

importFromESRI(): fix import of multi line MERCATOR SRS (closes #6523)

comment:6 Changed 18 months ago by Even Rouault

In 34267:

importFromESRI(): fix import of multi line MERCATOR SRS (closes #6523)

comment:7 Changed 18 months ago by Even Rouault

In 34268:

importFromESRI(): fix import of multi line MERCATOR SRS (closes #6523)

comment:8 Changed 18 months ago by Even Rouault

Milestone: 1.11.5

(As far as I can see, the issue existed since support for MERCATOR was introduced in r19742)

comment:9 Changed 18 months ago by Mike Toews

Perfect, thanks! I see the first line was PARAM_1, not PARAM_0.

Note: See TracTickets for help on using tickets.