vector.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: GDALDataset wrapper for vector operations
 * Author: Joseph Meyer
 * Date: June, 2025
 *
 ******************************************************************************/

 
#pragma once
 
#include <iostream>
#include <filesystem>
 
#include <gdal_priv.h>
#include <ogrsf_frmts.h>
#include <ogr_core.h>
#include <ogr_geometry.h>
#include <ogr_api.h>
#include <ogr_recordbatch.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
 
namespace sgs {
namespace vector {
 
namespace py = pybind11;
using namespace pybind11::literals;

class GDALVectorWrapper {
    private:
    GDALDatasetUniquePtr p_dataset;
    OGRSpatialReference srs; 
    bool haveSRS = false;
 
    public:  
    GDALVectorWrapper(std::string filename, std::string projDBPath) {
        //set proj.db search path to search for the proj.db file which is included in sgs package
        char **paths = nullptr;
        paths = CSLAddString(paths, projDBPath.c_str());
        OSRSetPROJSearchPaths(paths);
        CSLDestroy(paths);
        
        //must register drivers before trying to open a dataset
        GDALAllRegister();
 
        //dataset
        this->p_dataset = GDALDatasetUniquePtr(GDALDataset::Open(filename.c_str(), GDAL_OF_VECTOR));
        if (!this->p_dataset) {
            throw std::runtime_error("dataset pointer is null after initialization, dataset unable to be initialized.");
        }
    }
  
        GDALVectorWrapper(GDALDataset *p_dataset, std::string projection) {
        this->p_dataset = GDALDatasetUniquePtr(p_dataset);
        OGRErr err = this->srs.importFromWkt(projection.c_str());
        if (err) {
            throw std::runtime_error("unable to get Spatial Reference System from projection string.");
        }
        this->haveSRS = true;
    }

    GDALVectorWrapper(std::string bytes, std::string projection, std::string name, std::string projDBPath) {
        //set proj.db search path to search for the proj.db file which is included in sgs package
        
        std::cout << "in GDALVectorWrapper constructor" << std::endl;
        char **paths = nullptr;
        paths = CSLAddString(paths, projDBPath.c_str());
        OSRSetPROJSearchPaths(paths);
        CSLDestroy(paths);
        
        GDALAllRegister();
        
        //create dataset from geojson string to copy to new layer
        GDALDataset *p_indataset = GDALDataset::FromHandle(GDALOpenEx(bytes.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
        OGRLayer *p_inlayer = p_indataset->GetLayerByName("OGRGeoJSON");
    
        //set spatial reference
            OGRErr err = this->srs.importFromWkt(projection.c_str());
        if (err) {
            throw std::runtime_error("unable to get Spatial Reference System from projection string.");
        }
        this->haveSRS = true;
 
        //create dataset and layer with correct spatial reference and 
        GDALDriver *p_driver = GetGDALDriverManager()->GetDriverByName("MEM");
        if (!p_driver) {
            throw std::runtime_error("unable to create dataset driver.");
        }
        GDALDataset *p_dataset = p_driver->Create("", 0, 0, 0, GDT_Unknown, nullptr);
        if (!p_dataset) {
            throw std::runtime_error("unable to create dataset from driver.");
        }
        OGRLayer *p_outlayer = p_dataset->CreateLayer(name.c_str(), &this->srs, wkbUnknown, nullptr);
        if (!p_outlayer) {
            throw std::runtime_error("unable to create dataset layer.");
        }
        this->p_dataset = GDALDatasetUniquePtr(p_dataset);
 
        /* copy data over from dataset created using geojson string
         *
         * This is done, because when converting a geopandas geodataframe to json, the layer name,
         * and (more importantly) the spatial reference system may not be transferred correctly.
         * This way does include A LOT of copying data, however it ensures the srs is correct and
         * all of the data from the geodataframe is moved over effectively.
         */
 
        //copy over field definitions   
        OGRFeatureDefn *p_featdef = p_inlayer->GetLayerDefn();
        int fcount = p_featdef->GetFieldCount();
        for (int i = 0; i < fcount; i++) {
            OGRFieldDefn *p_fielddef = p_featdef->GetFieldDefn(i);
            OGRErr err = p_outlayer->CreateField(p_fielddef);
            if (err) {
                std::cout << "unable to copy field definition." << std::endl;
            }
        }
 
        //copy over features
        for (auto& p_infeature : p_inlayer) {
            OGRFeature *p_outfeature = OGRFeature::CreateFeature(p_featdef);
            for (int i = 0; i < fcount; i++) {
                p_outfeature->SetField(i, p_infeature->GetRawFieldRef(i));
            }
            OGRGeometry *p_geom = p_infeature->GetGeometryRef();
            if (!p_geom) {
                throw std::runtime_error("an invalid geometry was encountered in the vector file.");                
            }
            p_geom->assignSpatialReference(&this->srs);
            p_outfeature->SetGeometry(p_geom);
            p_outlayer->CreateFeature(p_outfeature);
            OGRFeature::DestroyFeature(p_outfeature);
        }
    }   
            
    GDALDataset *getDataset() {
        return this->p_dataset.get();
    }

    std::vector<std::string> getLayerNames() {
        std::vector<std::string> retval;
 
        for (OGRLayer *p_layer : this->p_dataset->GetLayers()) {
            retval.push_back(std::string(p_layer->GetName()));
        }
 
        return retval;
    }

    std::unordered_map<std::string, std::string> getLayerInfo(std::string layerName) {
        std::unordered_map<std::string, std::string> retval;
 
        OGRLayer *p_layer = this->p_dataset->GetLayerByName(layerName.c_str());
        std::unique_ptr<OGREnvelope> extent = std::unique_ptr<OGREnvelope>(new OGREnvelope);
        p_layer->GetExtent(extent.get());
 
        retval.emplace("feature_count", std::to_string(p_layer->GetFeatureCount()));
        retval.emplace("field_count", std::to_string(p_layer->GetLayerDefn()->GetFieldCount()));
        retval.emplace("geometry_type", OGRGeometryTypeToName(p_layer->GetGeomType()));
        retval.emplace("xmin", std::to_string(extent->MinX));
        retval.emplace("xmax", std::to_string(extent->MaxX));
        retval.emplace("ymin", std::to_string(extent->MinY));
        retval.emplace("ymax", std::to_string(extent->MaxY));
        if (!p_layer->GetSpatialRef()) {
            std::cout << "WARNING: cannot get spatial reference for layer " << layerName << "." << std::endl;
        }
        else {
            retval.emplace("crs", std::string(p_layer->GetSpatialRef()->GetName()));
        }
        return retval;
    }

    OGRLayer *getLayer(std::string layerName) {
        return this->p_dataset->GetLayerByName(layerName.c_str());
    }

    std::vector<std::vector<double>> getPoints(std::string layerName) {
        OGRLayer *p_layer = this->p_dataset->GetLayerByName(layerName.c_str());
        std::vector<double> xCoords;
        std::vector<double> yCoords;
 
        for (const auto& p_feature : *p_layer) {
            OGRGeometry *p_geometry = p_feature->GetGeometryRef();
            if (!p_geometry) {
                throw std::runtime_error(std::format("an invalid geometry was encountered in the {} layer.", layerName));
            }
            switch (wkbFlatten(p_geometry->getGeometryType())) {
                case OGRwkbGeometryType::wkbPoint: {
                    OGRPoint *p_point = p_geometry->toPoint();
                    xCoords.push_back(p_point->getX());
                    yCoords.push_back(p_point->getY());
                    break;
                }
                case OGRwkbGeometryType::wkbMultiPoint: {
                    for (const auto& p_point : *p_geometry->toMultiPoint()) {
                        xCoords.push_back(p_point->getX());
                        yCoords.push_back(p_point->getY());
                    }
                    break;
                }
                default:
                    throw std::runtime_error("encountered a geometry which was not of type Point or MultiPoint.");
            }
        }
 
        return {xCoords, yCoords};
    }

    std::vector<std::string> getPointsAsWkt(std::string layerName){
        OGRLayer *p_layer = this->p_dataset->GetLayerByName(layerName.c_str());
        std::vector<std::string> retval;
 
        for (const auto& p_feature : *p_layer) {
            OGRGeometry *p_geometry = p_feature->GetGeometryRef();
            if (!p_geometry) {
                throw std::runtime_error(std::format("an invalid geometry was encountered in the {} layer.", layerName));
            }
            switch (wkbFlatten(p_geometry->getGeometryType())) {
                case OGRwkbGeometryType::wkbPoint: {
                    OGRPoint *p_point = p_geometry->toPoint();
                    retval.push_back(p_point->exportToWkt());
                    break;
                }
                case OGRwkbGeometryType::wkbMultiPoint: {
                    for (const auto& p_point : *p_geometry->toMultiPoint()) {
                        retval.push_back(p_point->exportToWkt());
                    }
                    break;
                }
                default:
                    throw std::runtime_error("encountered a geometry which was not of type Point or MultiPoint.");
            }
        }
 
        return retval;
    }

    std::vector<std::vector<std::vector<double>>> getLineStrings(std::string layerName){
        OGRLayer *p_layer = this->p_dataset->GetLayerByName(layerName.c_str());
        std::vector<std::vector<std::vector<double>>> retval;
 
        for (const auto& p_feature : *p_layer) {
            OGRGeometry *p_geometry = p_feature->GetGeometryRef();
            if (!p_geometry) {
                throw std::runtime_error(std::format("an invalid geometry was encountered in the {} layer.", layerName));
            }
            switch (wkbFlatten(p_geometry->getGeometryType())) {
                case OGRwkbGeometryType::wkbLineString: {
                    std::vector<double> xCoords;
                    std::vector<double> yCoords;
                    for (const auto& p_point : *p_geometry->toLineString()) {
                        xCoords.push_back(p_point.getX());
                        yCoords.push_back(p_point.getY());
                    }
                    retval.push_back({xCoords, yCoords});
                    break;
                }
                case OGRwkbGeometryType::wkbMultiLineString: {
                    for (const auto& p_lineString : *p_geometry->toMultiLineString()) {
                        std::vector<double> xCoords;
                        std::vector<double> yCoords;
                        for (const auto& p_point : *p_lineString) {
                            xCoords.push_back(p_point.getX());
                            yCoords.push_back(p_point.getY());
                        }
                        retval.push_back({xCoords, yCoords});
                    }
                    break;
                }
                default:
                    throw std::runtime_error("encountered a point which was not of type LineString or MultiLineString");
            }
        }
 
        return retval;
    }

    void write(std::string filename) {
        std::filesystem::path filepath = filename;
        std::string extension = filepath.extension().string();
            
        GDALAllRegister();
        GDALDriver *p_driver; 
        if (extension == ".geojson") {
            p_driver = GetGDALDriverManager()->GetDriverByName("GeoJSON");
        }
        else if (extension == ".shp") {
            p_driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
        }
        else {
            throw std::runtime_error("file extension must be one of : .geojson, .shp");
        }
    
        GDALDataset *datasetCopy = p_driver->CreateCopy(
            filename.c_str(),
            this->p_dataset.get(),
            FALSE,
            nullptr,
            nullptr,
            nullptr
        );
    
        if (!datasetCopy) {
            std::cout << "failed to create dataset with filename " << filename << "." << std::endl;
        }
    
        CPLErr err = GDALClose(datasetCopy);
    
        if (err != CE_None) {
            std::cout << "failed to close dataset of file " << filename << ". The file output may not be correct." << std::endl;
        }
    }

    std::string getFullProjectionInfo() {
        char *p_proj;
        if (haveSRS) {
            srs.exportToPrettyWkt(&p_proj);
        }
        else {
            this->p_dataset->GetLayer(0)->GetSpatialRef()->exportToPrettyWkt(&p_proj);
        }
        
        std::string retval = std::string(p_proj);
        CPLFree(p_proj);
        return retval;
    }

    OGRSpatialReference *getSRS(void) {
        if (!this->haveSRS) {
            throw std::runtime_error("do not have OGRSpatialReference associated with GDALVectorWrapper.");
        }
 
        return &this->srs;
    }
};
 
} //namespace vector
} //namespace sgs