[gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability

[Andre Vautour]

On 2015-04-15 14:59, Even Rouault wrote:
> Le mercredi 15 avril 2015 18:21:51, xavier lhomme a écrit :
>> Dis you try with 102113 or 102100  instead of 3857 ? I had once to use
>> those code in order to be discovered as 3857 in arcgis. when i created
>> geotiff with gdal warp/transform.
> Xavier,
>
> thanks for the hint. I've just tried it and thinks it qualifies as yet another
> not completely satisfactory solution.
>
> Using -a_srs EPSG:102113 effectively results in a correctly placed geotiff in
> ArcGIS (with a warning about it being in GCS_WGS_1984_Major_Auxiliary_Sphere),
> but there are several points that are not that great :
> - this is not a standard EPSG code
> - due to being > 32767, it is not encoded in ProjectedCSTypeGeoKey, which was
> one of my "requirement"
> - when read by GDAL, the expanded GeoTIFF definition isn't properly identified
> as WebMercator, consequently the proj.4 string lacks the "+nadgrids=@null"
> stuff, which makes it a non-complete datum definition and cause issues when
> reprojecting such GeoTIFF with GDAL (lack of datum transformation from/to the
> WGS84 datum of WebMercator). That could admitedly being tweaked in GDAL.
>
> For the record, the GeoTIFF keys written by GDAL are:
> {{{
>        GTModelTypeGeoKey (Short,1): ModelTypeProjected
>        GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
>        GTCitationGeoKey (Ascii,22): "WGS_1984_Web_Mercator"
>        GeographicTypeGeoKey (Short,1): User-Defined
>        GeogCitationGeoKey (Ascii,151): "GCS Name =
> GCS_WGS_1984_Major_Auxiliary_Sphere|Datum = WGS_1984_Major_Auxiliary_Sphere|
> Ellipsoid = WGS_1984_Major_Auxiliary_Sphere|Primem = Greenwich|"
>        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): User-Defined
>        ProjectionGeoKey (Short,1): User-Defined
>        ProjCoordTransGeoKey (Short,1): CT_Mercator
>        ProjLinearUnitsGeoKey (Short,1): Linear_Meter
>        ProjNatOriginLongGeoKey (Double,1): 0
>        ProjNatOriginLatGeoKey (Double,1): 0
>        ProjFalseEastingGeoKey (Double,1): 0
>        ProjFalseNorthingGeoKey (Double,1): 0
>        ProjScaleAtNatOriginGeoKey (Double,1): 1
> }}}
>
> And GDAL reads that as :
> {{{
> PROJCS["WGS_1984_Web_Mercator",
>      GEOGCS["GCS_WGS_1984_Major_Auxiliary_Sphere",
>          DATUM["WGS_1984_Major_Auxiliary_Sphere",
>              SPHEROID["WGS_1984_Major_Auxiliary_Sphere",6378137,0]],
>          PRIMEM["Greenwich",0],
>          UNIT["degree",0.0174532925199433]],
>      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"]]]
> }}}
> (not sure why it lacks the latitude_of_origin,  but that's another minor
> issue)
>
> ArcCatalog identifies the SRS as:
> {{{
> WGS_1984_Web_Mercator
> WKID: 3785 Authority: EPSG
>
> 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_Major_Auxiliary_Sphere
> Angular Unit: Degree (0.0174532925199433)
> Prime Meridian: Greenwich (0.0)
> Datum: D_WGS_1984_Major_Auxiliary_Sphere
>    Spheroid: WGS_1984_Major_Auxiliary_Sphere
>      Semimajor Axis: 6378137.0
>      Semiminor Axis: 6378137.0
>      Inverse Flattening: 0.0
> }}}
>
> Even

The way I see it, there are two different ways to model Google Mercator:
1. As a WGS84 datum/ellipsoid with a custom Mercator-based projection 
which uses only the semi-major axis, as seems to be currently done in 
EPSG: http://epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857
2. As the Google spheroid with Mercator 1SP and no transformation when 
going to WGS 84.

I don't think #1 is currently an option given the list of coordinate 
transformation codes available in the GeoTIFF specification.

For #2, I would think/hope that if you specified the identity 
7-parameter transformation in the GeoTIFF header, ArcGIS would properly 
reproject the image:
http://trac.osgeo.org/geotiff/ticket/26

André