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

[Dmitry Baryshnikov]

15.04.2015 23:46, Even Rouault пишет:
> Le mercredi 15 avril 2015 22:14:37, Dmitry Baryshnikov a écrit :
>> As for me the main problem is, that gdal write wrong WKT for 3857
>> (http://trac.osgeo.org/gdal/ticket/3962).
>> The Geographic SRS use WGS84 ellipsoid
>> (SPHEROID["WGS_1984",6378137,298.257223563]]) instead sphere
>> (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"],) and all the data shifts.
> I do think that the above SPHEROID WKT and Proj.4 are both OK, even if they
> don't look consistent.
The ArcGIS (or QGIS in qpj) save 3857 in such prj file:

PROJCS["WGS 84 / Pseudo-Mercator",
     GEOGCS["WGS 84",
         DATUM["WGS_1984",
             SPHEROID["WGS 84",6378137,298.257223563,
                 AUTHORITY["EPSG","7030"]],
             AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
             AUTHORITY["EPSG","9122"]],
         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"]],
     AXIS["X",EAST],
     AXIS["Y",NORTH],
     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"]]

And GDAL opens such shape file and set the spheroid.
But if I save the shape file in 3857 via GDAL I get ellipsoid:

   PROJCS["WGS_84_Pseudo_Mercator",
     GEOGCS["GCS_WGS_1984",
         DATUM["D_WGS_1984",
             SPHEROID["WGS_1984",6378137,298.257223563]],
         PRIMEM["Greenwich",0],
         UNIT["Degree",0.017453292519943295]],
         PROJECTION["Mercator"],
         PARAMETER["central_meridian",0],
         PARAMETER["false_easting",0],
         PARAMETER["false_northing",0],
         UNIT["Meter",1],
         PARAMETER["standard_parallel_1",0.0]]

And both GDAL(QGIS) and ArcGIS open this shape file and set the ellipsoid!

> The datum is technically WGS84 datum (ellipsoidal), even if the projection
> only uses the semi major axis (spherical). This is all the hell of this
> WebMercator projection ! We use a hack of proj.4 to use classical mercator
> transformation with spherical parameters, while still being considered as a
> WGS84 datum. When reprojecting from a non-WGS 84 datum (let's say EPSG:4322)
> to EPSG:3857, the datum transformation should be from this non-WGS 84 to
> ellipsoid WGS 84, and then the WebMercator projection uses the spherical
> parameters to compute the easting/northing.
>
> Look:
>
> 1) Normal EPSG:3857 transformation (with +nadgrids=@null) :
>
> $ echo "2 49" | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:3857
> 222656.112419298 6274866.20215882 2.9538588412106
>
> --> correct. Same result as operating in 2 steps :
>
> $ echo "2 49" | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:4326 |
> gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
> 222656.112419298 6274866.20215883 2.9538588412106
>
> 2) Incorrect definition using +towgs84=0,0,0,0,0,0,0 :
>
> $ gdaltransform -s_srs EPSG:4322 -t_srs "+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 +towgs84=0,0,0,0,0,0,0"
> 2 49
> 222656.112419297 6242580.32725437 -12133.4419494262
>
> 2 49 is transformed from EPSG:4322 datum to a spherical datum, instead to the
> ellipsoid WGS84 datum.
>
> 3) Incorrect definition without +towgs84 or +nadgrids=@null : no datum
> transformation at all is done since its only an ellipsoid definition, and
> proj.4 will not do any ellipsoid transformation in that case
>
> $ gdaltransform -s_srs EPSG:4322 -t_srs "+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"
> 2 49
> 222638.981586547 6274861.39400658 0
>
>
>> Best regards,
>>       Dmitry
>>
>> 15.04.2015 19:15, Even Rouault пишет:
>>> Hi,
>>>
>>> I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of
>>> issues and findings when trying to write GeoTIFF files in EPSG:3857 in a
>>> "standard" way (using ProjectedCSTypeGeoKey = 3857), while making them
>>> correctly displayed by ArcGIS (and ideally by other independant
>>> implementations).
>>> The current situation is that there's no such way (AFAIK). I'd appreciate
>>> any feedback/suggestion related to that issue.
>>>
>>> Best regards,
>>>
>>> Even
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev

Best regards,
     Dmitry