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:
Include the necessary headers:
#include "ogr_api.h"
#include "ogr_srs_api.h"
#include "cpl_conv.h" // for CPLMalloc()
Initialize the GDAL/OGR library:
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");
exit(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);
exit(1);
}
Create the OGRGeometryH
object:
// Example: Create a point geometry
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();
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.