inaccurate rect reprojection/transformation for 28992 to 4326

[roelarents]

Expected behavior and actual behavior.

Expected the only feature in the dataset in the output, instead of an empty result set.

Steps to reproduce the problem.

Have a source in NetherlandsRD (EPSG:28992) with a single feature/geom.

Do not set an EXTENT in the WFS config in the mapfile.

Make a WFS request with srsName EPSG:4326.

QUERY_STRING=service=WFS&request=GetFeature&count=1&version=2.0.0&outputFormat=geojson&typeName=test&srsName=EPSG:4326

Analysis

(Disclaimer: I'm not very deep into MapServer and C++, but this is what I could come up with while debugging a bit:)

The rectangle from the source is taken and reprojected to the output SRS. When constructing a (proj) transformation for that, +epsgaxis=ne is added to the proj string for the output, even though EPSG:4326 already has north/easting cq lat/lon cq y/x axis. Somewhere along the line (I forgot the breakpoint, sorry), the WKT for EPSG:4326 then becomes some "unknown" thing without authority name and code. And the resulting transformation from proj does not use RDNAPTRANS but some ballparking thing (I think, I found it hard to debug the actual resulting projection operations).

Commenting out the line that adds +epsgaxis=ne line does result in one feature (with lon/lat order coordinates).

The unknown thing:

"GEOGCRS[\"unknown\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]"

vs the real thing:

"GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\",ID[\"EPSG\",1166]],MEMBER[\"World Geodetic System 1984 (G730)\",ID[\"EPSG\",1152]],MEMBER[\"World Geodetic System 1984 (G873)\",ID[\"EPSG\",1153]],MEMBER[\"World Geodetic System 1984 (G1150)\",ID[\"EPSG\",1154]],MEMBER[\"World Geodetic System 1984 (G1674)\",ID[\"EPSG\",1155]],MEMBER[\"World Geodetic System 1984 (G1762)\",ID[\"EPSG\",1156]],MEMBER[\"World Geodetic System 1984 (G2139)\",ID[\"EPSG\",1309]],ELLIPSOID[\"WGS84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7030]],ENSEMBLEACCURACY[2.0],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],USAGE[SCOPE[\"unknown\"],AREA[\"World.\"],BBOX[-90,-180,90,180]]]"

Related

• no crs output in geojson when reprojecting to epsg with inverted axes #6966 -> also due to +epsgaxis=ne (I think)
• rectObj wrong projection with PROJ > 6 #6478 -> also a rect reprojection going wrong

Attached simple test case

wfs-test.zip

# commands.txt
mapserv -nh "QUERY_STRING=map=wfs-test.map&service=WFS&version=2.0.0&request=GetFeature&typeName=test&outputFormat=geojson&srsName=EPSG:4326"

Operating system

ubuntu 20.04

MapServer version and installation method

mapserver 8.0.1 (from source)
proj 9.3 (from source)
gdal 3.6.2 (from source, also with proj 9.3)

[roelarents]

related: #6467