Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability

[Even Rouault]

This ticket has strong connection with #3962, but is focused in the issue of encoding in GeoTIFF.

The way GDAL currently encodes EPSG:3857 in GeoTIFF makes such GeoTIFFs being misplaced in a very significant way when opened in ArcGIS. Basically because ArcGIS must interpret it as Mercator_1SP with ellipsoidal definition of WGS84 instead of the the spherical formula.

GDAL currently writes the following keys:

      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,25): "WGS 84 / Pseudo-Mercator"
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      ProjectedCSTypeGeoKey (Short,1): Unknown-3857
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.

ArcGIS produces EPSG:3857 GeoTIFF with the following keys:

      GTModelTypeGeoKey (Short,1): User-Defined
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,50): "PCS Name = WGS_1984_Web_Mercator_Auxiliary_Sphere"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,13): "GCS_WGS_1984"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257223563    
      PCSCitationGeoKey (Ascii,443): "ESRI PE String = PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Mercator_Auxiliary_Sphere"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["Standard_Parallel_1",0.0],
PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]"
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter

(carriage returned added in ESRI PE string for display)

Since r27594 (GDAL 2.0dev), that brings some support for the "Mercator_Auxiliary_Sphere" projection found in the ESRI PE string in the exportToProj4() code, GDAL is able to correctly reproject such as raster, but reports the ESRI PE as WKT. But using such formulation as general purpose output by GDAL is probably not desirable and wouldn't likely be correctly by other packages (such as older GDAL, and even GDAL 2.0).

Several attempts to find an intermediate way to write ProjectedCSTypeGeoKey=3857 while being correctly registered by ArcGIS have been made, and 2 possible candidates have been found. They basically consist in using ProjectedCSTypeGeoKey=3857, but overriding the GCS related geokeys to use the spheroid instead of the ellipsoid. Which is kind of a hack since technically WebMercator refers to the ellipsoidal WGS84 but with spherical formula.

Formulation 1 :

      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeographicTypeGeoKey (Short,1): User-Defined
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): User-Defined
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6378137          
      GeogPrimeMeridianLongGeoKey (Double,1): 0                
      ProjectedCSTypeGeoKey (Short,1): Unknown-3857
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter

seen in ArcCatalog properties as :

WGS_1984_Pseudo_Mercator
Authority: Custom

Projection: Mercator
false_easting: 0.0
false_northing: 0.0
central_meridian: 0.0
standard_parallel_1: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: 
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_unknown
  Spheroid: Unknown
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6378137.0
    Inverse Flattening: 0.0

or Formulation 2 (the same without GeographicTypeGeoKey):

      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): User-Defined
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6378137          
      GeogPrimeMeridianLongGeoKey (Double,1): 0                
      ProjectedCSTypeGeoKey (Short,1): Unknown-3857
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter

seen in ArcCatalog as :

WGS_1984_Pseudo_Mercator
Authority: Custom

Projection: Mercator
false_easting: 0.0
false_northing: 0.0
central_meridian: 0.0
standard_parallel_1: 0.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: 6378137.0
    Inverse Flattening: 0.0

The difference is that Formulation 2 makes ArcGIS label the datum as WGS_1984, but still using a spherical definition. When using one of those formulations ArcGIS issue a warning about the raster opened being of a different GCS than the data frame it is imported into (even when the data frame uses a WGS84 datum and with formulation 2).

Both are read as the following starting with GDAL 1.9:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

Note the odd-datum named "unkwown", but inside a "WGS 84" GEOGCS with EPSG code 4326e, and above all the ellipsoid defined as a spheroid. When importing one or another of those formulations, we could imagine hacking the GTiff driver to return more traditional WKT - essentially ignoring the spheroidal override on read in this very specific case.

[Even Rouault]

trunk r28923 "GeoTIFF: on reading better deal with a few ESRI formulations of WebMercator (#5924)"

[Even Rouault]

trunk r30435 "morphToESRI(): use Mercator_Auxiliary_Sphere projection for EPSG:3857. morphFromESRI(): map Mercator_Auxiliary_Sphere to EPSG:3857 (#5924)"

trunk r30436 "GTiff: add GEOTIFF_KEYS_FLAVOR=ESRI_PE creation option to write EPSG:3857 in a ESRI compatible way (#5924)"

branches/2.0 r30438 "morphFromESRI(): map Mercator_Auxiliary_Sphere to EPSG:3857 (#5924)"

[Even Rouault]

branches/1.11 r30439 "GTiff: on reading attach PROJ4 extension node for Mercator_Auxiliary_Sphere / WebMercator (#5924)"

branches/1.11 r30440 "morphFromESRI(): map Mercator_Auxiliary_Sphere to EPSG:3857 (#5924)"

[Even Rouault]

trunk r30490, branches/2.0 r30491, branches/1.11 r30492 "Fix recently introduced memory leak in morphFromESRI() with Mercator Auxiliary Sphere (#5924) (CID 1324530)"

[Jukka Rahkonen]

What is missing? Someone to test with recent ArcMap/ArcGIS Pro version, or?

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.