Merge branch 'isciencesgh-12' into 'master'

Resolve isciences/exactextractr#12

See merge request isciences/exactextractr!23

R/NEWS.md

@@ -1,3 +1,7 @@
+# version 0.1.2
+
+- Generate zero-row data frame for polygons not intersecting the raster
+
 # version 0.1.1
 
 - Attempt to revert to static linking against GEOS if dynamic linking fails

R/R/exact_extract.R

@@ -119,6 +119,13 @@ setMethod('exact_extract', signature(x='Raster', y='sf'), function(x, y, fun=NUL
  sum(sapply(a, nchar) == 0)
 }
 
+emptyVector <- function(rast) {
+ switch(substr(raster::dataType(rast), 1, 3),
+ LOG=logical(),
+ INT=integer(),
+ numeric())
+}
+
 .exact_extract <- function(x, y, fun=NULL, ..., include_xy=FALSE, progress=TRUE, max_cells_in_memory=30000000) {
  if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
  warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
@@ -184,24 +191,41 @@ setMethod('exact_extract', signature(x='Raster', y='sf'), function(x, y, fun=NUL
  appfn(sf::st_as_binary(y), function(wkb) {
  ret <- CPP_exact_extract(x, wkb)
 
- vals <- raster::getValuesBlock(x,
- row=ret$row,
- col=ret$col,
- nrow=nrow(ret$weights),
- ncol=ncol(ret$weights))
+ if (length(ret$weights) > 0) {
+ vals <- raster::getValuesBlock(x,
+ row=ret$row,
+ col=ret$col,
+ nrow=nrow(ret$weights),
+ ncol=ncol(ret$weights))
 
- if(is.matrix(vals)) {
- vals <- as.data.frame(vals)
+ if(is.matrix(vals)) {
+ vals <- as.data.frame(vals)
+ } else {
+ vals <- data.frame(value=vals)
+ }
  } else {
- vals <- data.frame(value=vals)
+ # Polygon does not intersect raster.
+ # Construct a zero-row data frame with correct column names/types.
+ vals <- do.call(data.frame, lapply(seq_len(raster::nlayers(x)),
+ function(i) emptyVector(x[[i]])))
+ if (raster::nlayers(x) == 1) {
+ names(vals) <- 'value'
+ } else {
+ names(vals) <- names(x)
+ }
  }
 
  if (include_xy) {
+ if (nrow(vals) == 0) {
+ vals$x <- numeric()
+ vals$y <- numeric()
+ } else {
  x_coords <- raster::xFromCol(x, col=ret$col:(ret$col+ncol(ret$weights) - 1))
  y_coords <- raster::yFromRow(x, row=ret$row:(ret$row+nrow(ret$weights) - 1))
 
  vals$x <- rep.int(x_coords, times=nrow(ret$weights))
  vals$y <- rep(y_coords, each=ncol(ret$weights))
+ }
  }
 
  cov_fracs <- as.vector(t(ret$weights))

R/src/exact_extract.cpp

@@ -136,9 +136,33 @@ Rcpp::List CPP_exact_extract(Rcpp::S4 & rast, const Rcpp::RawVector & wkb) {
  }
  }
 
+ if (nrow > 0) {
+ size_t row_us = 1 + coverage_fractions.grid().row_offset(grid);
+ if (row_us > std::numeric_limits<int>::max()) {
+
+ }
+ }
+
+ int row = NA_INTEGER;
+ int col = NA_INTEGER;
+ if (nrow > 0) {
+ size_t row_us = (1 + coverage_fractions.grid().row_offset(grid));
+ if (row_us > std::numeric_limits<int>::max()) {
+ throw std::runtime_error("Cannot represent row offset as an R integer");
+ }
+ row = static_cast<int>(row_us);
+ }
+ if (ncol > 0) {
+ size_t col_us = (1 + coverage_fractions.grid().col_offset(grid));
+ if (col_us > std::numeric_limits<int>::max()) {
+ throw std::runtime_error("Cannot represent column offset as an R integer");
+ }
+ col = static_cast<int>(col_us);
+ }
+
  return Rcpp::List::create(
- Rcpp::Named("row") = nrow > 0 ? (1 + coverage_fractions.grid().row_offset(grid)) : NA_INTEGER,
- Rcpp::Named("col") = ncol > 0 ? (1 + coverage_fractions.grid().col_offset(grid)) : NA_INTEGER,
+ Rcpp::Named("row") = row,
+ Rcpp::Named("col") = col,
  Rcpp::Named("weights") = weights
  );
 }

R/tests/testthat/test_exact_extract.R

@@ -311,6 +311,13 @@ test_that('We get acceptable default values when processing a polygon that does
  )
  ), crs=sf::st_crs(rast)) # extent of Antarctica in Natural Earth
 
+ # RasterLayer
+ expect_equal(list(data.frame(value=numeric(), coverage_fraction=numeric())),
+ exact_extract(rast, poly))
+
+ expect_equal(list(data.frame(value=numeric(), x=numeric(), y=numeric(), coverage_fraction=numeric())),
+ exact_extract(rast, poly, include_xy=TRUE))
+
  expect_equal(0, exact_extract(rast, poly, function(x, c) sum(x)))
  expect_equal(0, exact_extract(rast, poly, 'count'))
  expect_equal(0, exact_extract(rast, poly, 'sum'))
@@ -321,6 +328,23 @@ test_that('We get acceptable default values when processing a polygon that does
  expect_equal(NA_real_, exact_extract(rast, poly, 'mean'))
  expect_equal(NA_real_, exact_extract(rast, poly, 'min'))
  expect_equal(NA_real_, exact_extract(rast, poly, 'max'))
+
+ # RasterStack
+ rast2 <- as.integer(rast)
+ raster::dataType(rast2) <- 'INT4S'
+
+ stk <- raster::stack(list(q=rast, xi=rast2, area=raster::area(rast)))
+
+ expect_equal(list(data.frame(q=numeric(), xi=integer(), area=numeric(), coverage_fraction=numeric())),
+ exact_extract(stk, poly))
+
+ expect_equal(list(data.frame(q=numeric(), xi=integer(), area=numeric(), x=numeric(), y=numeric(), coverage_fraction=numeric())),
+ exact_extract(stk, poly, include_xy=TRUE))
+
+ exact_extract(stk, poly, function(values, cov) {
+ expect_equal(values, data.frame(q=numeric(), xi=integer(), area=numeric()))
+ expect_equal(cov, numeric())
+ })
 })
 
 test_that('We can optionally get cell center coordinates included in our output', {