[PROJ] Pseudo-mercator

[Even Rouault]

Andrew,

Le 01/12/2021 à 21:51, Andrew Patterson a écrit :
> Even,
>
> That seems like a logical suggestion!
>
> This is the WKT for the map that's off by ~75km. That 
> 'latitude_of_center' seems to be contributing at the very least, but 
> as I said before, there appears to be more to it than that because 
> without it there's still a ~0.18deg latitude shift.
>
> PROJCS["WGS_84_Pseudo_Mercator",GEOGCS["WGS 
> 84",DATUM["wgs_1984",SPHEROID["WGS 
> 1984",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["false_easting",0],PARAMETER["false_northing",0],PARAMETER["latitude_of_center",6.64456744726493],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],UNIT["Meter",1,AUTHORITY["EPSG","9001"]]]
>
That one is a very odd beast. The 
PARAMETER["latitude_of_center",6.64456744726493] is completely 
inconsistant with the name "WGS_84_Pseudo_Mercator". I see GDAL 2.4 
honoured this latitude of center, as gdalsrsinfo on it returns 
+proj=merc +lon_0=0 +lat_ts=6.64456744726493 +x_0=0 +y_0=0 +ellps=WGS84 
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs.    Newer PROJ versions have a 
heustistics that only look at the name WGS_84_Pseudo_Mercator and ignore 
the projection parameters to decide to use the official EPSG:3857 "WGS 
84 / Pseudo-Mercator" definition. One could argue it is a bug to ignore 
the latitude_of_center parameter, but this definition is very confusing 
to start with. It is hard to know what the intention of the user was: 
use official EPSG:3857 "WGS 84 / Pseudo-Mercator", or use a modified 
version of it...

> This is the WKT for the map that's off by ~26km:
>
> PROJCS["Popular_Visualisation_CRS_Mercator_deprecated",GEOGCS["GCS_Popular 
> Visualisation 
> CRS",DATUM["D_Popular_Visualisation_Datum",SPHEROID["Popular_Visualisation_Sphere",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],PARAMETER["standard_parallel_1",0],UNIT["Meter",1]]
>
> Looking at it more closely, I'm not sure it actually is a 
> pseudo-mercator  map. I think I say 'populator visualisation' and 
> assumed it was. But it _is_ off by ~26km after moving from 2.4/6/1 -> 
> 3.3/8.1.1

As far as I can see, GDAL 2.4 can't transform that to a PROJ.4 string at 
all (it doesn't understand the "Mercator" projection), so I'm not sure 
why you get a result at all with it. With recent PROJ,

$ projinfo 
'PROJCS["Popular_Visualisation_CRS_Mercator_deprecated",GEOGCS["GCS_Popular 
Visualisation 
CRS",DATUM["D_Popular_Visualisation_Datum",SPHEROID["Popular_Visualisation_Sphere",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],PARAMETER["standard_parallel_1",0],UNIT["Meter",1]]'

returns

PROJ.4 string:
+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs 
+type=crs

WKT2:2019 string:
PROJCRS["Popular_Visualisation_CRS_Mercator_deprecated",
     BASEGEOGCRS["GCS_Popular Visualisation CRS",
         DATUM["Popular Visualisation Datum",
             ELLIPSOID["Popular_Visualisation_Sphere",6378137,0,
                 LENGTHUNIT["metre",1]],
             ID["EPSG",6055]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["Degree",0.0174532925199433]]],
     CONVERSION["unnamed",
         METHOD["Mercator (variant B)",
             ID["EPSG",9805]],
         PARAMETER["Latitude of 1st standard parallel",0,
             ANGLEUNIT["Degree",0.0174532925199433],
             ID["EPSG",8823]],
         PARAMETER["Longitude of natural origin",0,
             ANGLEUNIT["Degree",0.0174532925199433],
             ID["EPSG",8802]],
         PARAMETER["False easting",0,
             LENGTHUNIT["metre",1],
             ID["EPSG",8806]],
         PARAMETER["False northing",0,
             LENGTHUNIT["metre",1],
             ID["EPSG",8807]]],
     CS[Cartesian,2],
         AXIS["(E)",east,
             ORDER[1],
             LENGTHUNIT["metre",1,
                 ID["EPSG",9001]]],
         AXIS["(N)",north,
             ORDER[2],
             LENGTHUNIT["metre",1,
                 ID["EPSG",9001]]]]

Which is a reasonable enough interpretation of this WKT. It could have 
potentially been better recognized as the deprecated EPSG:3785 (which 
PROJ then interprets as EPSG;3857), but there's no heuristics for that 
particular formulation, so it is recognized as classic Mercator on a 
sphere.  In practice, it should give identical results to official 
EPSG:3857 "WGS 84 / Pseudo-Mercator", except when a datum transformation 
is involved. So if you transform between WGS 84 and this, you should get 
correct results. But results off by 26 km are beyond an ignored datum 
shift (differences should be off by a few hundred of meters max), so I'm 
not sure why you observe this.

>
> This is the WKT for a map that saw no change between GDAL/PROJ versions:
>
> 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["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc 
> +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m 
> +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

Yes that one is the canonical WKT1:GDAL formulation

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.