Add WFS and OAPIF support

[kannes]

@pathmapper and I were trying spatial queries on WFS layers and it might "just work". Not sure if every WFS will allow for it but it might be worth trying.

Example:
On https://www.wfs.nrw.de/geobasis/wfs_nw_dvg -> dvg:nw_dvg2_krs the filter SELECT * FROM nw_dvg1_krs WHERE ST_Intersects(msGeometry, ST_GeometryFromText('POINT(455600 5720000)')) works.

msGeometry is the name of the geometry column. We looked it up manually on https://www.wfs.nrw.de/geobasis/wfs_nw_dvg?SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=2.0.0&TYPENAMES=dvg:nw_dvg2_krs&TYPENAME=dvg:nw_dvg2_krs

nw_dvg1_krs is the name of the layer server-side.

It might be ugly/hard to get this via PyQGIS as layer.dataProvider().uri().geometryColumn() returns an empty string :(

Being able to filter (slow and quickly overloaded) WFS servers would be great!

[kannes]

The storage type is 'OGC WFS (Web Feature Service)'

OGC API Features would be useful too.

[kannes]

We need to check if the filter actually results in a filtering WFS query (server provides filtered data) or if it is just client-side filtering (full data (in bbox) is downloaded).

Does it work?
Does it work for WFS1, WFS2 and OAPIF?

[kannes]

I will look into this soonish.

[kannes]

It might be ugly/hard to get this via PyQGIS as layer.dataProvider().uri().geometryColumn() returns an empty string :(

Since QGIS 3.42/3.40.2 there is https://qgis.org/pyqgis/master/core/QgsVectorDataProvider.html#qgis.core.QgsVectorDataProvider.geometryColumnName

qgis/QGIS#59909

[pathmapper]

Since QGIS 3.42/3.40.2 there is

Nice find!

[pathmapper]

If this would result in WFS being supported by “Spatial Filter” only if the plugin is used on QGIS >= 3.40.2, here is an example of how to check the QGIS version: https://github.com/kreis-viersen/xplan-reader/blob/main/xplan_reader.py#L80-L94

[kannes]

Sadly geometryColumnName is not (yet?) implemented for the WFS and OAPIF providers, I filed a feature request at qgis/QGIS#60055

[kannes]

More findings and notes

Server-side layer/type name is available

✔️ The server-side layer/type name is available for both WFS and OAPIF via layer.dataProvider().uri().param('typename').

WFS: SQL is more elaborate and supports no ST_Transform but works

❌ ST_Intersects(msGeometry, ST_GeometryFromText('POINT(455600 5720000)')) does not work (all rows are returned)

✔️ But a filter string like SELECT * FROM "dvg:nw_dvg1_krs" WHERE ST_Intersects(msGeometry, ST_GeometryFromText('POINT(455600 5720000)')) gets turned into a FES filter and passed in the query to the WFS for server-side filtering by QGIS.
The typename must match (e.g. for this one without the dvg: prefix if using WFS 1 instead of WFS 2).

Note: QGIS will show a different query builder dialog once a filter string like this is set 🫠
Note: But it is still a normal layer.subsetString() 😃

ST_Transform does not seem to be available. It is also not listed in the functions available in the Query Builder.

WFS: Geometry column name unnecessary if we can use a QGIS expression instead of SQL

The geometry column name might not be necessary! At least if using a QGIS expression instead of a SQL-ish query.

https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#wfs-web-feature-service-data-provider-wfs says:

filter=string: QGIS expression or OGC/FES filter

and

The ‘filter’ key value can either be a QGIS expression or an OGC XML filter. If the value is set to a QGIS expression the driver will turn it into OGC XML filter before passing it to the WFS server. Beware the QGIS expression filter only supports” =, !=, <, >, <=, >=, AND, OR, NOT, LIKE, IS NULL” attribute operators, “BBOX, Disjoint, Intersects, Touches, Crosses, Contains, Overlaps, Within” spatial binary operators and the QGIS local “geomFromWKT, geomFromGML” geometry constructor functions.

https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#ogc-api-features-data-provider-oapif says

Use the special $geometry parameter to provide the layer geometry column as input into the spatial binary operators e.g intersects($geometry, geomFromWKT(‘POINT (5 6)’))

and this seems to actually be part of the WFS documentation too, accidentally only in the OAPIF section!

✔️ intersects($geometry, geomFromWKT('POINT (455600 5720000)')) does work for WFS (1 and 2), QGIS turns it into a FES filter for server-side filtering.

❌ SELECT * FROM "dvg:nw_dvg1_krs" WHERE ST_Intersects($geometry, ST_GeometryFromText('POINT(455600 5720000)')) does not work on WFS (1 and 2).

OAPIF: QGIS expression works, but only used as local (not server-side) filter

Tested with https://ogc-api.geobasis-bb.de/datasets/vg

https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#ogc-api-features-data-provider-oapif says

filter=string: QGIS expression (only datetime filtering is forwarded to the server)

and:

You can use various functions available in the QGIS Expression list, however the function must exist server side and have the same name and arguments to work.
Use the special $geometry parameter to provide the layer geometry column as input into the spatial binary operators e.g intersects($geometry, geomFromWKT(‘POINT (5 6)’))

🤷 intersects($geometry, geomFromWKT('POINT (394519 5816141)')) works but is apparently only applied locally. The amount of data downloaded matches an unfiltered layer and there is no sign of the filter in the QGIS network log entry.

TODO find out what "You can use various functions available in the QGIS Expression list, however the function must exist server side and have the same name and arguments to work." is even supposed to mean and if it actually applies to OAPIF or just WFS.
-> Asked at qgis-developer

OAPIF: SQL query seems impossible

The geometry column for bv on https://ogc-api.geobasis-bb.de/datasets/vg should be named geom according to https://www.wfs.nrw.de/geobasis/wfs_nw_dvg?SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=2.0.0&TYPENAMES=dvg:nw_dvg2_krs&TYPENAME=dvg:nw_dvg2_krs

• ❌ SELECT * FROM bv WHERE ST_Intersects(geom, ST_GeometryFromText('POINT(394519 5816141)')): Query Builder says "Error in query. The subset string could not be set."
• ❌ SELECT * FROM bv WHERE ST_Intersects($geometry, ST_GeometryFromText('POINT(394519 5816141)')): Query Builder says "Error in query. The subset string could not be set."
• ❌ ST_Intersects(geom, ST_GeometryFromText('POINT(394519 5816141)')): Query Builder says "Error in query. The subset string could not be set."
• ❌ ST_Intersects($geometry, ST_GeometryFromText('POINT(394519 5816141)')): Query Builder says "Error in query. The subset string could not be set."

I found no way to use an SQL query for OAPIF. The documentation at https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#ogc-api-features-data-provider-oapif kinda hints at a "no" since for WFS there is an extra sql parameter while for OAPIF there is only filter.

WFS and OAPIF expression "dialect" support

	QGIS Expression	"Full" SQL query
WFS 1.0.0	✔️ (turned into OGC Filter)	✔️ (turned into OGC Filter)
WFS 1.1.0	✔️ (turned into OGC Filter)	✔️ (turned into OGC Filter)
WFS 2.0.0	✔️ (turned into FES Filter)	✔️ (turned into FES Filter)
OAPIF	(✔️)(❌) (only client-side filtering)	❌

• A geometry transform function does not seem available in any of these options
• In QGIS Expressions $geometry can be used instead of the actual geometry column name

Example filters:

• QGIS Expression: intersects($geometry, geom_from_wkt('MultiPolygon (((9.2 50, 9.3 50.2, 9.5 49.9, 9.2 50)))'))
• SQL for WFS 1.*: SELECT * FROM "Bundeslaender" WHERE ST_Intersects(msGeometry,ST_GeometryFromText('MultiPolygon (((9.2 50, 9.3 50.2, 9.5 49.9, 9.2 50)))'))
• SQL for WFS 2.*: SELECT * FROM "ms:Bundeslaender" WHERE ST_Intersects(msGeometry, ST_GeometryFromText('MultiPolygon (((9.2 50, 9.3 50.2, 9.5 49.9, 9.2 50)))'))

Predicates

(ST_)Intersects, (ST_)Disjoint and (ST_)Within all seem supported with both QGIS and SQL expressions.