[postgis-users] PostGIS & PROJ 6 -- Y axis shift with 900913

[Shira Bezalel]

On Mon, Jul 20, 2020 at 11:25 AM Shira Bezalel <shira at sfei.org> wrote:

> We recently recompiled PostGIS in one of our Postgres 12.3 environments
> with PROJ 6 (6.3.2) from PROJ 4 (4.9.3).
>
> postgis_full_version() output:
> POSTGIS="3.0.1 ec2a9aa" [EXTENSION] PGSQL="120" GEOS="3.7.1-CAPI-1.11.1
> 27a5e771" PROJ="6.3.2" LIBXML="2.9.4" LIBJSON="0.12.1"
>
> We're now noticing that using ST_Transform(geometry, srid) to transform
> some points in California from SRID 4326 to 900913 is causing a shift along
> the Y axis of about 20-30 meters. The difference is evidenced by a
> comparison of the old and new ST_Y() values.
>
> If I force the PROJ 4 definition to be used, the transformation works fine
> with no shift. (To force the PROJ 4 definition, I used the method described
> in  http://blog.cleverelephant.ca/2019/02/proj4-postgis.html and set
> authname, authsrid and srtext to null so that the proj4text column is used.)
>
> Is this shift expected with PROJ 6? I realize 900913 is deprecated and we
> plan to migrate to 3857 as a longer term solution, but looking for
> a shorter term fix if possible since much data and code is impacted by
> this. Note we are seeing no shift issues with SRID 3857.
>
> Let me know if I can provide more information.
>
> Thanks much,
> Shira
>
>
Hello again. Wanted to follow up on this. We have found a workaround that
seems somewhat better than just setting the auth_name, auth_srid and srtext
columns to null especially since external applications rely on the srtext
column. If we replace the record for srid 900913 in the spatial_ref_sys
table with the one at epsg.io (https://epsg.io/900913) then the points are
correctly transforming without a shift with st_transform(geometry,900913).
As a related consequence of this, I am noticing a shift when exporting
these 900913 points to shapefile format with pgsql2shp (it's fine with the
old srtext definition), but we don't use pgsql2shp much so I think that's
okay. Our long-term plan is to move away from 900913, but this might be
suitable for the short-term. If I'm missing any pitfalls with this approach
(apart from the pgsql2shp one), please feel free to let me know.

New srtext value from epsg.io:

PROJCS["Google Maps Global 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.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_2SP"],PARAMETER["standard_parallel_1",0],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],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","900913"]]

Original srtext value in spatial_ref_sys:

PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular
Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular
Visualisation
Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]

Thanks,
Shira
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20200722/0f93737d/attachment.html>