rectObj wrong projection with PROJ > 6

[WouterVisscher]

in short

When executing a WMS/WFS request in a different projection as the source a inaccured rectObj is generated from the &BBOX param and used in the backend query. Resulting in 'missing' objects in a image or features in a FeatureCollection.

some background

We are preparing our stacks for a GDAL/PROJ upgrade when we noticed this behaviour.

the issue

So it seems that with our Mapserver stack + PROJ > 6 (and recent GDAL version) this behaviour manifests.
This doesn't seem to be the case when using PROJ < 6 (like PROJ 5.2.0). The type of backend also doesn't really seem to matter, OGR(GPKG/GeoJSON), PostGIS.

Interesting enough the results, image (for WMS) and features (for WFS), don't seem to be affected. In the sense that the objects in the PNG are in the right place and the features return by the WFS are also projected correctly. They are 'only' missing objects that should have been there.

So the green bbox is the BBOX in EPSG:3857 from the request (WMS or WFS)
And the grey bbox is the one generated by Mapserver and used in the backend for requesting the features.

Logging of the WMS/WFS request and generated rectObj

This effect creates this kind of object popping in a WMS. Because the deviation is always the same it's 'only' really noticeable in the lower zoomlevels. The effects of this miscasting of the BBOX also has a effect on WMS GetFeatureInfo requests.

To be honest this could stil maybe be a config/setting issue in our setup, but because the results that are given back are projected (EPSG:28992 > EPSG:3857) correctly we have a feeling something might be of with the creation of the rectObj.

Tech used

Mapserver main, commit 68de818
PROJ Rel. 8.2.1, January 1st, 2022
GDAL GDAL 3.4.1, released 2021/12/27
Debian bullseye

Also with older versions of proj (like 7.1.1) it seems to happen.
And with the latest Mapserver release https://github.com/MapServer/MapServer/tree/rel-7-6-4
Mapfile and data used in example ams.updated.zip

[geographika]

Possibly related - the pull request #6472 wouldn't apply to this issue, but it may be very similar and missing a reprojection call.

[jmckenna]

@WouterVisscher can you attach a small test case as a zip archive here? (mapfile snippet, data sample, commands.txt containing the commands to reproduce - you can use mapserv -nh QUERY_STRING="" for your problematic GetMap request and include it in commands.txt

I'm also about to report a PROJ >=6 issue, maybe all 3 reports/findings are the same.

[jmckenna]

Opps I see your zip above, thanks! Sorry I missed it.

[WouterVisscher]

For completion, the mapserver config file

CONFIG
  ENV
    MS_MAP_NO_PATH "true"
  END
  MAPS
    MAP1 "/srv/data/ams.map"
  END
END

query strings use for testing:
WFS:
QUERY_STRING="SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAMES=bag:pand&STARTINDEX=0&COUNT=1000&SRSNAME=urn:ogc:def:crs:EPSG::3857&BBOX=545287.2873572960962,6867556.049125162885,545689.3385149866808,6868025.580225903541,urn:ogc:def:crs:EPSG::3857"

WMS: QUERY_STRING="SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=545287.2873572960962,6867556.049125162885,545689.3385149866808,6868025.580225903541&CRS=EPSG:3857&WIDTH=995&HEIGHT=1162&LAYERS=pand&STYLES=&FORMAT=image/png&TRANSPARENT=TRUE"

[jmckenna]

Thanks @WouterVisscher, appreciated.
(for your next ticket please include the commands.txt in the zip)

[WouterVisscher]

@jmckenna I can make a new zip with config + commands.txt if you want (and data + mapfile of course)

[jmckenna]

@WouterVisscher might be good just to make it easier if someone else wants to test. Thanks.

[WouterVisscher]

@jmckenna updated the main post with a new zip

[jmckenna]

@WouterVisscher thanks, sorry to be a pain, ha. Thanks for this.

[rouault]

Fix in #6479 + OSGeo/PROJ#3055
A workaround without them might be to expand "+init=epsg:28992" to " +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs"

[WouterVisscher]

That's fast! 👍🏻
Yup, using "+proj=sterea +lat_0=...." seems to mitigate the behaviour.

[geographika]

@rouault - I'm currently fixing the MapServer Appveyor builds, which rely on Proj 9.1 in the GISInternals SDK. The test added to check this issue is failing. See https://ci.appveyor.com/project/SethG/mapserver/build/job/wm0wngk7mg169uoh#L610

The point seems to be reprojected to different coordinates than the lower-left point of the rectangle. The assertion is false so the test isn't skipped, but fails on assert rect.minx == pytest.approx(121711, abs=2)

def test_reprojection_rect_and_datum_shift():
    
        webmercator = mapscript.projectionObj("init=epsg:3857")
        epsg_28992 = mapscript.projectionObj("init=epsg:28992") # "Amersfoort / RD New"
    
        point = mapscript.pointObj(545287, 6867556)
        point.project(webmercator, epsg_28992)
        if point.x == pytest.approx(121685, abs=2):
            # Builds of PROJ >= 6 and < 8 with -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H
            # use pj_transform() but it doesn't work well with datum shift
            # This is a non-nominal configuration used in some CI confs when use a
            # PROJ 7 build to test PROJ.4 API and PROJ 6 API.
            pytest.skip('This test cannot run with PROJ [6,8[ in builds with ACCEPT_USE_OF_DEPRECATED_PROJ_API_H')
    
        rect = mapscript.rectObj(545287, 6867556, 545689, 6868025)
        assert rect.project(webmercator, epsg_28992) == 0
>       assert rect.minx == pytest.approx(121711, abs=2)
E       AssertionError: assert 121685.76222976402 == 121711 +- 2.0e+00
E        +  where 121685.76222976402 = <mapscript.rectObj; proxy of <Swig Object of type 'rectObj *' at 0x03B3A920> >.minx
E        +  and   121711 +- 2.0e+00 = <function approx at 0x0370F8F0>(121711, abs=2)
E        +    where <function approx at 0x0370F8F0> = pytest.approx

[rouault]

@geographika No idea comes to mind

[jmckenna]

linked to OSGeo/gdal#7179

[rouault]

linked to OSGeo/gdal#7179

I absoluetely see no connection between a PROJ releated issue and a MSSQL authentication one