Intersecting OGRGeometryH with Shapefile Layer in C

In C, get the intersection of a OGRGeometryH and the shapefile layer

To get the intersection of an OGRGeometryH and a shapefile layer in C, you can use the GDAL/OGR library. Below is a step-by-step guide and example code to achieve this:

1. Include the necessary headers:

   #include "ogr_api.h"
   #include "ogr_srs_api.h"
   #include "cpl_conv.h" // for CPLMalloc()
   

2. Initialize the GDAL/OGR library:

   OGRRegisterAll();
   

3. Open the shapefile:

   const char *shapefilePath = "path/to/your/shapefile.shp";
   OGRDataSourceH hDS = OGROpen(shapefilePath, FALSE, NULL);
   if (hDS == NULL) {
       fprintf(stderr, "Failed to open shapefile.\n");
       exit(1);
   }
   

4. Get the layer from the shapefile:

   OGRLayerH hLayer = OGR_DS_GetLayer(hDS, 0);
   if (hLayer == NULL) {
       fprintf(stderr, "Failed to get layer from shapefile.\n");
       OGR_DS_Destroy(hDS);
       exit(1);
   }
   

5. Create the OGRGeometryH object:

   // Example: Create a point geometry
   OGRGeometryH hGeom;
   OGR_G_CreateFromWkt("POINT (1 1)", NULL, &hGeom);
   

6. Iterate through the features in the layer and compute the intersection:

   OGRFeatureH hFeature;
   OGR_L_ResetReading(hLayer);
   while ((hFeature = OGR_L_GetNextFeature(hLayer)) != NULL) {
       OGRGeometryH hFeatureGeom = OGR_F_GetGeometryRef(hFeature);
       if (hFeatureGeom != NULL) {
           OGRGeometryH hIntersection = OGR_G_Intersection(hGeom, hFeatureGeom);
           if (hIntersection != NULL && !OGR_G_IsEmpty(hIntersection)) {
               // Do something with the intersection geometry
               char *wkt;
               OGR_G_ExportToWkt(hIntersection, &wkt);
               printf("Intersection: %s\n", wkt);
               CPLFree(wkt);
               OGR_G_DestroyGeometry(hIntersection);
           }
       }
       OGR_F_Destroy(hFeature);
   }
   

7. Clean up:

   OGR_G_DestroyGeometry(hGeom);
   OGR_DS_Destroy(hDS);
   

8. Unregister GDAL/OGR (optional):

   OGRCleanupAll();
   

Here is the complete example code:

#include "ogr_api.h"
#include "ogr_srs_api.h"
#include "cpl_conv.h" // for CPLMalloc()

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

    // Open the shapefile
    const char *shapefilePath = "path/to/your/shapefile.shp";
    OGRDataSourceH hDS = OGROpen(shapefilePath, FALSE, NULL);
    if (hDS == NULL) {
        fprintf(stderr, "Failed to open shapefile.\n");
        return 1;
    }

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

    // Create the OGRGeometryH object (example: a point)
    OGRGeometryH hGeom;
    OGR_G_CreateFromWkt("POINT (1 1)", NULL, &hGeom);

    // Iterate through the features in the layer and compute the intersection
    OGRFeatureH hFeature;
    OGR_L_ResetReading(hLayer);
    while ((hFeature = OGR_L_GetNextFeature(hLayer)) != NULL) {
        OGRGeometryH hFeatureGeom = OGR_F_GetGeometryRef(hFeature);
        if (hFeatureGeom != NULL) {
            OGRGeometryH hIntersection = OGR_G_Intersection(hGeom, hFeatureGeom);
            if (hIntersection != NULL && !OGR_G_IsEmpty(hIntersection)) {
                // Do something with the intersection geometry
                char *wkt;
                OGR_G_ExportToWkt(hIntersection, &wkt);
                printf("Intersection: %s\n", wkt);
                CPLFree(wkt);
                OGR_G_DestroyGeometry(hIntersection);
            }
        }
        OGR_F_Destroy(hFeature);
    }

    // Clean up
    OGR_G_DestroyGeometry(hGeom);
    OGR_DS_Destroy(hDS);

    // Unregister GDAL/OGR (optional)
    OGRCleanupAll();

    return 0;
}

Make sure to link against the GDAL library when compiling your program. For example:

gcc -o intersect_example intersect_example.c -lgdal

This example assumes you have a point geometry and a shapefile with geometries to intersect. Adjust the OGR_G_CreateFromWkt function to create the appropriate geometry type as needed.