ESRI morphing broken for Mercator with non-zero standard_parallel_1

[mikrit]

To reproduce, take the ESRI WKT for EPSG:3001 from ​http://help.arcgis.com/en/arcims/10.0/mainhelp/mergedProjects/ArcXMLGuide/elements/pcs.htm which (when pretty-printed) is

PROJCS["Batavia_NEIEZ",
    GEOGCS["GCS_Batavia",
        DATUM["D_Batavia",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Mercator"],
    PARAMETER["False_Easting",3900000.0],
    PARAMETER["False_Northing",900000.0],
    PARAMETER["Central_Meridian",110.0],
    PARAMETER["Standard_Parallel_1",4.45405154589751],
    UNIT["Meter",1.0]]

Save it in a file called "Batavia_NEIEZ.prj" and do

gdalsrsinfo -p "Batavia_NEIEZ.prj"

which gives

PROJ.4 : '+proj=merc +lon_0=110 +lat_ts=4.45405154589751 +x_0=3900000 +y_0=90000
0 +ellps=bessel +units=m +no_defs '

OGC WKT :
PROJCS["Batavia_NEIEZ",
    GEOGCS["GCS_Batavia",
        DATUM["Batavia",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["False_Easting",3900000.0],
    PARAMETER["False_Northing",900000.0],
    PARAMETER["Central_Meridian",110.0],
    PARAMETER["latitude_of_origin",4.45405154589751],
    UNIT["Meter",1.0]]

Note that the standard_parallel_1 has been renamed to latitude_of_origin. That's not right.

• A standard parallel at 4.454... degrees North means that the scale factor has been adjusted so that the scale is true at that parallel. But the projected point (False_Easting, False_Northing) should still coincide with the geographic point (Central_Meridian, Latitude_of_Origin), which would be (110dE, 0dN) here since the Latitude_of_Origin of an ESRI Mercator definition is always zero, at least this was true in 2008: ​http://lists.osgeo.org/pipermail/metacrs/2008-November/000178.html

• On the other hand, a latitude_of_origin at 4.454... degrees North without a specified standard_parallel_1 or scale_factor, should mean that the scale is true on the equator, but (False_Easting, False_Northing) should coincide with (Central_Meridian, latitude_of_origin) which would be (110dE, 4.454...dN) here. That's different.

The erroneous renaming of standard_parallel_1 to latitude_of_origin was introduced in http://trac.osgeo.org/gdal/changeset/7542 . Some consolation is that the renaming doesn't happen when the projection name in ESRI WKT is given as "Mercator_2SP" or "Mercator (2SP)", which perhaps would be a better projection name. But the ArcGis page where I got the ESRI WKT uses just "Mercator".

An analogous ticket for Equidistant Cylindrical is http://trac.osgeo.org/gdal/ticket/3036 .

I can add that the confusion, at least for Equidistant Cylindrical, may be caused by old and bad EPSG terminology, where the Equidistant Cylindrical used to be Operation Method 9823 and have a "Latitude of Natural Origin" parameter that was really intended to be the standard parallel. This was changed as part of EPSG Change Request 2009.023, in that Operation Method 9823 was replaced by Operation Method 1028 where the parameter has the better name "Latitude of 1st standard parallel".

[avautour]

I ran into this as well.

A couple of clarifications/observations/suggestions:

1. The latitude of natural origin for Mercator should always be zero (Equator) by definition, so having OGRSpatialReference::SetMercator() take in the center latitude is problematic. ESRI, EPSG and Proj4 are all consistent on that point as far as I can tell.
2. ESRI's Mercator is really Mercator 2SP, so it stands to reason that SRS_PT_MERCATOR_2SP should be mapped to it by simply renaming the projection's name.
3. The formulas for converting between the scale factor of 1SP to the first standard parallel of 2SP are as follows for the Mercator projection. See the EPSG Guidance Notes at ​http://www.iogp.org/pubs/373-07-2.pdf:

k0

scale factor at natural origin

sp

first standard parallel

e

ellipsoid's first eccentricity

2SP to 1SP

k0 = cos(sp) / sqrt( 1 - (e² * sin²(sp)) )

1SP to 2SP

sp = acos( sqrt( (1 - e²) / (1 / k0²) - e²) ) )

[avautour]

Rough math proof to convert from scale factor at natural origin to first standard parallel

[Even Rouault]

Sorry for my dubness, but on the example of that ticket what should be the corresponding OGC WKT: just changing "Mercator" to "Mercator_2SP" ? What would be the 2nd parallel: 0 ?

[avautour]

No worries, the 2SP convention confused me as well at the beginning.

For Mercator 2SP, both standard parallels are the same latitude in the different hemispheres. Since you can infer the second standard parallel simply by negating the first, it is usually implicitly defined and not specified as a parameter.

So, to answer your question, yes, I believe going from ESRI WKT to OGC WKT, you'd simply need to rename "Mercator" to SRS_PT_MERCATOR_2SP.

[Even Rouault]

Ah you're right now I see ​http://www.remotesensing.org/geotiff/proj_list/mercator_2sp.html.

Hopefully fixed by r29414 "OSR: correctly transform Mercator_2SP and _1SP to ESRI Mercator, and back from ESRI Mercator to Mercator_2SP (#4861)"

Thanks for the maths. I've checked on one sample with gdaltransform that the 1SP and 2SP variants were equivalent.

[Even Rouault]

trunk r29415 "Fix osr_esri_28 test to work on all platforms (#4861)"

[Even Rouault]

Backported in branches/2.0 in r29421

[Even Rouault]

In 34494:

OSR morphToESRI(): correctly compute standard_parallel_1 of Mercator(2SP) projection from scale factor of Mercator(1SP) (fixes #6456, #4861)

[Even Rouault]

In 34495:

OSR morphToESRI(): correctly compute standard_parallel_1 of Mercator(2SP) projection from scale factor of Mercator(1SP) (fixes #6456, #4861)

[Even Rouault]

In 34496:

OSR morphToESRI(): correctly compute standard_parallel_1 of Mercator(2SP) projection from scale factor of Mercator(1SP) (fixes #6456, #4861)

[Even Rouault]

Note: fix in 2.0 branch *post* 2.0.3RC1

[Even Rouault]

A 2.0.3RC2 has been pushed, so this will be in final 2.0.3