Filtering OGRLayer Features Using GEOSGeometry in C

In C, filter the features in a OGRLayer by a GEOSGeometry

To filter features in an OGRLayer by a GEOSGeometry in C, you can use the GDAL/OGR library along with the GEOS library. The process involves creating a spatial filter using the GEOSGeometry and applying it to the OGRLayer. Below is a step-by-step guide and example code to achieve this:

1. Include necessary headers: Make sure to include the headers for GDAL/OGR and GEOS.

2. Initialize GDAL/OGR: Initialize the GDAL/OGR library.

3. Open the data source and get the layer: Open the data source (e.g., a shapefile) and get the layer you want to filter.

4. Create a GEOSGeometry: Create or read the GEOSGeometry that you will use as the spatial filter.

5. Convert GEOSGeometry to OGRGeometry: Convert the GEOSGeometry to an OGRGeometry because OGR uses its own geometry classes.

6. Set the spatial filter on the layer: Apply the spatial filter to the layer using the OGRLayer's SetSpatialFilter method.

7. Iterate over the filtered features: Iterate over the features in the layer that match the spatial filter.

Here is an example code that demonstrates these steps:

#include <stdio.h>
#include <gdal.h>
#include <ogr_api.h>
#include <ogr_srs_api.h>
#include <ogr_geometry.h>
#include <geos_c.h>

void OGRGeometryFromGEOSGeometry(const GEOSGeometry *geosGeom, OGRGeometryH *ogrGeom) {
    char *wkt = GEOSGeomToWKT(geosGeom);
    OGR_G_CreateFromWkt(&wkt, NULL, ogrGeom);
    GEOSFree(wkt);
}

int main() {
    // Initialize GDAL/OGR
    GDALAllRegister();
    OGRRegisterAll();

    // Open the data source
    const char *dataSourcePath = "path/to/your/shapefile.shp";
    GDALDatasetH hDS = GDALOpenEx(dataSourcePath, GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (hDS == NULL) {
        fprintf(stderr, "Failed to open data source.\n");
        return 1;
    }

    // Get the layer
    OGRLayerH hLayer = GDALDatasetGetLayer(hDS, 0);
    if (hLayer == NULL) {
        fprintf(stderr, "Failed to get layer.\n");
        GDALClose(hDS);
        return 1;
    }

    // Create a GEOSGeometry (example: a point)
    GEOSContextHandle_t geosContext = GEOS_init_r();
    GEOSGeometry *geosGeom = GEOSGeomFromWKT_r(geosContext, "POINT (1 1)");

    // Convert GEOSGeometry to OGRGeometry
    OGRGeometryH ogrGeom;
    OGRGeometryFromGEOSGeometry(geosGeom, &ogrGeom);

    // Set the spatial filter on the layer
    OGR_L_SetSpatialFilter(hLayer, ogrGeom);

    // Iterate over the filtered features
    OGRFeatureH hFeature;
    OGR_L_ResetReading(hLayer);
    while ((hFeature = OGR_L_GetNextFeature(hLayer)) != NULL) {
        // Process the feature
        OGR_F_DumpReadable(hFeature, stdout);

        // Destroy the feature to free memory
        OGR_F_Destroy(hFeature);
    }

    // Clean up
    OGR_G_DestroyGeometry(ogrGeom);
    GEOSGeom_destroy_r(geosContext, geosGeom);
    GEOS_finish_r(geosContext);
    GDALClose(hDS);

    return 0;
}

Explanation:

1. Initialization:

   • GDALAllRegister() and OGRRegisterAll() initialize the GDAL/OGR library.
   • GDALOpenEx() opens the data source.

2. Layer and Geometry:

   • GDALDatasetGetLayer() retrieves the layer from the data source.
   • GEOSGeomFromWKT_r() creates a GEOSGeometry from a WKT string.

3. Conversion and Filtering:

   • OGRGeometryFromGEOSGeometry() converts the GEOSGeometry to an OGRGeometry.
   • OGR_L_SetSpatialFilter() sets the spatial filter on the layer.

4. Feature Iteration:

   • OGR_L_GetNextFeature() iterates over the features that match the spatial filter.
   • OGR_F_DumpReadable() prints the feature information.

5. Cleanup:

   • Clean up allocated resources to avoid memory leaks.

Make sure to link against the GDAL and GEOS libraries when compiling this code. For example:

gcc -o filter_features filter_features.c -lgdal -lgeos_c

This example assumes you have the GDAL and GEOS libraries installed and properly configured on your