GetFeatureCount one off with GPKG with wkt_filter?

[edzer]

This was raised here: r-spatial/sf#2248

Expected behavior and actual behavior.

I expect GetFeatureCount to return the number of features that can be read from a layer; here is an example where it seems to be one more than that, when reading a GPKG with a Spatial Filter set. When replacing the GPKG with the fgb file with identical contents it does not happen.

Steps to reproduce the problem.

The example files are found here: https://uni-muenster.sciebo.de/s/Mn0iuk2xSkIP7D6 ; both were created by GDAL (using R package sf) see the original issue for how this was done.

#include <stdio.h>
#include <ogrsf_frmts.h>

int main() {
    OGRRegisterAll();
    GDALDataset *poDS;
    // poDS = (GDALDataset *) GDALOpenEx("cmn.fgb", GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
    poDS = (GDALDataset *) GDALOpenEx("cmn.gpkg", GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
    if (poDS == NULL) {
        printf("NULL\n");
        exit(1);
    }
    OGRLayer *poLayer = poDS->GetLayer(0);
    const char *wkt = "POLYGON ((528381.4 5073212, 537634.4 5073212, 537634.4 5083671, 528381.4 5083671, 528381.4 5073212))";
    OGRGeometry *new_geom;
    OGRErr err = OGRGeometryFactory::createFromWkt((const char *) wkt, poLayer->GetSpatialRef(), &new_geom);
    if (err != OGRERR_NONE) {
            printf("error\n");
            exit(1);
    }
    poLayer->SetSpatialFilter(new_geom);
    OGRFeature *poFeature;
    int i = 0;
    double n_d = (double) poLayer->GetFeatureCount();
    printf("feature count: %g\n", n_d);
    while ((poFeature = poLayer->GetNextFeature()) != NULL) {
        i++;
    }
    printf("number of features found: %d\n", i);
    return 0;
}

compile & run:

$ g++ -o read -I /usr/include/gdal read.cpp -lgdal ; ./read
feature count: 20
number of features found: 19

Operating system

ubuntu 22.04

GDAL version and provenance

3.6.4, but also observed with 3.7.2

[rouault]

The GetFeatureCount() implementation in the GPKG driver only did bounding box intersection testing up to now. Extended to full geometry intersection when needed in #8626

[edzer]

The spatial filter in this case is a bounding box polygon; I tried using SetSpatialFilterRect with the four corners but that gave the same problem.

[rouault]

The spatial filter in this case is a bounding box polygon; I tried using SetSpatialFilterRect with the four corners but that gave the same problem.

yes, I know. But if your feature geometry aren't rectangle, the current/past behaviour of doing bounding box intersection only doesn't lead to the same result as doing full geometry intersection

Before the fix:
$ ogrinfo ~/Téléchargements/cmn.gpkg -spat 528381.4 5073212 537634.4 5083671 -al -so|grep Count
Feature Count: 20

After:
$ ogrinfo ~/Téléchargements/cmn.gpkg -spat 528381.4 5073212 537634.4 5083671 -al -so|grep Count
Feature Count: 19

[edzer]

Ah, I hadn't seen you fixed it already - fantastic, thanks!