Flipped axis order from postgis feature provider using EPSG:4269 as provider-crs.

[dblodgett-usgs]

Description
We have a QA-tier server up hosting a number of large EPSG:4269 (NAD83) tables.

https://labs-beta.waterdata.usgs.gov/api/fabric/pygeoapi/

config.yml for the server is here.

https://code.usgs.gov/wma/nhgf/pygeoapi/-/blob/main/configs/pygeoapi/fabric.yml?ref_type=heads

We've observed that the layers can be accessed via bbox, but axis order is flipped -- both in the bbox url field and in the database.

Steps to Reproduce

Go to: https://labs-beta.waterdata.usgs.gov/api/fabric/pygeoapi/collections/nhdplusv2-huc08/items?bbox=42,-94,43,-93

Observe that the order of coordinates makes no sense but the server is behaving normally.

Observe that no features are returned for: https://labs-beta.waterdata.usgs.gov/api/fabric/pygeoapi/collections/nhdplusv2-huc08/items?bbox=-94,42,-93,43

Expected behavior

When a postgres table is in 4269, it can be interacted with as normal lon/lat order CRS84 geojson.

Additional context

Related to #2010 ??

[dblodgett-usgs]

I've been working this issue with @webb-ben and @ptomasula --

We found that PostGIS WKT/WKB for EPSG:4269 is returned from a query in X/Y (lon/lat) order but pygeoapi's use of pyproj is assuming it is in the order that the coordinates appear in the CS block of the EPSG:4269 projection WKT.

From what I can tell, it comes down to the fact that postgis encodes WKB according to prevailing interpretation of the WKB specification (SFSQL) where axis order is always x,y and "longitude" is nominally/semantically an X coordinate. The WKB axis order note in the geopackage specification is instructive (see section 2.1.3 of the geopackage spec).

As a result, whenever pygeoapi retrieves geometry from postgis, it needs to honor the axis order the database delivers it in.

--- a digression ---

It's strange that CRS transformations and geojson serialization aren't delegated to PostGIS. By mixing PostGIS geometry and python-package operations, there is a near gauruntee that inconsistencies and bugs will be introduced.

--- end digression ---

To "fix" the bug we are experiencing, flipping always_xy to TRUE whenever pyproj.Transformer.from_crs [1] [2] is called in the context of a provider that works in WKT/WKB is in order. Both instances are in util.py so the calling code is going to have to specify the value of that method parameter.

[webb-ben]

This is not a perfect solution because it only addresses the manifestation and partially ignores https://proj.org/en/stable/faq.html#why-is-the-axis-ordering-in-proj-not-consistent). This could be controlled at runtime with an optional PYGEOAPI_ALWAYS_XY environment variable that would propagate to all invocations of PROJ transform, thereby allowing this type of granular level of control.

pygeoapi generally has avoided env variables for control but I am not sure there is a solution that will happily satisfy legacy systems that demand lon/lat regardless of the content crs while making use of axis ordering from the authority.

[mikemahoney218-usgs]

Could this be set in the server block of config.yml instead of requiring an env var?

[dblodgett-usgs]

We need to be very granular here.

PROJ respects the axis ordering as it was defined by the authority in charge of a given coordinate reference system.

And, "authority in charge of a given coordinate reference system" is NOT the encoding of the coordinate data in question.

By separating the encoding from the crs authority, you have to allow for encodings that don't change axis order based on the CRS. So, at some level, there needs to be logic that determines if a given set of coordinates are adherent to the CRS axis order or are always_xy = TRUE. The predominant case is with WKT/WKB -- it is always_xy = TRUE until someone decides that "latitude" can possibly be semantically interpreted as a X coordinate (never?).

[dblodgett-usgs]

Could this be set in the server block of config.yml instead of requiring an env var?

I think that just papers over the root of the issue. It's not impossible to get this right, you just have to follow SFSQL and know when you are dealing with a payload that's tightly coupled to the crs authority or not.