[PROJ] Pseudo-mercator

[Even Rouault]

Le 01/12/2021 à 22:57, Andrew Patterson a écrit :
> Even,
>
>     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...
>
>
> I should say up front that I have nothing invested in any of these 
> personally -- if they're wrong, that's an answer I can at least relay 
> back. It'll be ugly in some cases because we'll have umpteen hundred 
> maps using a questionable WKT that has to be remedied somehow, but at 
> least I can say there's no code solution at my end.
>
> That said, to play devil's advocate, what would the right way to 
> specify a customized pseudo-mercator look like? It sounds like you're 
> saying that the PROJCS name is the problem? I tried changing the name 
> to 'Unknown' and it had no effect, so either you're talking about a 
> PROJ version later than 8.1.1 or I'm misunderstanding what you mean.

Actually, I missed that latitude_of_center in PROJ >= 6 is completely 
ignored with Mercator_1SP. With Mercator_1SP, the parameter that matters 
is the scale_factor. With Mercator_2SP, the one that matters is 
standard_parallel_1.  So it looks GDAL 2.4 interpreted somehow this 
latitude_of_center as the +lat_ts parameter, but I'd say this was an 
accident. This WKT is definitely weird.


>
> Of the two cases, this is the one that I feel like there should be a 
> way to make it work properly. The next map is sounding more & more 
> like a lousy map to start with. This one is weird, but it should be 
> feasible (if strange) to tweak a pseudo-mercator like this, shouldn't it?
>
>     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.
>
>
> That actually explains something one of our QA told me. He said there 
> was an unnamed conversion in there but I wasn't seeing -- I guess he 
> was using WKT2!

The name of the CONVERSION when exporting as WKT2 is only informational. 
It has no computational consequence. When importing a WKT1 string, PROJ 
will use "unnamed" as the CONVERSION name as the information isn't present.

Even

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211201/c3acc29c/attachment-0001.html>