Opened 2 years ago

Last modified 23 months ago

#5924 new defect

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

Reported by: Even Rouault Owned by: warmerdam
Priority: normal Milestone:
Component: GDAL_Raster Version: unspecified
Severity: normal Keywords: gtiff 3857 webmercator
Cc:

Description

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.

Change History (4)

comment:1 Changed 2 years ago by Even Rouault

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

comment:2 Changed 23 months ago by 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)"

comment:3 Changed 23 months ago by 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)"

comment:4 Changed 23 months ago by 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)"

Note: See TracTickets for help on using tickets.