diff --git a/include/forcing/ForcingsEngineGriddedDataProvider.hpp b/include/forcing/ForcingsEngineGriddedDataProvider.hpp new file mode 100644 index 0000000000..a1e9ae5d41 --- /dev/null +++ b/include/forcing/ForcingsEngineGriddedDataProvider.hpp @@ -0,0 +1,107 @@ +#pragma once + +#include + +#if NGEN_WITH_PYTHON + +#include +#include +#include + +namespace data_access { + +struct GridSpecification +{ + //! Total number of rows (aka y) + std::uint64_t rows; + + //! Total number of columns (aka x) + std::uint64_t columns; + + //! Extent of the grid region (aka min-max corner points) + geojson::box_t extent; +}; + +struct GridMask { + //! Geographic extent of the mask region (aka min-max corner points) + geojson::box_t extent; + + //! Gridded extent of the mask region (aka min-max corner indices) + std::uint64_t rmin; + std::uint64_t cmin; + std::uint64_t rmax; + std::uint64_t cmax; + + constexpr std::uint64_t rows() const noexcept + { + return rmax - rmin + 1; + } + + constexpr std::uint64_t columns() const noexcept + { + return cmax - cmin + 1; + } + + constexpr std::uint64_t size() const noexcept + { + return rows() * columns(); + } + + static constexpr auto bad_position = static_cast(-1); + + //! Derive a discrete position along a real axis from a coordinate component + //! @param component Coordinate component value + //! @param lower_bound Lower bound of component's axis + //! @param upper_bound Upper bound of component's axis + //! @param cardinality Total number of elements along the discrete axis + //! @returns A position along the [0, cardinality) discrete axis. + static std::uint64_t position(double component, double lower_bound, double upper_bound, std::uint64_t cardinality) + { + if (component < lower_bound || component > upper_bound) { + return bad_position; + } + + // We can use a static_cast instead of std::floor because the position will never be negative, + // so truncating and flooring are equivalent here. + return static_cast( + (component - lower_bound) * (cardinality / (upper_bound - lower_bound)) + ); + } +}; + + +struct ForcingsEngineGriddedDataProvider final : + public ForcingsEngineDataProvider +{ + using base_type = ForcingsEngineDataProvider; + + ~ForcingsEngineGriddedDataProvider() override = default; + + ForcingsEngineGriddedDataProvider( + const std::string& init, + std::size_t time_begin_seconds, + std::size_t time_end_seconds, + geojson::box_t mask + ); + + data_type get_value( + const selection_type& selector, + data_access::ReSampleMethod m + ) override; + + std::vector get_values( + const selection_type& selector, + data_access::ReSampleMethod m + ) override; + + const GridMask& mask(); + + private: + int var_grid_id_; + GridSpecification var_grid_; + GridMask var_grid_mask_; +}; + +} // namespace data_access + +#endif // NGEN_WITH_PYTHON diff --git a/include/forcing/GridDataSelector.hpp b/include/forcing/GridDataSelector.hpp deleted file mode 100644 index d73ea7bdfa..0000000000 --- a/include/forcing/GridDataSelector.hpp +++ /dev/null @@ -1,259 +0,0 @@ -#pragma once - -#include -#include -#include - -#include -#include - -#include - -struct Cell { - std::uint64_t x = 0; - std::uint64_t y = 0; - std::uint64_t z = 0; - double value = NAN; - - friend inline bool operator<(const Cell& a, const Cell& b) { - // TODO: Not really meaningful - return a.x < b.x || a.y < b.y || a.z < b.z; - } -}; - -using box_t = boost::geometry::model::box; - -struct BoundingBox -{ - BoundingBox(box_t box) - : box_(std::move(box)) - {} - - double xmin() const noexcept - { - return box_.min_corner().get<0>(); - } - - double xmax() const noexcept - { - return box_.max_corner().get<0>(); - } - - double ymin() const noexcept - { - return box_.min_corner().get<1>(); - } - - double ymax() const noexcept - { - return box_.max_corner().get<1>(); - } - - geojson::polygon_t as_polygon() const noexcept { - boost::geometry::box_view view{box_}; - geojson::polygon_t poly; - poly.outer() = { view.begin(), view.end() }; - return poly; - } - - private: - box_t box_; -}; - -struct GridSpecification { - //! Total number of rows (aka y) - std::uint64_t rows; - - //! Total number of columns (aka x) - std::uint64_t columns; - - //! Extent of the grid region (aka min-max corner points) - BoundingBox extent; -}; - -struct SelectorConfig { - //! Initial time for query. - //! @todo Refactor to use std::chrono - time_t init_time; - - //! Duration for query, in seconds. - //! @todo Refactor to use std::chrono - long duration_seconds; - - //! Variable to return from query. - std::string variable_name; - - //! Units for output variable. - std::string variable_units; -}; - -struct GridDataSelector { - - //! Cell-based constructor - GridDataSelector(SelectorConfig config, boost::span cells) - : config_(std::move(config)) - , cells_(cells.begin(), cells.end()) { - } - - /** - * Point-based constructor - * - * Constructs a selector taking only the cells - * from @p grid that correspond to coordinates in @p points. - * - * @param config Selector configuration options - * @param grid Source grid specification - * @param points Target points used for extraction - * - * @todo Implementation - */ - GridDataSelector( - SelectorConfig config, - const GridSpecification& grid, - boost::span points - ) - : config_(std::move(config)) - { - // Using a std::set since points may be close enough that they are within - // the same grid cell. This ensures that each cell is uniquely indexed. - std::set cells; - for (const auto& point : points) { - cells.emplace(Cell{ - /*x=*/position_(point.get<0>(), grid.extent.xmin(), grid.extent.xmax(), grid.columns), - /*y=*/position_(point.get<1>(), grid.extent.ymin(), grid.extent.ymax(), grid.rows), - /*z=*/0UL, - /*value=*/NAN - }); - } - - cells_.assign(cells.begin(), cells.end()); - } - - /** - * Boundary-based constructor - * - * Constructs a selector taking only the cells - * from @p grid that intersect @p polygon. - * - * @param config Selector configuration options - * @param grid Source grid specification - * @param polygon Target polygon used as mask - */ - GridDataSelector( - SelectorConfig config, - const GridSpecification& grid, - const geojson::polygon_t& polygon - ) - : config_(std::move(config)) - { - const auto xmin = grid.extent.xmin(); - const auto xmax = grid.extent.xmax(); - const auto ymin = grid.extent.ymin(); - const auto ymax = grid.extent.ymax(); - const auto ydiff = static_cast(grid.rows) / (ymax - ymin); - const auto xdiff = static_cast(grid.columns) / (xmax - xmin); - - const auto bbox = BoundingBox{ boost::geometry::return_envelope(polygon) }; - for (double row = bbox.ymin(); row < bbox.ymax() - ydiff; row += ydiff) { - for (double col = bbox.xmin(); col < bbox.xmax() - xdiff; row += xdiff) { - const box_t cell_box = { - /*min_corner=*/{ col, row }, - /*max_corner=*/{ col + xdiff, row + ydiff } - }; - - if (boost::geometry::intersects(cell_box, polygon)) { - auto x = position_(col, xmin, xmax, grid.columns); - auto y = position_(row, ymin, ymax, grid.rows); - cells_.emplace_back(Cell{x, y, /*z=*/0UL, /*value=*/NAN}); - } - } - } - } - - /** - * Extent-based constructor - * - * @param config Selector configuration options - * @param grid Source grid specification - * @param extent Target bounding box used as mask - */ - GridDataSelector( - SelectorConfig config, - const GridSpecification& grid, - const BoundingBox& extent - ) - : config_(std::move(config)) - { - const auto col_min = position_(extent.xmin(), grid.extent.xmin(), grid.extent.xmax(), grid.columns); - const auto col_max = position_(extent.xmax(), grid.extent.xmin(), grid.extent.xmax(), grid.columns); - const auto row_min = position_(extent.ymin(), grid.extent.ymin(), grid.extent.ymax(), grid.rows); - const auto row_max = position_(extent.ymax(), grid.extent.ymin(), grid.extent.ymax(), grid.rows); - const auto ncells = (row_max - row_min) * (col_max - col_min); - - cells_.reserve(ncells); - for (auto row = row_min; row < row_max; row++) { - for (auto col = col_min; col < col_max; col++) { - cells_.emplace_back(Cell{/*x=*/col, /*y=*/row, /*z=*/0UL, /*value=*/NAN}); - } - } - } - - GridDataSelector() noexcept = default; - - virtual ~GridDataSelector() = default; - - time_t& initial_time() noexcept { - return config_.init_time; - } - - time_t initial_time() const noexcept { - return config_.init_time; - } - - long& duration() noexcept { - return config_.duration_seconds; - } - - long duration() const noexcept { - return config_.duration_seconds; - } - - std::string& variable() noexcept { - return config_.variable_name; - } - - const std::string& variable() const noexcept { - return config_.variable_name; - } - - std::string& units() noexcept { - return config_.variable_units; - } - - const std::string& units() const noexcept { - return config_.variable_units; - } - - boost::span cells() noexcept { - return cells_; - } - - boost::span cells() const noexcept { - return cells_; - } - - private: - static std::uint64_t position_(double position, double min, double max, std::uint64_t upper_bound) { - if (position < min || position > max) { - return static_cast(-1); - } - - return std::floor((position - min) * (static_cast(upper_bound) / (max - min))); - } - - //! General selector configuration - SelectorConfig config_; - - //! Cells to gather. - std::vector cells_; -}; diff --git a/include/forcing/GriddedDataSelector.hpp b/include/forcing/GriddedDataSelector.hpp new file mode 100644 index 0000000000..fe07893bdc --- /dev/null +++ b/include/forcing/GriddedDataSelector.hpp @@ -0,0 +1,10 @@ +#pragma once + +#include + +struct GriddedDataSelector { + std::string variable_name; + time_t init_time; + long duration; + std::string output_units; +}; diff --git a/include/geojson/JSONGeometry.hpp b/include/geojson/JSONGeometry.hpp index f76fd06bda..9146b57969 100644 --- a/include/geojson/JSONGeometry.hpp +++ b/include/geojson/JSONGeometry.hpp @@ -35,6 +35,8 @@ namespace geojson { typedef bg::model::multi_polygon multipolygon_t; + using box_t = boost::geometry::model::box; + /** * A two dimensional matrix of doubles */ @@ -176,4 +178,4 @@ namespace geojson { } } -#endif // GEOJSON_GEOMETRY_H \ No newline at end of file +#endif // GEOJSON_GEOMETRY_H diff --git a/src/forcing/ForcingsEngineGriddedDataProvider.cpp b/src/forcing/ForcingsEngineGriddedDataProvider.cpp new file mode 100644 index 0000000000..da2af72c63 --- /dev/null +++ b/src/forcing/ForcingsEngineGriddedDataProvider.cpp @@ -0,0 +1,143 @@ +#include + +namespace data_access { + +using Provider = ForcingsEngineGriddedDataProvider; +using BaseProvider = Provider::base_type; + +GridSpecification construct_grid_spec(bmi::Bmi* ptr, int grid) +{ + std::array shape = {-1, -1}; + ptr->GetGridShape(grid, shape.data()); + assert(shape[0] == shape[1]); // since the grid must be uniform rectilinear + + // Allocate a coordinate vector and reuse for X and Y + std::vector coordinate(shape[0]); + + // Get X bounds + ptr->GetGridX(grid, coordinate.data()); + auto xminmax = std::minmax_element(coordinate.begin(), coordinate.end()); + double xmin = *xminmax.first; + double xmax = *xminmax.second; + + // Get Y bounds + ptr->GetGridY(grid, coordinate.data()); + auto yminmax = std::minmax_element(coordinate.begin(), coordinate.end()); + double ymin = *yminmax.first; + double ymax = *yminmax.second; + + return { + /*rows=*/static_cast(shape[0]), + /*columns=*/static_cast(shape[1]), + /*extent=*/{ + /*min_corner=*/{xmin, ymin}, + /*max_corner=*/{xmax, ymax} + } + }; +} + +template +std::uint64_t position_helper( + const geojson::coordinate_t& mask, + const GridSpecification& grid +) +{ + return GridMask::position( + mask.get(), + grid.extent.min_corner().get(), + grid.extent.max_corner().get(), + I == 0 /*x*/ ? grid.columns : grid.rows + ); +} + +GridMask construct_grid_mask(geojson::box_t mask_extent, const GridSpecification& underlying_grid) +{ + auto cmin = position_helper<0>(mask_extent.min_corner(), underlying_grid); + auto cmax = position_helper<0>(mask_extent.max_corner(), underlying_grid); + auto rmin = position_helper<1>(mask_extent.min_corner(), underlying_grid); + auto rmax = position_helper<1>(mask_extent.max_corner(), underlying_grid); + + return { std::move(mask_extent), rmin, cmin, rmax, cmax }; +} + +Provider::ForcingsEngineGriddedDataProvider( + const std::string& init, + std::size_t time_begin_seconds, + std::size_t time_end_seconds, + geojson::box_t mask +) + : BaseProvider(init, time_begin_seconds, time_end_seconds) +{ + // FIXME: assert that var_grid_mask_ is (entirely) within var_grid_ + // (possibly, convert to polygon and use contains predicate) + var_grid_id_ = bmi_->GetVarGrid(get_available_variable_names()[0]); + var_grid_ = construct_grid_spec(bmi_.get(), var_grid_id_); + var_grid_mask_ = construct_grid_mask(std::move(mask), var_grid_); +} + +Provider::data_type Provider::get_value(const selection_type&, data_access::ReSampleMethod) +{ + throw std::runtime_error{"ForcingsEngineGriddedDataProvider::get_value() is not implemented"}; +} + +std::vector Provider::get_values( + const selection_type& selector, + data_access::ReSampleMethod m +) +{ + auto variable = ensure_variable(selector.variable_name); + + if (m != ReSampleMethod::SUM && m != ReSampleMethod::MEAN) { + throw std::runtime_error{ + "ForcingsEngineGriddedDataProvider::get_values(): " + "ReSampleMethod " + std::to_string(m) + " not implemented."}; + } + + const auto duration = std::chrono::seconds{selector.duration}; + + const auto start = clock_type::from_time_t(selector.init_time); + assert(start >= time_begin_); + + const auto end = start + duration; + assert(end <= time_end_); + + std::vector values; + values.reserve(var_grid_mask_.size()); + for (auto current = start; current < end; current += time_step_, bmi_->UpdateUntil((current - start).count())) { + // Get a span over the entire grid + boost::span full = { static_cast(bmi_->GetValuePtr(variable)), var_grid_.rows * var_grid_.columns }; + + // Iterate row by row over the grid, masking the grid columns in each row. + // For each row, we add the grid values to the masked grid values. + for (auto r = var_grid_mask_.rmin; r < var_grid_mask_.rmax; ++r) { + + // Get the starting index of the current row within the full span + // Equation: + ( * ) + const std::size_t row_address = var_grid_mask_.cmin + (r * var_grid_.columns); + + // Get the starting index (pointer address) of the current row within the masked grid + // Equation: + ( * ) + double* const mask_address = values.data() + ((r - var_grid_mask_.rmin) * var_grid_mask_.columns()); + + // Get a span over the current row index on the underlying grid + boost::span row = full.subspan(row_address, var_grid_mask_.columns()); + + // Get a mutable span over the current row index in the masked values + boost::span masked = { mask_address, var_grid_mask_.columns() }; + + // Add grid values to masked values + std::transform(row.begin(), row.end(), masked.begin(), masked.begin(), std::plus{}); + } + } + + if (m == ReSampleMethod::MEAN) { + auto steps = duration / time_step_; + for (auto& v : values) { + v /= steps; + } + } + + return values; +} + +} // namespace data_access