access.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: generate an access mask using vector and raster datasets
 * Author: Joseph Meyer
 * Date: October, 2025
 *
 ******************************************************************************/

 
#pragma once
 
#include "helper.h"
#include "raster.h"
#include "vector.h"
#include "gdal_utils.h"
 
namespace sgs {
namespace access {

struct Access {
    bool used = false;
    double area = -1;
    GDALDataset *p_dataset = nullptr;
    helper::RasterBandMetaData band;

    Access(vector::GDALVectorWrapper *p_vector,
           raster::GDALRasterWrapper *p_raster,
           std::string layerName, 
           double buffInner, 
           double buffOuter,
           bool largeRaster,
           std::string tempFolder,
           int xBlockSize,
           int yBlockSize) 
    {
        if (!p_vector) {
            return;
        }
 
        std::string rastProj = p_raster->getDataset()->GetProjectionRef();
        OGRSpatialReference rastSRS;
        rastSRS.importFromWkt(rastProj.c_str());
        OGRLayer *p_inputLayer = p_vector->getLayer(layerName.c_str());
        const OGRSpatialReference *p_vectSRS = p_inputLayer->GetSpatialRef();
        if (!rastSRS.IsSame(p_vectSRS)) {
            throw std::runtime_error("access vector and raster do not have the same spatial reference system.");
        }
 
        //step 1: create multipolygon buffers
        OGRMultiPolygon *buffInnerPolygons = new OGRMultiPolygon;
        OGRMultiPolygon *buffOuterPolygons = new OGRMultiPolygon;
        OGRGeometry *p_polygonMask;
    
        //occasionally, samples end up being placed just outside the accessible area,
        //typically when the buffer sizes are multiples of the pixel size. 
        //
        //Adjusting the buffers minorly fixes this problem.
        double pixelSize = std::min(p_raster->getPixelHeight(), p_raster->getPixelWidth());
        buffOuter = buffOuter - pixelSize / 50;
        buffInner = buffInner == 0 ? 0 : buffInner + pixelSize / 50;
 
        //step 2: add geometries to access polygon buffers
        for (const auto& p_feature : *p_vector->getLayer(layerName.c_str())) {
            OGRGeometry *p_geometry = p_feature->GetGeometryRef();
            if (!p_geometry) {
                throw std::runtime_error(std::format("an invalid geometry was encountered in the {} layer.", layerName));
            }
 
            OGRwkbGeometryType type = wkbFlatten(p_geometry->getGeometryType());
    
            switch (type) {
                case OGRwkbGeometryType::wkbLineString: {
                    OGRGeometry *p_outer = p_geometry->Buffer(buffOuter);
                    buffOuterPolygons->addGeometry(p_outer);
                    OGRGeometryFactory::destroyGeometry(p_outer);   
    
                    if (buffInner != 0) {
                        OGRGeometry *p_inner = p_geometry->Buffer(buffInner);
                        buffInnerPolygons->addGeometry(p_inner);
                        OGRGeometryFactory::destroyGeometry(p_inner);
                    }
                    break;
                }
                case OGRwkbGeometryType::wkbMultiLineString: {
                    for (const auto& p_lineString : *p_geometry->toMultiLineString()) {
                        OGRGeometry *p_outer = p_lineString->Buffer(buffOuter);
                        buffOuterPolygons->addGeometry(p_outer);
                        OGRGeometryFactory::destroyGeometry(p_outer);   
    
                        if (buffInner != 0) {
                            OGRGeometry *p_inner = p_lineString->Buffer(buffInner);
                            buffInnerPolygons->addGeometry(p_inner);
                            OGRGeometryFactory::destroyGeometry(p_inner);
                        }
                    }
    
                    break;
                }
                default: {
                    throw std::runtime_error("geometry type must be LineString or MultiLineString");
                }
            }
        }   
    
        //step 3: generate the polygon mask and free no longer used memory
        if (buffInner == 0) {
            p_polygonMask = buffOuterPolygons->UnionCascaded();
        }
        else {
            OGRGeometry *buffOuterUnion = buffOuterPolygons->UnionCascaded();
            OGRGeometry *buffInnerUnion = buffInnerPolygons->UnionCascaded();
            p_polygonMask = buffOuterUnion->Difference(buffInnerUnion);
            OGRGeometryFactory::destroyGeometry(buffOuterUnion);
            OGRGeometryFactory::destroyGeometry(buffInnerUnion);    
        }
        delete buffOuterPolygons;
        delete buffInnerPolygons;
 
        //update area to contain the total accessible are within the raster 
        this->area = 0;
        OGRPolygon extent(p_raster->getXMin(), p_raster->getYMin(), p_raster->getXMax(), p_raster->getYMax());
        OGRGeometry *p_polygonWithinExtent;
        OGRwkbGeometryType geomtype = p_polygonMask->getGeometryType();
        if (geomtype == wkbPolygon) {
            p_polygonWithinExtent = p_polygonMask->toPolygon()->Intersection(&extent);
        }
        else if (geomtype == wkbMultiPolygon) {
            p_polygonWithinExtent = p_polygonMask->toMultiPolygon()->Intersection(&extent);
        }
        else {
            throw std::runtime_error("p_polygonMask is not a polygon or a multipolygon!");
        }
 
        geomtype = p_polygonWithinExtent->getGeometryType();
        if (geomtype == wkbPolygon) {
            this->area += p_polygonWithinExtent->toPolygon()->toUpperClass()->get_Area();
        }
        else if (geomtype == wkbMultiPolygon) {
            for (const auto& p_polygon : *p_polygonWithinExtent->toMultiPolygon()) {
                this->area += p_polygon->toUpperClass()->get_Area();
            }
        }
        else {
            throw std::runtime_error("taking the intersection of polygons somehow resulted in something not a polygon! This is a bug.");
        }
 
        //invert polygon mask
        OGRGeometry *p_invertedMask = extent.Difference(p_polygonMask);
 
        //step 4: create new GDAL dataset to rasterize as access mask
        GDALDataset *p_accessPolygonDataset = GetGDALDriverManager()->GetDriverByName("MEM")->Create(
            "",
            0,
            0,
            0,
            GDT_Unknown,
            nullptr
        );
    
        OGRLayer *p_layer = p_accessPolygonDataset->CreateLayer(
            "access", 
            p_vector->getLayer(layerName.c_str())->GetSpatialRef(), 
            p_invertedMask->getGeometryType(), 
            nullptr
        );
        OGRFieldDefn field("index", OFTInteger);
        p_layer->CreateField(&field);
        
        //step 6: add the access polygon to the new layer
        OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
        p_feature->SetField("index", 0);
        p_feature->SetGeometry(p_invertedMask);
        p_layer->CreateFeature(p_feature); //error handling here???
        OGRFeature::DestroyFeature(p_feature);
    
        std::filesystem::path path = tempFolder;
        path = path / "access.tif";
 
        //step 9: generate options list for rasterization   
        char **argv = nullptr;
 
        //specify 'all touched'
        argv = CSLAddString(argv, "-at");
 
        //specify the burn value for the polygon
        argv = CSLAddString(argv, "-burn");
        argv = CSLAddString(argv, std::to_string(1).c_str());
    
        //specify the layer
        argv = CSLAddString(argv, "-l");
        argv = CSLAddString(argv, "access");
 
        //specify extent
        argv = CSLAddString(argv, "-te");
        argv = CSLAddString(argv, std::to_string(p_raster->getXMin()).c_str());
        argv = CSLAddString(argv, std::to_string(p_raster->getYMin()).c_str());
        argv = CSLAddString(argv, std::to_string(p_raster->getXMax()).c_str());
        argv = CSLAddString(argv, std::to_string(p_raster->getYMax()).c_str());
    
        //specify dimensions
        argv = CSLAddString(argv, "-ts");
        argv = CSLAddString(argv, std::to_string(p_raster->getWidth()).c_str());
        argv = CSLAddString(argv, std::to_string(p_raster->getHeight()).c_str());
 
        //specify output type
        argv = CSLAddString(argv, "-ot");
        argv = CSLAddString(argv, "Int8");
 
        GDALRasterizeOptions *options = GDALRasterizeOptionsNew(argv, nullptr);
 
        //step 10: rasterize vector creating in-memory dataset
        this->p_dataset = GDALDataset::FromHandle(GDALRasterize(
            path.string().c_str(),
            nullptr,
            p_accessPolygonDataset,
            options,
            nullptr
        ));
 
        this->band.p_band = this->p_dataset->GetRasterBand(1);
        this->band.type = GDT_Int8;
        this->band.size = sizeof(int8_t);
        this->band.p_band->GetBlockSize(&this->band.xBlockSize, &this->band.yBlockSize);
 
        //step 11: free dynamically allocated data
        GDALRasterizeOptionsFree(options);
        CSLDestroy(argv);
        GDALClose(GDALDataset::ToHandle(p_accessPolygonDataset));
        OGRGeometryFactory::destroyGeometry(p_polygonMask);
        OGRGeometryFactory::destroyGeometry(p_polygonWithinExtent);
        OGRGeometryFactory::destroyGeometry(p_invertedMask);
 
        this->used = true;
    }
 
    ~Access() {
        GDALClose(GDALDataset::ToHandle(this->p_dataset));
    }
};
 
} //namespace access
} //namespace sgs