Merge branch 'rasterize' into 'master'

Add rasterize_polygons See merge request isciences/exactextractr!60
isciences · Jun 27, 2022 · 6523ad5 · 6523ad5
2 parents 6d02446 + c5cddcc
commit 6523ad5
Show file tree

Hide file tree

Showing 12 changed files with 269 additions and 5 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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<Rcpp::CharacterVector> p_stat, Rcpp::Nullable<Rcpp::Function> 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
@@ -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
@@ -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<exactextract::bounded_extent> make_grid(const Rcpp::S4 & rast);
+exactextract::Grid<exactextract::bounded_extent> 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
@@ -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 <Rcpp.h>
+
+#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
@@ -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)
+  )
+})