From c5cddcccb0644fb7571d28c8dd45bb5c116e2203 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Tue, 14 Jun 2022 21:44:17 -0400 Subject: [PATCH] Add rasterize_polygons --- NAMESPACE | 1 + NEWS.md | 2 + R/RcppExports.R | 4 ++ R/rasterize.R | 83 ++++++++++++++++++++++++++++ src/Makevars.in | 2 +- src/Makevars.ucrt | 2 +- src/Makevars.win | 2 +- src/RcppExports.cpp | 17 ++++++ src/raster_utils.cpp | 8 ++- src/raster_utils.h | 3 +- src/rasterize.cpp | 54 +++++++++++++++++++ tests/testthat/test_rasterize.R | 96 +++++++++++++++++++++++++++++++++ 12 files changed, 269 insertions(+), 5 deletions(-) create mode 100644 R/rasterize.R create mode 100644 src/rasterize.cpp create mode 100644 tests/testthat/test_rasterize.R diff --git a/NAMESPACE b/NAMESPACE index 74b5b9c..e71eb01 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ exportMethods(coverage_fraction) exportMethods(exact_extract) exportMethods(exact_resample) +exportMethods(rasterize_polygons) import(raster) import(sf) importFrom(Rcpp,evalCpp) diff --git a/NEWS.md b/NEWS.md index 09cc570..e1614fa 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,8 @@ - Add `frac` and `weighted_frac` summary operations, for summarizing categorical rasters - Add `colname_fun` argument to `exact_extract`, for controlling constructed column names +- Add `rasterize_polygons` function to rasterize a polygon coverage, assigning cells to + polygon with greatest coverage area # version 0.8.2 diff --git a/R/RcppExports.R b/R/RcppExports.R index e504ccd..f99952b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,6 +13,10 @@ CPP_stats <- function(rast, weights, wkb, default_value, default_weight, coverag .Call('_exactextractr_CPP_stats', PACKAGE = 'exactextractr', rast, weights, wkb, default_value, default_weight, coverage_areas, p_area_method, stats, max_cells_in_memory, grid_compat_tol, quantiles) } +CPP_update_max_coverage <- function(extent, res, max_coverage, max_coverage_index, tot_coverage, wkb, index) { + invisible(.Call('_exactextractr_CPP_update_max_coverage', PACKAGE = 'exactextractr', extent, res, max_coverage, max_coverage_index, tot_coverage, wkb, index)) +} + CPP_resample <- function(rast_in, rast_out, p_stat, p_fun, coverage_area, area_method) { .Call('_exactextractr_CPP_resample', PACKAGE = 'exactextractr', rast_in, rast_out, p_stat, p_fun, coverage_area, area_method) } diff --git a/R/rasterize.R b/R/rasterize.R new file mode 100644 index 0000000..8d92a8a --- /dev/null +++ b/R/rasterize.R @@ -0,0 +1,83 @@ +# Copyright (c) 2022 ISciences, LLC. +# All rights reserved. +# +# This software is licensed under the Apache License, Version 2.0 (the "License"). +# You may not use this file except in compliance with the License. You may +# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +setGeneric("rasterize_polygons", function(x, y, ...) + standardGeneric("rasterize_polygons")) + +#' Create a raster approximation of a polygon coverage +#' +#' Returns a raster whose values indicate the index of the polygon covering +#' each cell. Where multiple polygons cover the same cell, the index of the +#' polygon covering the greatest area will be used, with the lowest index +#' returned in the case of ties. Cells that are not covered by any polygon, +#' or whose total covered fraction is less than \code{min_coverage}, will be +#' set to \code{NA}. +#' +#' @param x a \code{sf} or \code{sfc} object with polygonal geometries +#' +#' @param y a (possibly empty) \code{RasterLayer} whose resolution and +#' extent will be used for the generated \code{RasterLayer}. +#' @param min_coverage minimum fraction of a cell that must be covered by +#' polygons to be included in the output +#' @return a \code{RasterLayer} or \code{SpatRaster}, consistent with the type of \code{y} +#' @name rasterize_polygons +NULL + +.rasterize_polygons <- function(x, y, min_coverage = 0) { + num_rows <- terra::nrow(y) + num_cols <- terra::ncol(y) + + if (min_coverage == 1.0) { + # account for error in sum of single-precision plots + min_coverage <- (min_coverage - 1e-6) + } + + max_coverage_values <- matrix(0.0, nrow = num_rows, ncol = num_cols) + max_coverage_indexes <- matrix(NA_integer_, nrow = num_rows, ncol = num_cols) + tot_coverage <- matrix(0.0, nrow = num_rows, ncol = num_cols) + + dest_extent <- .extent(y) + dest_res <- .res(y) + + geoms <- sf::st_geometry(x) + + for (i in seq_along(geoms)) { + wkb <- sf::st_as_binary(geoms[[i]], EWKB=TRUE) + + CPP_update_max_coverage(dest_extent, dest_res, max_coverage_values, max_coverage_indexes, tot_coverage, wkb, i) + } + + max_coverage_indexes[tot_coverage < min_coverage] <- NA + + if (inherits(y, 'SpatRaster')) { + ret <- terra::rast(y) + terra::values(ret) <- max_coverage_indexes + } else { + ret <- raster::raster(y) + raster::values(ret) <- max_coverage_indexes + } + + return(ret) +} + +#' @import sf +#' @import raster +#' @useDynLib exactextractr +#' @rdname rasterize_polygons +#' @export +setMethod('rasterize_polygons', signature(x='sf', y='RasterLayer'), .rasterize_polygons) + +#' @useDynLib exactextractr +#' @rdname rasterize_polygons +#' @export +setMethod('rasterize_polygons', signature(x='sf', y='SpatRaster'), .rasterize_polygons) diff --git a/src/Makevars.in b/src/Makevars.in index 8a23051..bc737bd 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -6,4 +6,4 @@ PKG_LIBS = @PKG_LIBS@ EE = exactextract/src SOURCES = $(EE)/area.cpp $(EE)/box.cpp $(EE)/cell.cpp $(EE)/coordinate.cpp $(EE)/floodfill.cpp $(EE)/geos_utils.cpp $(EE)/grid.cpp $(EE)/perimeter_distance.cpp $(EE)/raster_cell_intersection.cpp $(EE)/side.cpp $(EE)/traversal.cpp $(EE)/traversal_areas.cpp $(EE)/weighted_quantiles.cpp -OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o resample.o $(SOURCES:.cpp=.o) +OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o rasterize.o resample.o $(SOURCES:.cpp=.o) diff --git a/src/Makevars.ucrt b/src/Makevars.ucrt index bdf4761..33d5f40 100644 --- a/src/Makevars.ucrt +++ b/src/Makevars.ucrt @@ -6,7 +6,7 @@ PKG_CXXFLAGS = -std=c++14 EE = exactextract/src SOURCES = $(EE)/area.cpp $(EE)/box.cpp $(EE)/cell.cpp $(EE)/coordinate.cpp $(EE)/floodfill.cpp $(EE)/geos_utils.cpp $(EE)/grid.cpp $(EE)/perimeter_distance.cpp $(EE)/raster_cell_intersection.cpp $(EE)/side.cpp $(EE)/traversal.cpp $(EE)/traversal_areas.cpp $(EE)/weighted_quantiles.cpp -OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o resample.o $(SOURCES:.cpp=.o) +OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o rasterize.o resample.o $(SOURCES:.cpp=.o) all: clean diff --git a/src/Makevars.win b/src/Makevars.win index d2b8297..166afb1 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -9,7 +9,7 @@ PKG_CXXFLAGS = -std=c++14 -I$(RWINLIB)/include/geos EE = exactextract/src SOURCES = $(EE)/area.cpp $(EE)/box.cpp $(EE)/cell.cpp $(EE)/coordinate.cpp $(EE)/floodfill.cpp $(EE)/geos_utils.cpp $(EE)/grid.cpp $(EE)/perimeter_distance.cpp $(EE)/raster_cell_intersection.cpp $(EE)/side.cpp $(EE)/traversal.cpp $(EE)/traversal_areas.cpp $(EE)/weighted_quantiles.cpp -OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o resample.o $(SOURCES:.cpp=.o) +OBJECTS = RcppExports.o exact_extract.o raster_utils.o coverage_fraction.o rasterize.o resample.o $(SOURCES:.cpp=.o) all: clean winlibs diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index fcac7a6..e6c691a 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -71,6 +71,22 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// CPP_update_max_coverage +void CPP_update_max_coverage(Rcpp::NumericVector& extent, Rcpp::NumericVector& res, Rcpp::NumericMatrix& max_coverage, Rcpp::IntegerMatrix& max_coverage_index, Rcpp::NumericMatrix& tot_coverage, const Rcpp::RawVector& wkb, int index); +RcppExport SEXP _exactextractr_CPP_update_max_coverage(SEXP extentSEXP, SEXP resSEXP, SEXP max_coverageSEXP, SEXP max_coverage_indexSEXP, SEXP tot_coverageSEXP, SEXP wkbSEXP, SEXP indexSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector& >::type extent(extentSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericVector& >::type res(resSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix& >::type max_coverage(max_coverageSEXP); + Rcpp::traits::input_parameter< Rcpp::IntegerMatrix& >::type max_coverage_index(max_coverage_indexSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix& >::type tot_coverage(tot_coverageSEXP); + Rcpp::traits::input_parameter< const Rcpp::RawVector& >::type wkb(wkbSEXP); + Rcpp::traits::input_parameter< int >::type index(indexSEXP); + CPP_update_max_coverage(extent, res, max_coverage, max_coverage_index, tot_coverage, wkb, index); + return R_NilValue; +END_RCPP +} // CPP_resample Rcpp::S4 CPP_resample(Rcpp::S4& rast_in, Rcpp::S4& rast_out, Rcpp::Nullable p_stat, Rcpp::Nullable p_fun, bool coverage_area, std::string area_method); RcppExport SEXP _exactextractr_CPP_resample(SEXP rast_inSEXP, SEXP rast_outSEXP, SEXP p_statSEXP, SEXP p_funSEXP, SEXP coverage_areaSEXP, SEXP area_methodSEXP) { @@ -92,6 +108,7 @@ static const R_CallMethodDef CallEntries[] = { {"_exactextractr_CPP_coverage_fraction", (DL_FUNC) &_exactextractr_CPP_coverage_fraction, 3}, {"_exactextractr_CPP_exact_extract", (DL_FUNC) &_exactextractr_CPP_exact_extract, 17}, {"_exactextractr_CPP_stats", (DL_FUNC) &_exactextractr_CPP_stats, 11}, + {"_exactextractr_CPP_update_max_coverage", (DL_FUNC) &_exactextractr_CPP_update_max_coverage, 7}, {"_exactextractr_CPP_resample", (DL_FUNC) &_exactextractr_CPP_resample, 6}, {NULL, NULL, 0} }; diff --git a/src/raster_utils.cpp b/src/raster_utils.cpp index ba8372a..2631b87 100644 --- a/src/raster_utils.cpp +++ b/src/raster_utils.cpp @@ -1,4 +1,4 @@ -// Copyright (c) 2018-2020 ISciences, LLC. +// Copyright (c) 2018-2022 ISciences, LLC. // All rights reserved. // // This software is licensed under the Apache License, Version 2.0 (the "License"). @@ -23,6 +23,12 @@ Grid make_grid(const Rcpp::S4 & rast) { Rcpp::NumericVector extent = extentFn(rast); Rcpp::NumericVector res = resFn(rast); + return make_grid(extent, res); +} + +Grid make_grid(const Rcpp::NumericVector & extent, + const Rcpp::NumericVector & res) +{ return {{ extent[0], extent[1], diff --git a/src/raster_utils.h b/src/raster_utils.h index 7fd5e47..7aeb143 100644 --- a/src/raster_utils.h +++ b/src/raster_utils.h @@ -1,4 +1,4 @@ -// Copyright (c) 2018-2021 ISciences, LLC. +// Copyright (c) 2018-2022 ISciences, LLC. // All rights reserved. // // This software is licensed under the Apache License, Version 2.0 (the "License"). @@ -23,6 +23,7 @@ // Construct a grid corresponding to 'rast' exactextract::Grid make_grid(const Rcpp::S4 & rast); +exactextract::Grid make_grid(const Rcpp::NumericVector & extent, const Rcpp::NumericVector & res); // Return the number of layers in 'rast' int get_nlayers(Rcpp::S4 & rast); diff --git a/src/rasterize.cpp b/src/rasterize.cpp new file mode 100644 index 0000000..c6a4661 --- /dev/null +++ b/src/rasterize.cpp @@ -0,0 +1,54 @@ +// Copyright (c) 2022 ISciences, LLC. +// All rights reserved. +// +// This software is licensed under the Apache License, Version 2.0 (the "License"). +// You may not use this file except in compliance with the License. You may +// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// [[Rcpp::plugins("cpp14")]] +#include + +#include "geos_r.h" +#include "raster_utils.h" + +#include "exactextract/src/raster_cell_intersection.h" + +// [[Rcpp::export]] +void CPP_update_max_coverage(Rcpp::NumericVector & extent, + Rcpp::NumericVector & res, + Rcpp::NumericMatrix & max_coverage, + Rcpp::IntegerMatrix & max_coverage_index, + Rcpp::NumericMatrix & tot_coverage, + const Rcpp::RawVector & wkb, + int index) +{ + GEOSAutoHandle geos; + + auto grid = make_grid(extent, res); + + auto coverage_fraction = exactextract::raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get()); + + auto ix = grid.row_offset(coverage_fraction.grid()); + auto jx = grid.col_offset(coverage_fraction.grid()); + + for (size_t i = 0; i < coverage_fraction.rows(); i++) { + for (size_t j = 0; j < coverage_fraction.cols(); j++) { + auto cov = coverage_fraction(i, j); + if (cov > 0) { + tot_coverage(i + ix, j + jx) += cov; + if (cov > max_coverage(i + ix, j + jx)) { + max_coverage(i + ix, j + jx) = cov; + max_coverage_index(i + ix, j + jx) = index; + } + } + } + } +} + + diff --git a/tests/testthat/test_rasterize.R b/tests/testthat/test_rasterize.R new file mode 100644 index 0000000..197f9c9 --- /dev/null +++ b/tests/testthat/test_rasterize.R @@ -0,0 +1,96 @@ +# Copyright (c) 2022 ISciences, LLC. +# All rights reserved. +# +# This software is licensed under the Apache License, Version 2.0 (the "License"). +# You may not use this file except in compliance with the License. You may +# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +context('rasterize_polygons') + +test_that('value is assigned to polygon with greatest coverage area', { + polys <- + st_as_sf( + data.frame(id = 1:3, + geom = c( + 'POLYGON ((10 0, 10 5, 5 5, 10 0))', + 'POLYGON ((0 0, 10 0, 5 5, 1 10, 0 10, 0 0))', + 'POLYGON ((5 5, 10 5, 10 10, 1 10, 5 5))' + )), + wkt = 'geom' + ) + + rt <- terra::rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, res = 2) + + r <- rasterize_polygons(polys, rt) + + # lower-right is a tie, so it goes to the first feature encountered + expect_equal( + extract(r, cbind(9, 1))[[1]], + 1) + + # center tell is touched by all three, goes to polygon that covers the greatest area + expect_equal( + extract(r, cbind(5, 5))[[1]], + 2) +}) + +test_that('min_coverage excludes cells with small coverage area', { + rt <- terra::rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, res = 1) + circ <- make_circle(5, 5, 3.5, crs=st_crs(rt)) + circ_pieces <- st_sf(st_intersection(circ, st_make_grid(circ, 1))) + + cfrac <- coverage_fraction(rt, circ, crop = FALSE)[[1]] + + # by default, all touched cells are included in output + r <- rasterize_polygons(circ_pieces, rt) + expect_equal( + values(cfrac) > 0, + !is.na(values(r)) + ) + + # min_coverage excludes cells with small coverage area + r <- rasterize_polygons(circ_pieces, rt, min_coverage = 0.5) + expect_equal( + values(cfrac) > 0.5, + !is.na(values(r)) + ) +}) + +test_that('input type is preserved', { + rt <- terra::rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, res = 2) + rr <- raster::raster(rt) + + circ <- st_sf(make_circle(5, 5, 3.5, crs=st_crs(rt))) + + r <- rasterize_polygons(circ, rt) + expect_s4_class(r, 'SpatRaster') + + r <- rasterize_polygons(circ, rr) + expect_s4_class(r, 'RasterLayer') +}) + +test_that('no error when polygon does not intersect raster', { + rt <- terra::rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, res = 2, crs=NA) + + circ <- st_sf(make_circle(500, 500, 3.5, crs=st_crs(rt))) + + r <- rasterize_polygons(circ, rt) + + expect_true(all(is.na(values(r)))) +}) + +test_that('no error when polygon partially intersects raster', { + rt <- terra::rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, res = 2, crs=NA) + + circ <- st_sf(make_circle(10, 5, 3.5, crs=st_crs(rt))) + + expect_invisible( + r <- rasterize_polygons(circ, rt) + ) +})