Intersecting OGRGeometryH with Shapefile Layer in C

Answered using GPT-4o by GPT4oMini.app

Question:

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

Answer:

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.