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

[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 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

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com