Transform regression for N,E CRS

[rcoup]

#4842 (in v3.1.2) fixed a bug in the behaviour change made in #4748 (v3.0.3), but has introduced a wider regression.

in v3.1.1 and earlier (expected):

# SELECT st_asewkt(st_transform('SRID=2193;POINT(1680238 5863854)'::geometry, 4326));
SRID=4326;POINT(174.863597538742 -36.785298415230315)

v3.1.2:

# select st_asewkt(st_transform('SRID=2193;POINT(1766289 5927325)'::geometry, 4326));
SRID=4326;POINT(-117.47905633499929 -50.9729929889798)

Adding ST_FlipCoordinates() returns the right output, so it's definitely a CRS Axis Order issue:

# select st_asewkt(st_transform(st_flipcoordinates('SRID=2193;POINT(1766289 5927325)'::geometry), 4326));
SRID=4326;POINT(174.863597538742 -36.785298415230315)

NZ Transverse Mercator (NZTM) (​EPSG:2193) is defined with N, E axis order due to some quirk of history, but it's the main CRS used in New Zealand mapping, not an obscure oddball one.

PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980", 6378137, 298.257222101, AUTHORITY["EPSG", "7019"]],
            TOWGS84[0, 0, 0, 0, 0, 0, 0],
            AUTHORITY["EPSG", "6167"]],
        PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
        UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]],
        AUTHORITY["EPSG", "4167"]],
    UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin", 0],
    PARAMETER["central_meridian", 173],
    PARAMETER["scale_factor", 0.9996],
    PARAMETER["false_easting", 1600000],
    PARAMETER["false_northing", 10000000],
    AUTHORITY["EPSG", "2193"],
    AXIS["Northing", NORTH],
    AXIS["Easting", EAST]]

[rcoup]

Ah, managed to mess up my copying and pasting. Should be:

---

in v3.1.1 and earlier (expected):

# SELECT st_asewkt(st_transform('SRID=2193;POINT(1766289 5927325)'::geometry, 4326));
SRID=4326;POINT(174.863597538742 -36.785298415230315)

[pramsey]

I can't tell whether your example point is in (northing, easting) order...

'SRID=2193;POINT(1766289 5927325)'

The code as it exists is pretty simple: if the CRS is in geographics (axis 1 = Latitude) then swap coordinates, otherwise report/ingest coordinates in their EPSG order. For EPSG:2193, that order, in the EPSG database, is northing,easting.

Extracting the problem from PostGIS altogether.

cs2cs -v
Rel. 6.3.3, July 1st, 2020

echo -36.785298415230315 174.863597538742 | cs2cs EPSG:4326 EPSG:2193 

5927325.00	1766289.00

So, is 1766289 the easting, or the northing? EPSG says that (a) it's the easting and (b) it should come *second* in the axis order. Proj does that behaviour, and as far as I can tell so does the current state of postgis:

select st_asewkt(st_transform('SRID=4326;POINT(174.863597538742 -36.785298415230315)'::geometry, 2193));

                  st_asewkt                  
---------------------------------------------
 SRID=2193;POINT(5927325 1766289.0000000002)

}}}

[rcoup]

The (1766289 5927325) coordinate is in Easting,Northing.

AFAIK X/Lon, Y/Lat has always been postgis' axis ordering? I'm assuming that assumption hasn't changed - certainly isn't noted in the release notes, and seems a fairly critical change for a patch bump...

I know Proj' behaviour has changed with v6 but I assumed PostGIS would be keeping it's existing behaviour (roughly ala ​GDAL's "traditional GIS" axis ordering)

[pramsey]

I see that under the old Proj 4 regime the output was easting/northing.

echo 174.863597538742 -36.785298415230315 | ./cs2cs +init=EPSG:4326 +to +init=EPSG:2193

1766289.00	5927325.00

So my guess is correct, the "old way", pre-proj6 and the internal EPSG database, was to output EPSG:2193 in easting/northing order, notwithstanding the EPSG definition.

[pramsey]

For history tracking, note that this issue is bound up with #4748, #4844, and #4842.

[rcoup]

PostGIS/etc has always done E/N with 2193 (& a group of associated local circuit CRS) since I started working with it >15 years ago.

With CRS WKT1 the axes were never really were visible (even now SELECT srtext FROM spatial_ref_sys WHERE srid=2193; doesn't show it), everything was always treated as X/E/Lng, Y/N/Lat. Same as EPSG:4326 ;-) And proj4 didn't do axis ordering AFAIK.

Maybe becomes a larger question (PostGIS 4?):

1. Should postgis support axis ordering ala the CRS defintion? Fundamentally, probably yes.
2. Should postgis continue to store data in XYZM order? In which case, seems like axis ordering is an operation to apply only at input/output/transform time.
3. Or should it be storing data in the CRS order? In which case, seems like many beahviour/definitions need to change wrt functions? eg: ST_XMax() returns "first axis" not X/E/Lng. And the IO-time/SetSRID operations need to change to flip as appropriate.
4. What's the path for a user to migrate from what they have now? Options at a function/geometry/column/user/session/database level?

WRT a simpler fix for the regression though...

Should ST_Transform() be calling proj_normalize_for_visualization()(​ref)? Then we can say everything is X/Lng/E, Y/Lat/N ala all other PostGIS versions, and continue to ignore the CRS axis ordering.

[pramsey]

I've got a fix cadged from GDAL that seems to do everything you want while not breaking the polar projections. I'm not sure what version of proj proj_normalize_for_visualization() showed up in, but I think it's a little recent, though it might be a nice improvement to make? Trouble is one of consistency over time (as we're demonstrating here).

[Paul Ramsey <pramsey@…>]

In 6f37390/git:

Flip N/E systems to E/N and geodetic systems to Lon/Lat, while leaving Polar systems as-is, references #4949

[Paul Ramsey <pramsey@…>]

In 8baf0b0/git:

Flip N/E systems to E/N and geodetic systems to Lon/Lat, while leaving Polar systems as-is, references #4949, 3.1 branch

[Paul Ramsey <pramsey@…>]

In 6a831ee1/git:

Flip N/E systems to E/N and geodetic systems to Lon/Lat, while leaving Polar systems as-is, references #4949, 3.0 branch

[rcoup]

Awesome, thanks for the fix.

According to ​https://github.com/OSGeo/PROJ/commit/6a7e24dce79f93b73f4919f267df2fdf3ee95713 proj_normalize_for_visualization() was introduced in Proj v6.1.0. Ironically ​the issue references postgis proj_crs_is_swapped() as the inspiration ;-)

[Mike Taves]

Yes, ironic to see the context! Agree that proj_normalize_for_visualization() is a more sustainable method for POSTGIS_PROJ_VERSION > 61. Pyproj's ​Transformer class has an always_xy option that leverages this built-in functionality.

See more here: ​https://proj.org/development/quickstart.html

[robe]

This was miscategorized