Axis order wrong in PostGIS 3.1.2 (may relates to #4842)

[tjay]

Since 3.1.2 ST_Transform swaps the coordinates of EPSG:31466 on output. In all cases Proj uses proj-datumgrid-1.8.zip. cs2cs does not show these swaps…

Postgis 3.1.2 and Proj 8.0.1

POSTGIS="3.1.2 cbe925d" [EXTENSION] PGSQL="120" GEOS="3.9.1-CAPI-1.14.2" SFCGAL="1.3.8" PROJ="8.0.1" GDAL="GDAL 3.2.3, released 2021/04/27" LIBXML="2.9.1" LIBJSON="0.11"
TOPOLOGY RASTER

Unexpected result

test=# select ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=25832;POINT(343806.1686 5676555.6896)'),31466),4);
                st_asewkt
-----------------------------------------
 SRID=31466;POINT(5676286.92 2553378.32)
(1 Zeile)

Postgis 3.1.2 and Proj 7.2.1

POSTGIS="3.1.2 cbe925d" [EXTENSION] PGSQL="120" GEOS="3.9.1-CAPI-1.14.2" SFCGAL="1.3.8" PROJ="7.2.1" GDAL="GDAL 3.2.3, released 2021/04/27" LIBXML="2.9.1" LIBJSON="0.11"
TOPOLOGY RASTER

Unexpected result

test=# select ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=25832;POINT(343806.1686 5676555.6896)'),31466),4);
                st_asewkt
-----------------------------------------
 SRID=31466;POINT(5676286.92 2553378.32)
(1 Zeile)

Postgis 3.1.1 and Proj 7.2.1

POSTGIS="3.1.1 aaf4c79" [EXTENSION] PGSQL="120" GEOS="3.9.1-CAPI-1.14.2" SFCGAL="1.3.8" PROJ="7.2.1" GDAL="GDAL 3.2.3, released 2021/04/27" LIBXML="2.9.1"
LIBJSON="0.11" TOPOLOGY RASTER

Expected result

select ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=25832;POINT(343806.1686 5676555.6896)'),31466),4);
                st_asewkt
-----------------------------------------
 SRID=31466;POINT(2553378.32 5676286.92)
(1 Zeile)

[fnadler]

I suspect this behaviour has been introduced by #4842, which states, that in case first column does not include "Lat", the EPSG database order is chosen.

In case of EPSG 31466, proj's internal database (sqlite) outputs the following axis order:

select axis.name, axis.coordinate_system_order
from projected_crs
         inner join axis on projected_crs.coordinate_system_code = axis.coordinate_system_code
where projected_crs.code = 31466
order by coordinate_system_order

--------------------------------------------------------------------------------------------
Northing
Easting

[pramsey]

It sure *looks* like the EPSG definition wants a northing/easting order.

​http://epsg.io/31466

[pramsey]

In fact, the EPSG Remarks say:

"Remarks: Zone width 3 degrees. Also used offshore between 4.5°E and 7.5°E by State Geological Surveys. Not used offshore in oil industry. See CRS code 5676 for variant with axes order reversed to easting before northing for use in GIS applications."

So, use EPSG:5676.

SELECT ST_AsEWKT(
  ST_Transform(
    ST_GeomFromEWKT('SRID=25832;POINT(343806.1686 5676555.6896)'), 
    5676),
  4);

               st_asewkt                
----------------------------------------
 SRID=5676;POINT(2553378.32 5676286.92)

[tjay]

with #4949, the mirrored axis of EPSG:31466 is mirrored back. Is the definition of Proj / EPSG.IO wrong or is PostGIS now implementing the axis orientation incorrectly?

[pramsey]

We are deliberately doing axis order "incorrectly" relative to EPSG guidance, because our user base seems to like it that way. The "desired" orderings are now in regression tests, so the behaviour should stop flip flopping around. For modern versions of proj we are in fact delegating the decision to proj via the proj_normalize_for_visualization() function.