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@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Best regards,
    Dmitry

_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to