Inconsistencies in features count retrieved while consuming a WFS webservice

Question

GDAL approach

I'm retrieving data using GDAL with the following

ogr2ogr -f GeoJSON basol.geojson \
        --config OGR_WFS_PAGING_ALLOWED ON \
        --config OGR_WFS_PAGE_SIZE 800 \
        --config GDAL_HTTP_UNSAFESSL yes \
        WFS:"https://georisques.gouv.fr/services" \
        "ms:SSP_INSTR_GE" \
        -nln "ssp_instr_ge"

In the log, I got

ERROR 1: HTTP error code : 400
ERROR 1: Error returned by server : HTTP error code : 400 (0)

But it writes features, then on this same retrieved GeoJSON, I do the following

ogrinfo -so basol.geojson ssp_instr_ge

I get # Feature Count: 8282

QGIS approach

Then, I open QGIS, add https://georisques.gouv.fr/services as new WFS endpoint and add the same layer e.g ms:SSP_INSTR_GE. In this case, the features count is 7229 when opening the attribute table. In the metadata from the same layer in QGIS, it returns 1000

WFS resultType=hits approach

When using WFS resultType=hits https://georisques.gouv.fr/services?request=GetFeature&service=WFS&typeName=ms:SSP_INSTR_GE&version=1.1.0&resultType=hits), the value returned to count is 1000. I suppose it's "normal" as it could be relative to paging. I also suppose that's why QGIS returns 1000 in the metadata panel.

My main issue is why there is a difference between QGIS (7229 features) and GDAL (8282 features)?

Additional clues

I've also looked at geometry type to see if QGIS was doing something behind the scene but the results do not match any other numbers too (see below code)

ogrinfo basol.geojson -dialect SQLite -sql "SELECT GeometryType(GEOMETRY) AS geomtype, count(*) FROM ssp_instr_ge GROUP BY GeometryType(GEOMETRY)"

With return messages

INFO: Open of `basol.geojson'
      using driver `GeoJSON' successful.

Layer name: SELECT
Geometry: None
Feature Count: 2
Layer SRS WKT:
(unknown)
geomtype: String (0.0)
count(*): Integer (0.0)
OGRFeature(SELECT):0
  geomtype (String) = MULTIPOLYGON
  count(*) (Integer) = 823

OGRFeature(SELECT):1
  geomtype (String) = POLYGON
  count(*) (Integer) = 7459

I've also looked at geometry validity with

ogrinfo basol.geojson -dialect SQLite -sql "SELECT IsValid(GEOMETRY) as validity, count(*) FROM ssp_instr_ge GROUP BY IsValid(GEOMETRY)"

and it returns the following

Feature Count: 2
Layer SRS WKT:
(unknown)
validity: Integer (0.0)
count(*): Integer (0.0)
OGRFeature(SELECT):0
  validity (Integer) = 0
  count(*) (Integer) = 356

OGRFeature(SELECT):1
  validity (Integer) = 1
  count(*) (Integer) = 7926

I've also looked if repeated records with

ogrinfo -so basol.geojson -dialect SQLite -sql "SELECT DISTINCT code_metier, nom_usuel, adresse, code_postal, code_insee, nom_commune, nature_localisation, precision_centroide, precision_adresse, statut_id, statut_instruction, type_instr, id_inventaire_instruction, date_saisie_commune, polygone_wgs84, polygone_2154, polygone_wgs84_json, superficie, complement_adresse, polygone_2975, polygone_2972, GEOMETRY FROM ssp_instr_ge"

that returns Feature Count: 7441 amongst other infos

PS: I've also thrown below way to consume the WFS with GDAL. I mention it as it could be a clue?!

ogr2ogr -f GPKG basol.gpkg \
         --config OGR_WFS_PAGING_ALLOWED ON \
         --config OGR_WFS_PAGE_SIZE 800 \
         --config GDAL_HTTP_UNSAFESSL yes \
         WFS:"https://georisques.gouv.fr/services" \
         "ms:SSP_INSTR_GE" \
         -nln "ssp_instr_ge"

as it was returning the following error message

ERROR 1: HTTP error code : 400
ERROR 1: Error returned by server : HTTP error code : 400 (0)
ERROR 1: failed to execute insert : NOT NULL constraint failed: ssp_instr_ge.gml_id
ERROR 1: Unable to write feature 0 from layer ms:SSP_INSTR_GE.
ERROR 1: Terminating translation prematurely after failed
translation of layer ms:SSP_INSTR_GE (use -skipfailures to skip errors)
ERROR 1: sqlite3_exec(CREATE TRIGGER "trigger_insert_feature_count_ssp_instr_ge" AFTER INSERT ON "ssp_instr_ge" BEGIN UPDATE gpkg_ogr_contents SET feature_count = feature_count + 1 WHERE lower(table_name) = lower('ssp_instr_ge'); END;) failed: no such table: main.ssp_instr_ge
ERROR 1: sqlite3_exec(CREATE TRIGGER "trigger_delete_feature_count_ssp_instr_ge" AFTER DELETE ON "ssp_instr_ge" BEGIN UPDATE gpkg_ogr_contents SET feature_count = feature_count - 1 WHERE lower(table_name) = lower('ssp_instr_ge'); END;) failed: no such table: main.ssp_instr_ge

Answer 1

It may be hard to deal properly with that WFS server because it sends invalid data. You can make a test with a browser https://georisques.gouv.fr/services?SERVICE=WFS&VERSION=2.0.0&REQUEST=GetFeature&TYPENAMES=ms:SSP_INSTR_GE&STARTINDEX=0&COUNT=8

The response even contains this message that is meant to be a hint for those who administrate the service:

<!-- WARNING: No featureid defined for typename 'SSP_INSTR_GE'. Output will not validate. -->

In GML the gml:id attribute is mandatory and it must be unique for the whole dataset. Because the server in georisques.gouv.fr is not configured right it gives the same gml:id .1 for all features.

<gml:Polygon gml:id=".1" srsName="urn:ogc:def:crs:EPSG::4326">

By the error message there may even be features without any gml:id. Because of the gml:id issues the conversion into GeoPackage will fail. First because gml:id is a "not null" attribute the conversion will stop whengml:id is missing, and if you could go beyond that then the same gml:ids will break the UNIQUE constraint.

GeoJSON is a flexible format and GDAL can successfully insert the GML data into that format. The data seems to be OK - despite the non-unique gml:id - and QGIS can open the GeoJSON file that GDAL has created and both of them find the same number of features that is 8282. I would trust that GDAL does it right.

By the way, you will get a lot of information about what GDAL is doing in the background by adding --debug on into the command

ogr2ogr -f geojson basol.geojson --config OGR_WFS_PAGING_ALLOWED ON --config OGR_WFS_PAGE_SIZE 800 --config GDAL_HTTP_UNSAFESSL yes  WFS:"https://georisques.gouv.fr/services"   "ms:SSP_INSTR_GE"   -nln "ssp_instr_ge" --debug on

That reveals also that the http 400 error comes from a DescribeCoverage request. That does not seem to be critical for GDAL.

The WFS driver of QGIS seems to be flexible as well because it accepts invalid GML but I can't say why it finds less features than GDAL. I would open the GeoJSON from GDAL into the same QGIS project and compare. Perhaps missing features are outside the bbox of the QGIS map view.

I would also connect georisques.gouv.fr and ask them to fix the gml:id issue.

EDIT

The reason for different number of features in QGIS vs. GDAL is that there seems to be duplicates in the data. Probably they are not duplicates in the source data but because gml:id is not correctly configured they may appear as total duplicates. For example this query, when run against the GeoJSON data

select code_metier, adresse
 from ssp_instr_ge
where code_metier=
'SSP001083301'
order by adresse;

prints these features:

SSP001083301    1 Av Micheline
SSP001083301    1 Av Micheline
SSP001083301    183 r beranger
SSP001083301    183 r beranger
SSP001083301    2 Av Micheline
SSP001083301    2 Av Micheline
SSP001083301    3 Av Micheline
SSP001083301    3 Av Micheline
SSP001083301    3 B Av Micheline
SSP001083301    3 B Av Micheline
SSP001083301    5 Av Micheline
SSP001083301    5 Av Micheline
SSP001083301    7 Av Micheline
SSP001083301    7 Av Micheline

In the QGIS attribute table I can see only one appearance of each code_metier+adresse pair. It seems that QGIS is filtering the duplicates. That is normal behavior for the WFS driver because when it is used in the mode that is continuously requesting new features by the bounding box it is normal that already existing features are coming again from the service.

Only solid way to resolve the problem is to fix the WFS server to provide unique fids.