Features endpoint swap x and y coordinates when using EPSG CRS

[mikemahoney218-usgs]

Description

When providing no CRS or using OGC CRS strings, pygeoapi returns coordinates in (the correct) x, y order:

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&limit=1&crs=http://www.opengis.net/def/crs/OGC/1.3/CRS84

{
    "type":"FeatureCollection",
    "features":[
        {
            "type":"Feature",
            "geometry":{
                "type":"Point",
                "coordinates":[
                    5.092297735992092,
                    52.17198007698048
                ]
            },
            "properties":{ <snip> },
            "id":"Molens.1"
        }
    ],
    "links":[ <snip> ],
    "timeStamp":"2025-05-09T16:56:46.236376Z"
}

Using an EPSG CRS (which the docs list) flips these to y, x order:

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&limit=1&crs=http://www.opengis.net/def/crs/EPSG/0/4326 

{
    "type":"FeatureCollection",
    "features":[
        {
            "type":"Feature",
            "geometry":{
                "type":"Point",
                "coordinates":[
                    52.17198007698048,
                    5.092297735992092
                ]
            },
            "properties":{ <snip> },
            "id":"Molens.1"
        }
    ],
    "links":[ <snip> ],
    "timeStamp":"2025-05-09T16:57:17.131157Z"
}

Steps to Reproduce
Steps to reproduce the behavior:

Query a Features collection with crs set to an EPSG CRS.

Expected behavior
A clear and concise description of what you expected to happen.

Coordinates in GeoJSON objects should be in x, y order.

Environment

• pygeoapi version: 0.20.dev0

[mikemahoney218-usgs]

Possibly related, bbox-crs doesn't return the correct results when using an EPSG CRS:

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http://www.opengis.net/def/crs/OGC/1.3/CRS84

{
    "type":"FeatureCollection",
    "features":[
        {
            "type":"Feature",
            "geometry":{
                "type":"Point",
                "coordinates":[
                    5.092297735992092,
                    52.17198007698048
                ]
            },
            "properties":{ <snip > },
            "id":"Molens.1"
        },
        {
            "type":"Feature",
            "geometry":{
                "type":"Point",
                "coordinates":[
                    5.057482747807393,
                    52.23651108372093
                ]
            },
<snip>
}

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http://www.opengis.net/def/crs/EPSG/0/4326

{
    "type":"FeatureCollection",
    "features":[],
    "links":[
        {
            "type":"application/geo+json",
            "rel":"self",
            "title":"This document as GeoJSON",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326"
        },
        {
            "rel":"alternate",
            "type":"application/ld+json",
            "title":"This document as RDF (JSON-LD)",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=jsonld&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326"
        },
        {
            "type":"text/html",
            "rel":"alternate",
            "title":"This document as HTML",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=html&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326"
        },
        {
            "type":"application/json",
            "title":"Windmills within The Netherlands",
            "rel":"collection",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills"
        }
    ],
    "timeStamp":"2025-05-09T17:08:37.672879Z"
}

That said, they do return identical results when resulttype=hits:

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http://www.opengis.net/def/crs/EPSG/0/4326&resulttype=hits

{
    "type":"FeatureCollection",
    "numberMatched":1872,
    "features":[],
    "links":[
        {
            "type":"application/geo+json",
            "rel":"self",
            "title":"This document as GeoJSON",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326&resulttype=hits"
        },
        {
            "rel":"alternate",
            "type":"application/ld+json",
            "title":"This document as RDF (JSON-LD)",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=jsonld&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326&resulttype=hits"
        },
        {
            "type":"text/html",
            "rel":"alternate",
            "title":"This document as HTML",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=html&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326&resulttype=hits"
        },
        {
            "type":"application/geo+json",
            "rel":"next",
            "title":"Items (next)",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills/items?offset=10&bbox=3.37,50.75,7.21,53.47&bbox-crs=http%3A%2F%2Fwww.opengis.net%2Fdef%2Fcrs%2FEPSG%2F0%2F4326&resulttype=hits"
        },
        {
            "type":"application/json",
            "title":"Windmills within The Netherlands",
            "rel":"collection",
            "href":"https://demo.pygeoapi.io/master/collections/dutch_windmills"
        }
    ],
    "timeStamp":"2025-05-09T17:09:10.170692Z"
}

[mikemahoney218-usgs]

Okay, looking into this more, I think I'm having trouble understanding the intention of Features Part 2 here.

The proximate issue is that these EPSG CRS are lat/lon ordered rather than lon/lat:

>>> import pyproj
>>> url ="https://www.opengis.net/def/crs/OGC/1.3/CRS84"
>>> authority, code = url.rsplit("/", maxsplit=3)[1::2]
>>> crs = pyproj.CRS.from_authority(authority, code)
>>> crs
<Geographic 2D CRS: OGC:CRS84>
Name: WGS 84 (CRS84)
Axis Info [ellipsoidal]:
- Lon[east]: Geodetic longitude (degree)
- Lat[north]: Geodetic latitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> epsg_url='https://www.opengis.net/def/crs/EPSG/0/4326'
>>> epsg_authority, epsg_code = epsg_url.rsplit("/", maxsplit=3)[1::2]
>>> pyproj.CRS.from_authority(epsg_authority, epsg_code)
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Because always_xy defaults to False, I suppose coordinates are getting flipped:
https://github.com/geopython/pygeoapi/blob/master/pygeoapi/util.py#L765

This then makes the output invalid GeoJSON (and also messes up the map display in the HTML).

Features Part 2 obliquely suggests this might be an issue for HTML responses:

HTML does not have any provisions for spatial geometries and coordinate reference systems. However, note that schema.org that is embedded in HTML only supports WGS 84 in the axis order latitude/longitude, so any coordinates in schema.org markup will have to be in that coordinate reference system, independent of the requested coordinate reference system.

But there's no discussion of how coordinate ordering should be handled in GeoJSON. This despite GeoJSON explicitly requiring x/y ordering in Section 3.1.1:

The first two elements are longitude and latitude, or easting and northing, precisely in that order and using decimal numbers.

I personally don't read the provisions of Section 4 of the GeoJSON spec as allowing alternate coordinate orderings. It also seems like users would have some challenge with trying to process GeoJSON with flipped axis ordering using standard tools.

All that to say, this might not be a bug with how Part 2 has been implemented, but does seem like a less usable implementation than preserving x/y ordering. I'm wondering if there's any examples -- or any interest -- in following either Case 2 or Case 4 of "Revision to Axis Order Policy and Recommendations" https://portal.opengeospatial.org/files/?artifact_id=76024 (word doc):

Case 2: Payload encoding explicitly overrides the Axis Order as specified in the CRS metadata: In this case the payload explicitly encodes additional metadata in order to determine axis order. For example, in the draft ISO 19107 (Draft International Standard), axis "reordering" may be included as a CS transformation in payload comments. The included CS transformation must reference the axis names actually defined by the CRS metadata.

Case 4: The Encoding specification (document) explicitly overrides the Axis Order as specified in the CRS metadata: In this case, a specification (document) specifies a default CRS but then overrides the axis order of the CRS as specified in the CRS metadata. A typical example is to specify WGS 84 (EPSG 4326) in the specification but then state in the document that the axis order is the opposite of the CRS metadata. The payloads would then be encoded using the axis order as specified in the spec document and not in the CRS metadata. The OGC does not recommend implementing this use case.

Either of these approaches -- adding metadata announcing x/y ordering or relying on the GeoJSON standard to override the axis order -- seems like it would help address this issue. I understand if that's too far a deviation from Part 2 as written, though -- in which case I think the server I'm responsible for probably won't support non x/y ordered CRS.

[mikemahoney218-usgs]

One last comment, to expand upon what I mean by "challenge with trying to process GeoJSON with flipped axis ordering using standard tools".

Running ogr2ogr windmills.gpkg https://demo.pygeoapi.io/master/collections/dutch_windmills/items\?f\=json\&crs\=http://www.opengis.net/def/crs/EPSG/0/4326 places the /dutch_windmills collection in the orange spot below, vs the purple spot if you don't provide a CRS argument:

Similarly, programming languages which use GDAL for IO wind up corrupting the coordinates:

library(sf)
windmills <- read_sf(
  "https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&crs=http://www.opengis.net/def/crs/EPSG/0/4326"
)
bounds <- read_sf("~/Downloads/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")
plot(sf::st_geometry(bounds))
plot(sf::st_geometry(windmills), col = "red", fill = "red", add = TRUE)

[justb4]

@mikemahoney218-usgs @tomkralidis Axis ordering (x,y, or lon,lat) as in the examples above can be dependent on the format of the projection/CRS string even in GeoJSON as returned by pygeoapi. We went through great length to follow the standards when implementing CRS-support in pygeoapi. See for example:

• CRS support PR comments
• these unit tests

So for example:

http://www.opengis.net/def/crs/OGC/1.3/CRS84 --> lon, lat
http://www.opengis.net/def/crs/EPSG/0/4326':   --> lat, lon
http://www.opengis.net/def/crs/EPSG/0/28992   --> x,y

[mikemahoney218-usgs]

Yeah, apologies, I only realized that the CRS were specifying the axis ordering in the third extremely long comment on this issue #2010 (comment)

My question now is if, when a CRS specifies y/x ordering, it's better to return invalid GeoJSON (as pygeoapi currently does), or look to the "Revision to Axis Order Policy and Recommendations" document (linked in the bibliography of part 2) for a way that would reorder coordinates to return valid GeoJSON. What gets returned now violates section 3.1.1:

The first two elements are longitude and latitude, or easting and northing, precisely in that order and using decimal numbers.

The interoperability issues listed above stem from this.

Separately, I think there might still be an issue with bbox-crs as discussed in the second comment, where different numbers of results are matched when resulttype=hits vs not.

[justb4]

@mikemahoney218-usgs no need for apologies! We're all in the same (rocky) boat. In The Netherlands it is even common geo-knowledge that when your data ends up near the coast of Somalia that there is an Axis ordering issue :-(. Often a decision is a Catch-22 like the GeoJSON case. I don't have a good answer on that.

But the bbox-crs with/without hits in your second comment may be a bug (as I think 'hits' should be 0, right, CRS with lat, lon thus Somalia?):

curl https://demo.pygeoapi.io/master/collections/dutch_windmills/items?f=json&bbox=3.37,50.75,7.21,53.47&bbox-crs=http://www.opengis.net/def/crs/EPSG/0/4326&resulttype=hits

Can you open a new issue for the `bbox-crs hits problem' ?

[mikemahoney218-usgs]

Opened in #2027