From a820b06a8fc28d8de0a0e0a3cf158c6b1ab36281 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Fri, 16 Aug 2024 21:33:13 -0400 Subject: [PATCH] Processor: expose grid_compat_tol --- python/src/exactextract/exact_extract.py | 3 +++ python/src/pybindings/processor_bindings.cpp | 1 + python/tests/test_exact_extract.py | 19 ++++++++++++++++ src/exactextract.cpp | 6 +++++ src/feature_sequential_processor.cpp | 2 +- src/grid.h | 16 +++++++------ src/operation.cpp | 4 ++-- src/operation.h | 2 +- src/processor.h | 6 +++++ src/raster_sequential_processor.cpp | 2 +- test/test_grid.cpp | 24 ++++++++++++++++++++ 11 files changed, 73 insertions(+), 12 deletions(-) diff --git a/python/src/exactextract/exact_extract.py b/python/src/exactextract/exact_extract.py index f4ebce94..ee3a5da0 100644 --- a/python/src/exactextract/exact_extract.py +++ b/python/src/exactextract/exact_extract.py @@ -266,6 +266,7 @@ def exact_extract( include_geom=False, strategy: str = "feature-sequential", max_cells_in_memory: int = 30000000, + grid_compat_tol: float = 1e-3, output: str = "geojson", output_options: Optional[Mapping] = None, progress=False, @@ -304,6 +305,7 @@ def exact_extract( a cost of higher memory usage. max_cells_in_memory: Indicates the maximum number of raster cells that should be loaded into memory at a given time. + grid_compat_tol: require value and weight grids to align within ``grid_compat_tol`` times the smaller of the two grid resolutions output: An :py:class:`OutputWriter` or one of the following strings: - "geojson" (the default): return a list of GeoJSON-like features @@ -350,6 +352,7 @@ def exact_extract( if include_geom: processor.add_geom() processor.set_max_cells_in_memory(max_cells_in_memory) + processor.set_grid_compat_tol(grid_compat_tol) if progress: processor.show_progress(True) diff --git a/python/src/pybindings/processor_bindings.cpp b/python/src/pybindings/processor_bindings.cpp index 7294ce25..bf51b26c 100644 --- a/python/src/pybindings/processor_bindings.cpp +++ b/python/src/pybindings/processor_bindings.cpp @@ -33,6 +33,7 @@ bind_processor(py::module& m) .def("add_col", &Processor::include_col) .def("add_geom", &Processor::include_geometry) .def("process", &Processor::process) + .def("set_grid_compat_tol", &Processor::set_grid_compat_tol) .def("set_max_cells_in_memory", &Processor::set_max_cells_in_memory, py::arg("n")) .def("set_progress_fn", [](Processor& self, py::function fn) { std::function wrapper = [fn](double frac, std::string_view message) { diff --git a/python/tests/test_exact_extract.py b/python/tests/test_exact_extract.py index f1758e28..63a4dbb2 100644 --- a/python/tests/test_exact_extract.py +++ b/python/tests/test_exact_extract.py @@ -1241,3 +1241,22 @@ def test_progress(): squares = [square] * 10 exact_extract(rast, squares, "count", progress=True) + + +def test_grid_compat_tol(): + + values = NumPyRasterSource( + np.arange(9).reshape(3, 3), xmin=0, xmax=1, ymin=0, ymax=1 + ) + weights = NumPyRasterSource( + np.arange(9).reshape(3, 3), xmin=1e-3, xmax=1 + 1e-3, ymin=0, ymax=1 + ) + + square = make_rect(0.1, 0.1, 0.9, 0.9) + + with pytest.raises(RuntimeError, match="Incompatible extents"): + exact_extract(values, square, "weighted_mean", weights=weights) + + exact_extract( + values, square, "weighted_mean", weights=weights, grid_compat_tol=1e-2 + ) diff --git a/src/exactextract.cpp b/src/exactextract.cpp index d282f1fe..81b0ec18 100644 --- a/src/exactextract.cpp +++ b/src/exactextract.cpp @@ -58,6 +58,7 @@ main(int argc, char** argv) bool progress = false; bool nested_output = false; bool include_geom = false; + double grid_compat_tol = std::numeric_limits::quiet_NaN(); app.add_option("-p,--polygons", poly_descriptor, "polygon dataset")->required(true); app.add_option("-r,--raster", raster_descriptors, "raster dataset")->required(true); @@ -72,6 +73,7 @@ main(int argc, char** argv) app.add_flag("--nested-output", nested_output, "nested output"); app.add_option("--include-col", include_cols, "columns from input to include in output"); app.add_flag("--include-geom", include_geom, "include geometry in output"); + app.add_flag("--grid-compat-tol", grid_compat_tol, "grid compatibility tolerance"); app.add_flag("--progress", progress); app.set_config("--config"); @@ -147,6 +149,10 @@ main(int argc, char** argv) proc->include_geometry(); } + if (!std::isnan(grid_compat_tol)) { + proc->set_grid_compat_tol(grid_compat_tol); + } + proc->set_max_cells_in_memory(max_cells_in_memory); proc->show_progress(progress); if (progress) { diff --git a/src/feature_sequential_processor.cpp b/src/feature_sequential_processor.cpp index dfc8b8aa..da53d9e8 100644 --- a/src/feature_sequential_processor.cpp +++ b/src/feature_sequential_processor.cpp @@ -71,7 +71,7 @@ FeatureSequentialProcessor::process() Box feature_bbox = exactextract::geos_get_box(m_geos_context, geom); - auto grid = common_grid(m_operations.begin(), m_operations.end()); + auto grid = common_grid(m_operations.begin(), m_operations.end(), m_grid_compat_tol); if (feature_bbox.intersects(grid.extent())) { // Crop grid to portion overlapping feature diff --git a/src/grid.h b/src/grid.h index 5085606d..a431f320 100644 --- a/src/grid.h +++ b/src/grid.h @@ -19,6 +19,8 @@ #include "box.h" +constexpr double DEFAULT_GRID_COMPAT_TOL = 1e-6; + namespace exactextract { struct infinite_extent { @@ -274,7 +276,7 @@ class Grid } template - Grid common_grid(const Grid& b, double tol = 1e-6) const + Grid common_grid(const Grid& b, double tol = DEFAULT_GRID_COMPAT_TOL) const { if (!compatible_with(b, tol)) { throw std::runtime_error("Incompatible extents."); @@ -303,7 +305,7 @@ class Grid } template - Grid overlapping_grid(const Grid& b, double tol = 1e-6) const + Grid overlapping_grid(const Grid& b, double tol = DEFAULT_GRID_COMPAT_TOL) const { if (!compatible_with(b, tol)) { throw std::runtime_error("Incompatible extents."); @@ -368,19 +370,19 @@ subdivide(const Grid& grid, size_t max_size); template Grid -common_grid(T begin, T end) +common_grid(T begin, T end, double tol = DEFAULT_GRID_COMPAT_TOL) { if (begin == end) { return Grid::make_empty(); } else if (std::next(begin) == end) { - return (*begin)->grid(); + return (*begin)->grid(tol); } return std::accumulate( std::next(begin), end, - (*begin)->grid(), - [](auto& acc, auto& op) { - return acc.common_grid(op->grid()); + (*begin)->grid(tol), + [tol](auto& acc, auto& op) { + return acc.common_grid(op->grid(tol), tol); }); } diff --git a/src/operation.cpp b/src/operation.cpp index fa0e0bc6..162ef29c 100644 --- a/src/operation.cpp +++ b/src/operation.cpp @@ -480,10 +480,10 @@ Operation::intersects(const Box& box) const } Grid -Operation::grid() const +Operation::grid(double compat_tol) const { if (weighted()) { - return values->grid().common_grid(weights->grid()); + return values->grid().common_grid(weights->grid(), compat_tol); } else { return values->grid(); } diff --git a/src/operation.h b/src/operation.h index 4fc25b12..968a9e76 100644 --- a/src/operation.h +++ b/src/operation.h @@ -132,7 +132,7 @@ class Operation /// Returns a newly constructed `Grid` representing the common grid between /// the value and weighting rasters - Grid grid() const; + Grid grid(double ncompat_tol) const; /// Returns the type of the `Operation` result virtual Feature::ValueType result_type() const = 0; diff --git a/src/processor.h b/src/processor.h index e79be0fd..970ecc06 100644 --- a/src/processor.h +++ b/src/processor.h @@ -89,6 +89,11 @@ class Processor m_max_cells_in_memory = n; } + void set_grid_compat_tol(double tol) + { + m_grid_compat_tol = tol; + } + void show_progress(bool val) { m_show_progress = val; @@ -145,6 +150,7 @@ class Processor std::vector m_include_cols; + double m_grid_compat_tol = DEFAULT_GRID_COMPAT_TOL; size_t m_max_cells_in_memory = 1000000L; std::function m_progress_fn; diff --git a/src/raster_sequential_processor.cpp b/src/raster_sequential_processor.cpp index fa73ff73..dfd521b3 100644 --- a/src/raster_sequential_processor.cpp +++ b/src/raster_sequential_processor.cpp @@ -51,7 +51,7 @@ RasterSequentialProcessor::process() read_features(); populate_index(); - auto grid = common_grid(m_operations.begin(), m_operations.end()); + auto grid = common_grid(m_operations.begin(), m_operations.end(), m_grid_compat_tol); auto subgrids = subdivide(grid, m_max_cells_in_memory); for (std::size_t i = 0; i < subgrids.size(); i++) { diff --git a/test/test_grid.cpp b/test/test_grid.cpp index d4a42992..41e6a6bd 100644 --- a/test/test_grid.cpp +++ b/test/test_grid.cpp @@ -252,6 +252,30 @@ TEST_CASE("Grid compatibility with reduced tolerance") CHECK(a.compatible_with(b, tol)); CHECK(b.compatible_with(a, tol)); + + CHECK_NOTHROW(a.common_grid(b, tol)); + + struct DummyOp + { + public: + DummyOp(const Grid& grid) + : g(grid) + { + } + + const Grid& grid(double) const + { + return g; + } + + private: + const Grid& g; + }; + + std::vector> grids; + grids.push_back(std::make_unique(a)); + grids.push_back(std::make_unique(b)); + CHECK_NOTHROW(common_grid(grids.begin(), grids.end(), tol)); } TEST_CASE("Grid compatibility (empty grid)", "[grid]")