CRS Support for OGC API Feature pygeoapi Provider

[justb4]

Overview

See PRs #1155 and #1149.
This PR combines these two as together they provide support
for the complete standard OGC API - Features - Part 2: Coordinate Reference Systems by Reference.

Related Issue / Discussion

Issue #1128 - See PRs 1155 and #1149.

Additional Information

See related PRs above.

Contributions and Licensing

(as per https://github.com/geopython/pygeoapi/blob/master/CONTRIBUTING.md#contributions-and-licensing)

• I'd like to contribute this feature to pygeoapi. I confirm that my contributions to pygeoapi will be compatible with the pygeoapi license guidelines at the time of contribution.
• I have already previously agreed to the pygeoapi Contributions and Licensing Guidelines

[justb4]

The CI is failing from what I think is a Proj/GDAL version issue. The test that fails is the OGR Provider which uses internal GDAL calls for the Transform (i.s.o. using Shapely/pyproj as in api.py) when a CrsTransform is passed on get() or query().

On my local system this unit test succeeds with versions:

GDAL 3.2.2, released 2021/03/05
pyproj==3.2.1
libproj 7.2.1, January 1st, 2021

The CI uses Ubuntu 20.04 without UbuntuGIS so may have different versions of these.

[MTachon]

Yes, reading about it, it definitely sounds like a GDAL/PROJ version issue. GDAL seems to throw an error when calling osr.CoordinateTransformation() using the WKT representation of CRSs exported by the given version of pyproj (when creating a CrsTransformWkt instance).

[justb4]

Probably best is to upgrade main.yml from Ubuntu 20.04 to 22.04. This is also inline with the Docker Image. As a start.
Test runs fine in Docker Image. I opened #1176 as it has quite some impact.

platform linux -- Python 3.10.6, pytest-6.2.5, py-1.10.0, pluggy-0.13.0
rootdir: /pygeoapi, configfile: pytest.ini
collected 17 items

tests/test_ogr_shapefile_provider.py .................                   [100%]

But interesting, 2 errors in the WFS tests with the GeoSolutions WFS:

platform linux -- Python 3.10.6, pytest-6.2.5, py-1.10.0, pluggy-0.13.0
rootdir: /pygeoapi, configfile: pytest.ini
collected 17 items

tests/test_ogr_shapefile_provider.py .................                   [100%]

[justb4]

There is some confusion about axis-ordering in the GeoJSON response with the crs parameter.
I think that the ordering should follow the definition of the CRS URI. Others say that the GeoJSON convention (lon, lat) needs to be followed, or at least are not sure.
I tested on the ldproxy demo and see that it follows 'my' convention:

https://demo.ldproxy.net/daraa/collections/AgriculturePnt/items?f=json&crs=https://www.opengis.net/def/crs/EPSG/0/4326 (lat, lon) , would be lat, lon for e.g. 4258 (ETRS89 INSPIRE)
https://demo.ldproxy.net/daraa/collections/AgriculturePnt/items?f=json&crs=https://www.opengis.net/def/crs/OGC/1.3/CRS84 (lon, lat)

Think pygeoapi needs to follow CRS Axis ordering in GeoJSON responses. Looks like all CRS is enforced xy order (using pyproj+Proj):
https://apitestbed.geonovum.nl/pygeoapi/collections/AddressesNL/items?crs=https://www.opengis.net/def/crs/OGC/1.3/CRS84&f=json (lon, lat) OK

https://apitestbed.geonovum.nl/pygeoapi/collections/AddressesNL/items?crs=https://www.opengis.net/def/crs/EPSG/0/4258&f=json (lon, lat) should be reversed IMO.

https://apitestbed.geonovum.nl/pygeoapi/collections/AddressesNL/items?crs=https://www.opengis.net/def/crs/EPSG/0/4326&f=json (lon, lat) should be reversed IMO.

https://apitestbed.geonovum.nl/pygeoapi/collections/AddressesNL/items?crs=https://www.opengis.net/def/crs/EPSG/0/28992&f=json x,y Dutch OK

For bbox-crs we also follow the CRS URI conventions.

[MTachon]

I agree with @justb4 , that the axis ordering in GeoJSON responses should follow the definition of the CRS, and not forced to lon/lat for geographic coordinate systems.

The section Coordinate reference system information independent of the feature encoding introduces the Content-Crs header as a solution for "...the continued confusion caused by the axis order of coordinates, ...". It defines it as "a mechanism for a server to clearly and unambiguously assert the CRS being used in a response document independent of the requested output format.". "Natively", GeoJSON does not support other CRSs than WGS84 lon/lat (OGC:CRS84), but the Part 2 of the standard OGC API - Features mentions the "the prior arrangement provision in the second paragraph of section 4 of the GeoJSON standard", which "allows other coordinate systems to be used". It also states: "An explicit request by a client with a query parameter crs establishes a prior arrangement. It is the responsibility of the client that submits the request to handle the coordinates in the response correctly. ". By handling the coordinates correctly, I would assume that axis ordering of the features' coordinates in the response should follow that of the queried CRS.

As for a solution, pyproj can perfectly handle the axis ordering for different CRSs with the always_xy parameter:

# Input CRS with lat/lon axis order
crs_in = pyproj.CRS.from_espg(4326)

# Output CRS with lon/lat axis order
crs_out = pyproj.CRS.from_authority("OGC", "CRS84")

# Create a pyproj.Transformer instance which forces lon/lat ordering
transformer_1 = pyproj.Transformer.from_crs(crs_in, crs_out, always_xy=True)

# Create a pyproj.Transformer instance which takes care of changing axis order if needed
transformer_2 = pyproj.Transformer.from_crs(crs_in, crs_out, always_xy=False)

>>> transformer_1.transform(5.0, 6.0)
(5.0, 6.0)

>>> transformer_2.transform(5.0, 6.0)
(6.0, 5.0)

On the other hand, when looking at the implementation of the OGRProvider, I see that it is using the SetAxisMappingStrategy method of the OGRSpatialReference to set the axis order to osr.OAMS_TRADITIONAL_GIS_ORDER, which forces the lon/lat axis order for coordinates transformation. In other words, the transformation will expect coordinates in lon/lat order as input (independently of the data source and its support for different coordinate axis order for different CRSs), and spit out coordinates in the same order. As I understand it, the GDAL/OGR OGRSpatialReference.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) is the same as the pyproj always_xy=True. If we decide to change one of them, the other should be changed accordingly.

If we decide to follow @justb4 's convention, setting always_xy=False in util.get_transform_from_crs and util.transform_bbox, as well as removing the occurences of OGRSpatialReference.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) in the OGRProvider may be enough?

As for data stored in PostGIS, it looks like the ST_Transform function switches the coordinates axis order if needed, so order should matter when inserting/getting geometries from there.

[MTachon]

Something else that came to my mind when reading the Part 2 of the standard OGC API - Features, is the implementation of pygeoapi when serving data in JSON for Linking Data format (JSON-LD). Currently, the linked_data.geojson2jsonld function uses the X/Y ordering, stored in shapely geometrical object. Coordinates order should be switched in this function, as coordinates in shema.org markup (GeoShape/GeoCoordinates) only supports WGS 84 in the axis order latitude/longitude. Part 2 of the standard also says "...any coordinates in schema.org markup will have to be in that coordinate reference system, independent of the requested coordinate reference system.". Depending on which crs is queried, or which storage_crs is used, if the user chooses the JSON-LD data format, the coordinates in the schema.org markup may need be transformed to EPSG:4326...

[justb4]

@MTachon wrote:

If we decide to follow @justb4 's convention, setting always_xy=False in util.get_transform_from_crs and > util.transform_bbox, as well as removing the occurences of > OGRSpatialReference.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) in the OGRProvider may be enough?

Yes, this is my current plan for this PR/branch, and get all tests working, maybe add a few more. Like with BBOX-CRS I have several to check for CRS-compliant Axis ordering.

[justb4]

Somehow the CRS Axis ordering when using OGR Transforms is not honoured it looks like. Maybe has to do as the it projection is loaded from WKT as passed to get() and query()`? It seems to work when loaded from an EPSG Code. For now disabling native transform until sorted out.

[MTachon]

I have added some changes regarding JSON-LD and schema.org markup in branch https://github.com/MTachon/pygeoapi/tree/crs-features-ogc

[MTachon]

I do not think that it has something to do with the projection being loaded from WKT, but more with osr.OAMS_AUTHORITY_COMPLIANT vs. osr.OAMS_TRADITIONAL_GIS_ORDER in the _get_spatial_reference(). See my comment for ogr.py file.

[justb4]

I do not think that it has something to do with the projection being loaded from WKT, but more with osr.OAMS_AUTHORITY_COMPLIANT vs. osr.OAMS_TRADITIONAL_GIS_ORDER in the _get_spatial_reference(). See my comment for ogr.py file.

Have you seen the latest changes in ogr.py, on line 79:
CRS_AXIS_ORDER = osgeo_osr.OAMS_AUTHORITY_COMPLIANT
This setting is used for all Projection loading, EPSG (except 4326) used with the source/target EPSG config and WKT as for the Transformation object passed in get() and query(). Still the Projection/Transformation using the latter WKT loading does not honour CRS Axis ordering, e.g. for 4258. For now I have disabled Transformation for ogr.py and use @crs_transform until sorted out. I think the WKT loading also caused the CI main.yml to fail (vague error pointing to 'WKT1 2019' or something). Unless someone can fix this, I propose to open a new issue/PR. Maybe the CrsTransformation Object can also specify/embed two EPSG code-ints, such that ogr.py does loading from EPSG i.s.o. WKT as the latter seems 'brittle' or maybe even the full Projection URIs, maybe that works as well. Other Providers may be helped with this as well.

Only 'gotcha' now may be storageCRS which should be set now to the 'target' EPSG, as in some cases there may be two reprojections.

[MTachon]

So if I understand correctly, importing a CRS (into a OGRSpatialReference object) from WKT gives a different result than when importing a CRS from EPSG? Or do you get an error when mixing source_srs/target_srs and crs/storage_crs? What dataset are you using for testing (and which CRS and axis ordering does it have)? What is the output CRS wanted? The CI main.yml may be related to a WKT version issue, resulting from different pyproj's PROJ and GDAL's PROJ versions. Specifying the WKT version can be done with pyproj when exporting CRSs (and constructing the CrsTransformWkt instance). If changing this works, it may be an easier fix than fixing main.yml.

[justb4]

So if I understand correctly, importing a CRS (into a OGRSpatialReference object) from WKT gives a different result than when importing a CRS from EPSG? Or do you get an error when mixing source_srs/target_srs and crs/storage_crs?

Yes, that looks like it: importing from WKT gives either wrong (Axis ordering) result in Docker deploy or failure in main.yml CI. Also when disabling source_srs/target_srs. Happens with any storage_crs: https://www.opengis.net/def/crs/OGC/1.3/CRS84 (stored in lon,lat) or https://www.opengis.net/def/crs/EPSG/0/28992 (always x,y).

What dataset are you using for testing (and which CRS and axis ordering does it have)? What is the output CRS wanted? The CI main.yml may be related to a WKT version issue, resulting from different pyproj's PROJ and GDAL's PROJ versions. Specifying the WKT version can be done with pyproj when exporting CRSs (and constructing the CrsTransformWkt instance). If changing this works, it may be an easier fix than fixing main.yml.

It follows from the unit tests. Using OGR with GPKG with several projections. You can see the testfile/URL here:
https://apitestbed.geonovum.nl/test/ . All CRS, BBOX-CRS tests here work as expected now, as the OGR Provider delegates transforms to core (via @crs_transform decorator). The pygeoapi instance runs from the repo branch crs-features-ogc-adr as Docker instance. Find data and config here.

My personal opinion is that WKT loading seems brittle among pyproj, Proj, GDAL. Users will run different versions. Will be hard to grasp this. That is why I suggested to to pass a Transform object with also the original EPSG code or maybe URI to get() and query(). At least within OGR Driver internal Projection instantiation is uniform, both via source/target and crs parameter, i,e, load from EPSG with osgeo_osr.OAMS_AUTHORITY_COMPLIANT.

[justb4]

With commit fd16159 All OGC CITE tests for Part 1 and Part 2 are now succeeding!

[MTachon]

@justb4 That's great, nice work!! :)

[tomkralidis]

Does this decorate mean an extra hop for all provider query and get functions?

[MTachon]

The decorator iterates through the features returned by the providers' query and get functions, transform the geometries' coordinates if needed (meaning that crs_transform_spec != None in query() and get() calls in api.py) before the Feature/FeatureCollection is assigned to the content variable in api.py. The decorator approach was used to avoid code duplication (fewer lines of code if the transformation logic is taken out of the providers' implementation), for simplicity and for supporting coordinates transformation for all providers. PostgreSQLProvider has its own implementation and does not use the @crs_transform decorator. The implementation for the PostgreSQL provider is more optimized as it does not re-iterate through the returned features, and we may want to optimize the coordinates transformation implementation of other providers later, especially if these providers have native support for coordinates transformation. In practice, re-iterating through some thousands of features would not result in much overhead.

[justb4]

Yes, agree. Also the OGRProvider does its own transforms. Our plan is to first have #1174 merged, having a stable version, CITE and INSPIRE compliant. Then do refactoring in separate PRs:

• move all CRS-dedicated/related functionality and constants spread over util.py and api.py into a new module pygeoapi/crs.py
• enhance CrsTransformSpec with direct-to-call transform functions as partials (e.g. for Shapely and GeoJSON dicts) that Providers without native CRS-support can call in their Feature-Collection loop
• upgrade Providers to using the CrsTransformSpec-based transforms, removing the Decorator, maybe even completely

[GeoSander]

@justb4 Now that #1152 has been merged, please let me know if you need a hand resolving conflicts.

[justb4]

Thanks @GeoSander ! Yes, will be quite a merge from master, but seen these changes from #1152 before, IntelliJ will be my friend. I'll let know.

[GeoSander]

FYI @justb4: I did a working tree diff with the current master this morning and did not find anything wrong or missing (during the merge), so in that sense it looks good to me 👍
Yet it cannot be merged now due to more conflicts..?

[justb4]

@GeoSander that is good to hear! Hopefully @tomkralidis finds some time to look at the last resolutions from the review.

[tomkralidis]

CI failures are unrelated (https://lists.osgeo.org/pipermail/mapserver-users/2023-April/083094.html).