systematic.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: C++ implementation of systematic sampling
 * Author: Joseph Meyer
 * Date: July, 2025
 *
 ******************************************************************************/

 
#include <iostream>
#include <random>
#include <optional>
 
#include "utils/access.h"
#include "utils/existing.h"
#include "utils/helper.h"
#include "utils/raster.h"
#include "utils/vector.h"
 
namespace sgs {
namespace systematic {

OGRGeometry *getAccessPolygon(vector::GDALVectorWrapper *p_access, std::string layerName, double buffInner, double buffOuter) {
    //step 1: create multipolygon buffers
    OGRMultiPolygon *buffInnerPolygons = new OGRMultiPolygon;
    OGRMultiPolygon *buffOuterPolygons = new OGRMultiPolygon;
    OGRGeometry *p_polygonMask;
 
    //step 2: add geometries to access polygon buffers
    for (const auto& p_feature : *p_access->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_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;
            }
            default: {
                throw std::runtime_error("access polygon 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;
 
    return p_polygonMask;
}

inline bool
checkExtent(double x, double y, double xMin, double xMax, double yMin, double yMax) {
    return (x >= xMin && x <= xMax && y >= yMin && y <= yMax); 
}

inline bool
checkAccess(OGRPoint *p_point, OGRGeometry *p_geometry) {
    return !p_geometry || p_point->Within(p_geometry);
}

inline bool
checkExisting(double x, double y, existing::Existing& existing) {
    return !existing.used || !existing.containsCoordinates(x, y);
}

inline bool
checkNotNan(raster::GDALRasterWrapper *p_raster, double *IGT, double xCoord, double yCoord, bool force) {
    if (force) {
        int x = static_cast<int>(IGT[0] + xCoord * IGT[1] + yCoord * IGT[2]);
        int y = static_cast<int>(IGT[3] + xCoord * IGT[4] + yCoord * IGT[5]);
    
        for (int i = 0; i < p_raster->getBandCount(); i++) {
            GDALRasterBand *p_band = p_raster->getRasterBand(i);
            
            double val;
            OGRErr err = p_band->RasterIO(GF_Read, x, y, 1, 1, &val, 1, 1, GDT_Float64, 0, 0);
            if (err) {
                throw std::runtime_error("could not read value from raster.");
            }
 
            if (val == p_band->GetNoDataValue() || std::isnan(val)) {
                return false;
            }   
        }   
    }
 
    return true;
}   

std::tuple<
    vector::GDALVectorWrapper *, //GDALDataset containing sample points
    std::vector<std::vector<double>>, //array of samples to plot
    std::vector<std::vector<std::vector<double>>> //array of grid to plot
>
systematic(
    raster::GDALRasterWrapper *p_raster,
    double cellSize,
    std::string shape,
    std::string location,
    vector::GDALVectorWrapper *p_existing,
    vector::GDALVectorWrapper *p_access,
    std::string layerName,
    double buffInner,
    double buffOuter,
    bool force,
    bool plot,
    std::optional<uint64_t> random_state,
    std::string filename)
{
    GDALAllRegister();
    
    double *GT;     //geotransform
    double IGT[6];  //inverse geotransform
    GT = p_raster->getGeotransform();
    GDALInvGeoTransform(GT, IGT);   
 
    //determine raster extent
    double xMin, xMax, yMin, yMax;
    xMin = p_raster->getXMin();
    xMax = p_raster->getXMax();
    yMin = p_raster->getYMin();
    yMax = p_raster->getYMax();
    
    //generate random number generator
    double xDiff = xMax - xMin;
    double yDiff = yMax - yMin;
    double rngMax = yDiff * xDiff;
    std::mt19937 engine(random_state.value_or(std::random_device{}()));
    std::uniform_real_distribution<double> dist(0.0, rngMax);
    auto rng = [&] { return dist(engine); };
 
    //determine random rotation angle within extent polygon
    double rotation = rng() / (rngMax / 180);
 
    //determine grid creation function
    std::string gridFunction;
    if (shape == "square") {
        gridFunction = "ST_SquareGrid";
    }
    else { //shape == "hexagon" 
        gridFunction = "ST_HexagonalGrid";
    }
 
    //create sql query using extent polygon and grid function
    std::string extentPolygon = "'POLYGON (( " 
        + std::to_string(xMin) + " " + std::to_string(yMin) + ", "
        + std::to_string(xMin) + " " + std::to_string(yMax) + ", "
        + std::to_string(xMax) + " " + std::to_string(yMax) + ", "
        + std::to_string(xMax) + " " + std::to_string(yMin) + ", "
        + std::to_string(xMin) + " " + std::to_string(yMin) + " ))'";
 
    std::string queryString = "SELECT RotateCoords(" 
        + gridFunction 
        + "(RotateCoords("
        + "ST_GeomFromText("
        + extentPolygon
        + "), " + std::to_string(-rotation) + "), "
        + std::to_string(cellSize) + "), "
        + std::to_string(rotation) + ")";
 
    //query to create grid
    OGRLayer *p_grid = p_raster->getDataset()->ExecuteSQL(queryString.c_str(), nullptr, "SQLITE");
 
    //create output dataset before anything that might take a long time (in case creation of dataset fails)
    GDALDriver *p_driver = GetGDALDriverManager()->GetDriverByName("MEM");
    if(!p_driver) {
        throw std::runtime_error("unable to create output sample dataset driver.");
    }
    GDALDataset *p_sampleDataset = p_driver->Create("", 0, 0, 0, GDT_Unknown, nullptr);
    if (!p_sampleDataset) {
        throw std::runtime_error("unable to create output dataset with driver.");
    }
 
    vector::GDALVectorWrapper *p_wrapper = new vector::GDALVectorWrapper(p_sampleDataset, std::string(p_raster->getDataset()->GetProjectionRef()));
    OGRLayer *p_sampleLayer = p_sampleDataset->CreateLayer("samples", p_wrapper->getSRS(), wkbPoint, nullptr);    
    if (!p_sampleLayer) {
        throw std::runtime_error("unable to create output dataset layer.");
    }
 
    //get access polygon if access is given
    OGRGeometry *access = nullptr;
    if (p_access) {
        //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;
 
        access = getAccessPolygon(p_access, layerName, buffInner, buffOuter);
    }
 
    //coordinate representations the samples as vectors only if PLOT is true, returned to the (Python) caller
    std::vector<double> xCoords, yCoords;
 
    //create existing struct
    existing::Existing existing(p_existing, p_raster, GT, p_raster->getWidth(), p_sampleLayer, plot, xCoords, yCoords);
    helper::Field fieldExistingFalse("existing", 0);
 
    //grid represents the grid used to create the sampels only if PLOT is true, and is returned to the (Python) caller
    std::vector<std::vector<std::vector<double>>> grid;
 
    //iterate through the polygons in the grid and populate the ruturn data depending on user inputs
    for (const auto& p_feature : *p_grid) {
        OGRGeometry *p_geometry = p_feature->GetGeometryRef();
        if (!p_geometry) {
            throw std::runtime_error("An invalid geometry was encountered in the grid. This is a bug and should be reported to the sgspy github page.");
        }
        for (const auto& p_polygon : *p_geometry->toMultiPolygon()) {
            //generate sample depending on 'location' parameter
            OGRPoint point;
            OGRPoint secondPoint;
            if (location == "centers") {
                p_polygon->Centroid(&point);
 
                double x = point.getX();
                double y = point.getY();
                if (checkExtent(x, y, xMin, xMax, yMin, yMax) &&
                    checkAccess(&point, access) &&
                    checkExisting(x, y, existing) &&
                    checkNotNan(p_raster, IGT, x, y, force)) 
                {
                    existing.used ?
                        helper::addPoint(&point, p_sampleLayer, &fieldExistingFalse) :
                        helper::addPoint(&point, p_sampleLayer);
 
                    if (plot) {
                        xCoords.push_back(x);
                        yCoords.push_back(y);
                    }   
                }
            }
            else if (location == "corners") {
                (*p_polygon->begin())->getPoint(0, &point);
                (*p_polygon->begin())->getPoint(1, &secondPoint);
 
                double x = point.getX();
                double y = point.getY();
                if (checkExtent(x, y, xMin, xMax, yMin, yMax) &&
                    checkAccess(&point, access) &&
                    checkExisting(x, y, existing) &&
                    checkNotNan(p_raster, IGT, x, y, force)) 
                {
                    existing.used ?
                        helper::addPoint(&point, p_sampleLayer, &fieldExistingFalse) :
                        helper::addPoint(&point, p_sampleLayer);
 
                    if (plot) {
                        xCoords.push_back(x);
                        yCoords.push_back(y);
                    }   
                }
 
                x = secondPoint.getX();
                y = secondPoint.getY();
                if (checkExtent(x, y, xMin, xMax, yMin, yMax) &&
                    checkAccess(&secondPoint, access) &&
                    checkExisting(x, y, existing) &&
                    checkNotNan(p_raster, IGT, x, y, force)) 
                {
                    existing.used ?
                        helper::addPoint(&point, p_sampleLayer, &fieldExistingFalse) :
                        helper::addPoint(&secondPoint, p_sampleLayer);
 
                    if (plot) {
                        xCoords.push_back(x);
                        yCoords.push_back(y);
                    }   
                }
 
            }
            else { //location == "random"
                OGREnvelope envelope;
                p_polygon->getEnvelope(&envelope);
                double xMinEnv = envelope.MinX;
                double xMaxEnv = envelope.MaxX;
                double xDiffEnv = xMaxEnv - xMinEnv;
                double yMinEnv = envelope.MinY;
                double yMaxEnv = envelope.MaxY;
                double yDiffEnv = yMaxEnv - yMinEnv;
 
                int tries = 0;
                bool found = false;
 
                while (tries < 10 && !found) {
                    point.setX(xMinEnv + rng() / (rngMax / xDiffEnv));
                    point.setY(yMinEnv + rng() / (rngMax / yDiffEnv));
                    while (!p_polygon->Contains(&point)) {
                        point.setX(xMinEnv + rng() / (rngMax / xDiffEnv));
                        point.setY(yMinEnv + rng() / (rngMax / yDiffEnv));
                    }
 
                    double x = point.getX();
                    double y = point.getY();
                    if (checkExtent(x, y, xMin, xMax, yMin, yMax) &&
                            checkAccess(&point, access) &&
                            checkExisting(x, y, existing) &&
                        checkNotNan(p_raster, IGT, x, y, force)) 
                    {
                        found = true;
                        existing.used ?
                            helper::addPoint(&point, p_sampleLayer, &fieldExistingFalse) :
                            helper::addPoint(&point, p_sampleLayer);
 
                        if (plot) {
                            xCoords.push_back(x);
                            yCoords.push_back(y);
                        }   
                    }
 
                    tries++;
                }
            }
    
            //set grid vector to be plot
            if (plot) {
                grid.push_back({}); //add new polygon to grid
                grid.back().push_back({}); //add new x vector to polygon grid
                grid.back().push_back({}); //add new y vector to polygon grid
                for (const auto& p_linearRing : *p_polygon) {
                    for (const auto& p_point : *p_linearRing) {
                        grid.back()[0].push_back(p_point.getX());
                        grid.back()[1].push_back(p_point.getY());
                    }
                }
            }
        }
    }
 
    if (filename != "") {
        try {
            p_wrapper->write(filename);
        }
        catch (const std::exception& e) {
            std::cout << "Exception thrown trying to write file: " << e.what() << std::endl;
        }
    }   
 
    //clean up memory
    p_raster->getDataset()->ReleaseResultSet(p_grid);
    if (p_access) {
        OGRGeometryFactory::destroyGeometry(access);
    }
 
    return {p_wrapper, {xCoords, yCoords}, grid};
}
 
} //namespace systematic
} //namespace sgs