Allow antimeridian bbox on OGC Features API

[geographika]

As reported by @jratike80 in the mailing list - https://lists.osgeo.org/pipermail/mapserver-dev/2023-October/017012.html
MapServer currently fails one of the OGC Teamengine tests for OGC Features API:

For WGS 84 longitude/latitude the bounding box is in most cases the sequence of minimum longitude, minimum latitude, maximum longitude and maximum latitude. However, in cases where the box spans the anti-meridian the first value (west-most box edge) is larger than the third value (east-most box edge).

Example 6. The bounding box of the New Zealand Exclusive Economic Zone
The bounding box of the New Zealand Exclusive Economic Zone in WGS 84 (from 160.6°E to 170°W and from 55.95°S to 25.89°S) would be represented in JSON as [ 160.6, -55.95, -170, -25.89 ] and in a query as bbox=160.6,-55.95,-170,-25.89.

A couple of items to resolve:

1. Should this switch only be valid for bounding-boxes in EPSG:4326, or would other coordinate systems have similar issues?
2. The current test data uses Minnesota counties so no results are returned (but there is no longer an error). It might be worth adding a JSON world cities dataset for testing this in a new Mapfile?

[rouault]

Unfortunately the fix would be much more complicated than just switching minx and maxx.

We should actually issue 2 queries:

• one with [bbox->minx, 180] range
• another one with [-180, bbox->maxx] range

and merge the results

[geographika]

Ah of course, this is going to get the completely wrong results.. I presume this issue wasn't addressed in WFS queries? Will look now.

[rouault]

I presume this issue wasn't addressed in WFS queries?

no, WFS hasn't any such special processing for antimeridian crossing

[rouault]

• Should this switch only be valid for bounding-boxes in EPSG:4326, or would other coordinate systems have similar issues?

mostly/only for geographic CRS

Some projected CRS like EPSG:3857 could also have wraparound (at eastings +/- 20037508.3427892 for EPSG:3857), but that seems rather esoteric in a OGCAPI Features context

[geographika]

I've found discussions on this test at opengeospatial/ets-ogcapi-features10#187 and the pull request that resolves it in pygeoapi - geopython/pygeoapi#837
Although the pull request stops the service from throwing an error, I'm not sure the rest of the providers in pygeoapi will handle getting the correct results.

[geographika]

From @msmitherdc - splitting the geometries using GDAL (see script below) the cross the antimeridian might be an option. Python code using GDAL API below for reference:

def ogr_meridian_split(geojson: 'Union[str, gdal.DataSet]') -> 'ogr.Geometry':
    """
    Splits a polygon which crosses the antimeridian.
 
    Args:
        geojson (Union[str, gdal.DataSet]): Geojson string, or gdal dataset.
 
    Returns:
        ogr.Geometry:OGR geometry
    """
 
    if isinstance(geojson, str):
        in_ds = gdal.OpenEx(geojson)
    else:
        in_ds = geojson
    out_ds = gdal.VectorTranslate("/vsimem/out.json",in_ds, format="GeoJSON", layerCreationOptions=["RFC7946=YES", "WRAPDATELINE=YES", "WRITE_BBOX=YES"])
    del out_ds
    newds = gdal.OpenEx("/vsimem/out.json")
    lyr = newds.GetLayer(0)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    wkt = geom.ExportToWkt()
    split_ogr_poly = ogr.CreateGeometryFromWkt(wkt)
    del newds
    return split_ogr_poly

[geographika]

Closing this and moving to an issue - #6981