Raster.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <array>
#include <stdexcept>
#include <utility>
 
#include "BaseLib/Error.h"
#include "MathLib/Point3d.h"
 
namespace GeoLib
{

struct RasterHeader final
{
    std::size_t n_cols;       // width
    std::size_t n_rows;       // height
    std::size_t n_depth;      // depth (for 3d image)
    MathLib::Point3d origin;  // lower left corner
    double cell_size;         // edge length of each pixel
    double no_data;           // no data value
 
    bool operator==(RasterHeader const&) const = default;
};

class Raster final
{
public:
    template <typename InputIterator>
    Raster(RasterHeader header, InputIterator begin, InputIterator end)
        : _header(std::move(header)),
          _raster_data(_header.n_cols * _header.n_rows)
    {
        unsigned long const number_of_input_values =
            static_cast<unsigned long>(std::distance(begin, end));
        if (number_of_input_values != _header.n_cols * _header.n_rows)
        {
            throw std::out_of_range(
                "Number of raster data mismatch, need " +
                std::to_string(_header.n_cols * _header.n_rows) +
                " values, but got " + std::to_string(number_of_input_values));
        }
        std::copy(begin, end, _raster_data.begin());
    }
 
    Raster(Raster const&) = default;
    Raster(Raster&&) = default;
    Raster& operator=(Raster const&) = default;
    Raster& operator=(Raster&&) = default;
    ~Raster() = default;

    RasterHeader const& getHeader() const { return _header; }

    void refineRaster(std::size_t scaling);

    std::vector<double>::const_iterator begin() const
    {
        return _raster_data.begin();
    }

    std::vector<double>::const_iterator end() const
    {
        return _raster_data.end();
    }
 
    double const* data() const { return _raster_data.data(); }

    double const& operator()(std::size_t const row, std::size_t const col) const
    {
        if (row >= _header.n_rows || col >= _header.n_cols)
        {
            OGS_FATAL(
                "Raster pixel ({}, {}) doesn't exist. Raster size is {} x {}.",
                row, col, _header.n_rows, _header.n_cols);
        }
        return _raster_data[(_header.n_rows - 1 - row) * _header.n_cols + col];
    }
    double& operator()(std::size_t const row, std::size_t const col)
    {
        return const_cast<double&>(std::as_const(*this)(row, col));
    }

    double getValueAtPoint(const MathLib::Point3d& pnt) const;

    double interpolateValueAtPoint(const MathLib::Point3d& pnt) const;

    bool isPntOnRaster(MathLib::Point3d const& pnt) const;
 
    bool operator==(Raster const&) const = default;
 
private:
    void setCellSize(double cell_size);
    void setNoDataVal(double no_data_val);
 
    GeoLib::RasterHeader _header;
    std::vector<double> _raster_data;
};
 
struct NamedRaster
{
    std::string raster_name;
    std::unique_ptr<GeoLib::Raster> raster;
};
 
}  // namespace GeoLib