Reading Shapefiles in C with GDAL and GEOS Libraries

In C, read shapefile using GEOS library

The GEOS (Geometry Engine - Open Source) library is primarily used for geometric operations and does not directly handle shapefiles. To read shapefiles in C, you typically use the Shapefile C Library (shapelib) or the GDAL (Geospatial Data Abstraction Library). Once you have the geometries from the shapefile, you can use GEOS for geometric operations.

Here's a basic example of how to read a shapefile using the GDAL library and then use GEOS for some geometric operations:

1. Install GDAL and GEOS: Make sure you have both GDAL and GEOS installed on your system. You can usually install them via package managers like apt on Ubuntu or brew on macOS.

2. Include the necessary headers: Include the headers for GDAL and GEOS in your C program.

3. Read the shapefile using GDAL: Use GDAL to read the shapefile and extract geometries.

4. Perform geometric operations using GEOS: Use GEOS to perform geometric operations on the extracted geometries.

Here's a simple example to illustrate these steps:

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

void geos_message_handler(const char *message, void *userdata) {
    fprintf(stderr, "GEOS: %s\n", message);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "Usage: %s <shapefile>\n", argv[0]);
        return 1;
    }

    const char *shapefile = argv[1];

    // Initialize GDAL
    GDALAllRegister();

    // Open the shapefile
    GDALDatasetH hDS = GDALOpenEx(shapefile, GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (hDS == NULL) {
        fprintf(stderr, "Failed to open shapefile: %s\n", shapefile);
        return 1;
    }

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

    // Initialize GEOS
    initGEOS(geos_message_handler, geos_message_handler);

    // Iterate over features in the layer
    OGRFeatureH hFeature;
    OGR_L_ResetReading(hLayer);
    while ((hFeature = OGR_L_GetNextFeature(hLayer)) != NULL) {
        OGRGeometryH hGeometry = OGR_F_GetGeometryRef(hFeature);
        if (hGeometry != NULL) {
            char *wkt;
            OGR_G_ExportToWkt(hGeometry, &wkt);

            // Create GEOS geometry from WKT
            GEOSWKTReader *reader = GEOSWKTReader_create();
            GEOSGeometry *geosGeom = GEOSWKTReader_read(reader, wkt);
            if (geosGeom != NULL) {
                // Perform some GEOS operations
                double area = GEOSArea(geosGeom);
                printf("Geometry area: %f\n", area);

                // Destroy GEOS geometry
                GEOSGeom_destroy(geosGeom);
            } else {
                fprintf(stderr, "Failed to create GEOS geometry from WKT\n");
            }

            // Clean up
            GEOSWKTReader_destroy(reader);
            CPLFree(wkt);
        }

        OGR_F_Destroy(hFeature);
    }

    // Clean up
    GDALClose(hDS);
    finishGEOS();

    return 0;
}

Explanation:

1. Initialization:

   • GDAL is initialized using GDALAllRegister().
   • GEOS is initialized using initGEOS().

2. Reading the Shapefile:

   • The shapefile is opened using GDALOpenEx().
   • The first layer is retrieved using GDALDatasetGetLayer().

3. Iterating Over Features:

   • Features are read in a loop using OGR_L_GetNextFeature().
   • For each feature, the geometry is extracted using OGR_F_GetGeometryRef().

4. Using GEOS:

   • The geometry is converted to WKT (Well-Known Text) using OGR_G_ExportToWkt().
   • A GEOS geometry is created from the WKT using GEOSWKTReader_read().
   • A geometric operation (calculating the area) is performed using GEOSArea().
   • The GEOS geometry is destroyed using GEOSGeom_destroy().

5. Cleanup:

   • The GDAL dataset is closed using GDALClose().
   • GEOS is finalized using finishGEOS().

Compilation:

To compile this program, you need to link against the GDAL and GEOS libraries. For example:

gcc -o read_shapefile read_shapefile.c -lgdal -lgeos_c

Make sure the GDAL and GEOS libraries are installed and accessible to your compiler.