existing.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 "raster.h"
#include "vector.h"
#include "helper.h"
 
#include <boost/unordered/unordered_flat_map.hpp>
#include <gdal_priv.h>
#include <ogrsf_frmts.h>
#include <ogr_core.h>
 
namespace sgs {
namespace existing {

struct Existing {
    bool used;
    boost::unordered::unordered_flat_map<int64_t, OGRPoint> samples;
    double IGT[6];
    int64_t width;

    Existing(
        vector::GDALVectorWrapper *p_vect,
            raster::GDALRasterWrapper *p_rast, 
        double *GT, 
        int64_t width, 
        OGRLayer *p_samples, 
        bool plot,
            std::vector<double>& xCoords,
        std::vector<double>& yCoords,
        bool addAllPoints = true) 
    {
        if (!p_vect) {
            this->used = false;
            return;
        }   
 
        OGRFieldDefn existingField("existing", OFTInteger);
        OGRErr err = p_samples->CreateField(&existingField);
        if (err) {
            throw std::runtime_error("cannot create 'existing' field in output layer.");
        }
 
        this->width = width;
 
        std::vector<std::string> layerNames = p_vect->getLayerNames();
        if (layerNames.size() > 1) {
            throw std::runtime_error("the file containing existing sample points must have only a single layer.");
        }
 
        //check to ensure spatial reference system of raster and vector match
        std::string rastProj = p_rast->getDataset()->GetProjectionRef();
        OGRSpatialReference rastSRS;
        rastSRS.importFromWkt(rastProj.c_str());
        const OGRSpatialReference *p_sampSRS = p_samples->GetSpatialRef();
        if (!rastSRS.IsSame(p_sampSRS)) {
            throw std::runtime_error("existing sample vector and raster do not have the same spatial reference system.");   
        }
 
        //invert geotransform so we can use IGT to convert from point to indexes
        GDALInvGeoTransform(GT, this->IGT);
 
        std::string name = layerNames[0];
        OGRLayer *p_layer = p_vect->getLayer(name);
        helper::Field fieldExistingTrue("existing", 1);
 
        for (const auto& p_feature : *p_layer) {
            OGRGeometry *p_geometry = p_feature->GetGeometryRef();
            switch (wkbFlatten(p_geometry->getGeometryType())) {
                case OGRwkbGeometryType::wkbPoint: {
                    OGRPoint *p_point = p_geometry->toPoint();
                    int64_t index = helper::point2index<int64_t>(p_point->getX(), p_point->getY(), IGT, width);
                    this->samples.emplace(index, *p_point);
                    if (addAllPoints) {
                        helper::addPoint(p_point, p_samples, &fieldExistingTrue);
                        if (plot) {
                            xCoords.push_back(p_point->getX());
                            yCoords.push_back(p_point->getY());
                        }
                    }
                    break;
                }
                case OGRwkbGeometryType::wkbMultiPoint: {
                    for (const auto& p_point : *p_geometry->toMultiPoint()) {
                        int64_t index = helper::point2index<int64_t>(p_point->getX(), p_point->getY(), IGT, width);
                        this->samples.emplace(index, *p_point);
                        if (addAllPoints) {
                            helper::addPoint(p_point, p_samples, &fieldExistingTrue);
                            if (plot) {
                                xCoords.push_back(p_point->getX());
                                yCoords.push_back(p_point->getY());
                            }
                        }
                    }
                    break;
                }
                default:
                    throw std::runtime_error("the file containing existing sample points must have only Point or MultiPoint geometries.");
            }       
        }
 
        this->used = true;
    }

    inline bool
    containsIndex(int64_t x, int64_t y) {
        int64_t index = y * this->width + x;
        return this->samples.find(index) != this->samples.end();    
    }

    inline OGRPoint
    getPoint(int64_t x, int64_t y) {
        int64_t index = y * this->width + x;
        return this->samples.find(index)->second;
    }

    inline bool
    containsCoordinates(double xCoord, double yCoord) {
        int64_t index = helper::point2index<int64_t>(xCoord, yCoord, this->IGT, this->width);
        return this->samples.find(index) != this->samples.end();
    }

    inline size_t
    count() {
        return this->samples.size();
    }
};
 
} //namespace existing
} //namespace sgs