Mercator 2SP modified to Mercator 1SP by GTiff

[Mike Taves]

Mercator (2SP) projections ​EPSG:3388 and ​EPSG:5641 are not properly stored by the GTiff driver, having different PROJ.4 parameters. The outputs are missing the +lat_ts/standard_parallel_1 data, and insert a +k=1 scale factor. The WKT shows that the projection was modified from Mercator_2SP to Mercator_1SP.

See test function test_prj_roundtrip

import difflib
from osgeo import gdal
from osgeo import osr
gdal.UseExceptions()
osr.UseExceptions()

def test_prj_roundtrip(code, show_proj4_diff=True, show_wkt_diff=False):
    ret = 'EPSG:' + str(code) + ' '
    srs = osr.SpatialReference()
    try:
        srs.ImportFromEPSG(code)
    except RuntimeError as e:
        return ret + str(e).replace('\n', ' ')
    srs_wkt = srs.ExportToWkt()
    try:
        srs_proj4 = srs.ExportToProj4()
    except RuntimeError as e:
        return ret + str(e).replace('\n', ' ')
        if reason in proj4_errors:
            proj4_errors[reason].add(code)
        else:
            proj4_errors[reason] = set([code])
        srs_proj4 = None
    fname = '/vsimem/tmp' + str(code)
    drv = gdal.GetDriverByName('GTiff')
    ds = drv.Create(fname, 2, 3)
    assert ds.SetProjection(srs_wkt) == 0
    # Close and re-open
    ds = None
    ds = gdal.Open(fname)
    srs_out = osr.SpatialReference()
    srs_out_wkt = ds.GetProjection()
    ds = None
    gdal.Unlink(fname)
    if srs_out:
        try:
            assert srs_out.ImportFromWkt(srs_out_wkt) == 0
        except RuntimeError as e:
            return ret + (str(e).replace('\n', ' ') + ' with: ' + srs_out_wkt)
    else:
        return ret + 'output does not have any projection info'

    # Show differences
    if show_proj4_diff:
        try:
            srs_out_proj4 = srs_out.ExportToProj4()
        except RuntimeError as e:
            return ret + str(e).replace('\n', ' ')
        if not srs_proj4 and not srs_out_proj4:
            ret += 'no PROJ.4 to compare'
        elif not srs_proj4 and srs_out_proj4:
            ret += 'only output has PROJ.4 : ' + srs_out_proj4
        elif srs_proj4 and not srs_out_proj4:
            ret += 'only source has PROJ.4 : ' + srs_proj4
        elif srs_proj4 != srs_out_proj4:
            ret += '\nDifferences in PROJ.4:\n'
            ret += ''.join(difflib.ndiff(
                [srs_proj4 + '\n'],
                [srs_out_proj4 + '\n'])).rstrip()
    if show_wkt_diff:
        try:
            srs_wkt = srs.ExportToPrettyWkt()
        except RuntimeError as e:
            return ret + str(e).replace('\n', ' ')
        srs_out_wkt = srs_out.ExportToPrettyWkt()
        if srs_wkt != srs_out_wkt:
            ret += ('\nDifferences in WKT:\n')
            ret += ''.join(difflib.ndiff(
                srs_wkt.splitlines(True),
                srs_out_wkt.splitlines(True))).rstrip()
    srs_out = None
    return ret

And test code:

for code in [3388, 5641]:
    print(test_prj_roundtrip(code, True, True))

has output:

EPSG:3388 
Differences in PROJ.4:
- +proj=merc +lon_0=51 +lat_ts=42 +x_0=0 +y_0=0 +ellps=krass +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12 +units=m +no_defs 
?                       ^^^^^^ ^^
+ +proj=merc +lon_0=51 +k=1 +x_0=0 +y_0=0 +ellps=krass +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12 +units=m +no_defs 
?                       ^ ^
Differences in WKT:
  PROJCS["Pulkovo 1942 / Caspian Sea Mercator",
      GEOGCS["Pulkovo 1942",
          DATUM["Pulkovo_1942",
-             SPHEROID["Krassowsky 1940",6378245,298.3,
?                                                    ^
+             SPHEROID["Krassowsky 1940",6378245,298.2999999999985,
?                                                    ^^^^^^^^^^^^^
                  AUTHORITY["EPSG","7024"]],
              TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
              AUTHORITY["EPSG","6284"]],
-         PRIMEM["Greenwich",0,
+         PRIMEM["Greenwich",0],
?                             +
-             AUTHORITY["EPSG","8901"]],
-         UNIT["degree",0.0174532925199433,
+         UNIT["degree",0.0174532925199433],
?                                         +
-             AUTHORITY["EPSG","9122"]],
          AUTHORITY["EPSG","4284"]],
-     PROJECTION["Mercator_2SP"],
?                          ^
+     PROJECTION["Mercator_1SP"],
?                          ^
-     PARAMETER["standard_parallel_1",42],
      PARAMETER["central_meridian",51],
+     PARAMETER["scale_factor",1],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
      AUTHORITY["EPSG","3388"]]
EPSG:5641 
Differences in PROJ.4:
- +proj=merc +lon_0=-43 +lat_ts=-2 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
?                        ^^^^^^ ^^
+ +proj=merc +lon_0=-43 +k=1 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
?                        ^ ^
Differences in WKT:
  PROJCS["SIRGAS 2000 / Brazil Mercator",
      GEOGCS["SIRGAS 2000",
          DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
-             SPHEROID["GRS 1980",6378137,298.257222101,
+             SPHEROID["GRS 1980",6378137,298.2572221010002,
?                                                      ++++
                  AUTHORITY["EPSG","7019"]],
              TOWGS84[0,0,0,0,0,0,0],
              AUTHORITY["EPSG","6674"]],
-         PRIMEM["Greenwich",0,
+         PRIMEM["Greenwich",0],
?                             +
-             AUTHORITY["EPSG","8901"]],
-         UNIT["degree",0.0174532925199433,
+         UNIT["degree",0.0174532925199433],
?                                         +
-             AUTHORITY["EPSG","9122"]],
          AUTHORITY["EPSG","4674"]],
-     PROJECTION["Mercator_2SP"],
?                          ^
+     PROJECTION["Mercator_1SP"],
?                          ^
-     PARAMETER["standard_parallel_1",-2],
      PARAMETER["central_meridian",-43],
+     PARAMETER["scale_factor",1],
      PARAMETER["false_easting",5000000],
      PARAMETER["false_northing",10000000],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
-     AXIS["X",EAST],
-     AXIS["Y",NORTH],
      AUTHORITY["EPSG","5641"]]

[Even Rouault]

Fixed in upstream libgeotiff / trunk internal libgeotiff.

[Even Rouault]

Milestone renamed