Merge pull request #54 from ecohealthalliance/RC1.1.0

Rc1.1.0

.github/workflows/R-CMD-check.yaml

@@ -2,7 +2,7 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main]
+    branches: [main, RC1.1.0]
   pull_request:
     branches: [main]
 

DESCRIPTION

@@ -20,13 +20,12 @@ Description: Provides a drop-in replacement for rasterize() from the 'raster'
 License: MIT + file LICENSE
 URL: https://github.com/ecohealthalliance/fasterize
 BugReports: https://github.com/ecohealthalliance/fasterize/issues
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Suggests: 
     testthat,
     microbenchmark,
     knitr,
     rmarkdown,
-    sf,
     spelling,
     geos
 Depends:

R/fasterize.R

@@ -47,16 +47,15 @@ make_sf <- function(x, attr = NULL) {
 #'   14-16, 1967, Fall Joint Computer Conference. AFIPS '67 (Fall).
 #'   \doi{10.1145/1465611.1465619}
 #' @examples
-#' library(sf)
+#' library(wk)
 #' library(fasterize)
-#' p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
-#' hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
-#' p1 <- list(p1, hole)
-#' p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
-#' p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
-#' pols <- st_sf(value = rep(1,3),
-#'               geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
-#' r <- raster(pols, res = 1)
+#' p123 <- c(paste0("POLYGON ((-180 -20, -140 55, 10 0, -140 -60, -180 -20),", 
+#'                  "(-150 -20, -100 -10, -110 20, -150 -20))"), 
+#'             "POLYGON ((-10 0, 140 60, 160 0, 140 -55, -10 0))", 
+#'             "POLYGON ((-125 0, 0 60, 40 5, 15 -45, -125 0))")
+#' pols <- data.frame(value = seq_along(p123), geometry = wk::as_wkt(p123))
+#' ex <- as.numeric(wk_bbox(pols))[c(1, 3, 2, 4)]
+#' r <- raster::raster(raster::extent(ex), res = 1)
 #' r <- fasterize(pols, r, field = "value", fun="sum")
 #' plot(r)
 #' @export

README.Rmd

@@ -19,7 +19,7 @@ knitr::opts_chunk$set(
 
 # fasterize
 
-Fast sf-to-raster conversion
+Fast polygon-to-raster conversion, burn polygon shapes and/or values into pixels. 
 
 
 
@@ -67,39 +67,47 @@ A `raster()` and `plot()` methods for rasters are re-exported from the [raster p
 ```{r example-1, message=FALSE}
 library(raster)
 library(fasterize)
-library(sf)
-p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
-hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
-p1 <- list(p1, hole)
-p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
-p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
-pols <- st_sf(value = c(1,2,3),
-             geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
-r <- raster(pols, res = 1)
+library(wk)
+library(fasterize)
+p123 <- c(paste0("POLYGON ((-180 -20, -140 55, 10 0, -140 -60, -180 -20),", 
+                  "(-150 -20, -100 -10, -110 20, -150 -20))"), 
+             "POLYGON ((-10 0, 140 60, 160 0, 140 -55, -10 0))", 
+             "POLYGON ((-125 0, 0 60, 40 5, 15 -45, -125 0))")
+pols <- data.frame(value = seq_along(p123), geometry = wk::as_wkt(p123))
+ex <- as.numeric(wk_bbox(pols))[c(1, 3, 2, 4)]
+r <- raster::raster(raster::extent(ex), res = 1)
 r <- fasterize(pols, r, field = "value", fun="sum")
 plot(r)
 ```
 
 ## Performance
 
-Let's compare `fasterize()` to `raster::rasterize()`:
+Let's compare `fasterize()` to  `terra::rasterize()`:
 
 ```{r benchmark, cache=TRUE}
-pols_r <- as(pols, "Spatial")
+pols_t <- terra::vect(p123)
+pols_t$value <- 1:3
+#pols_r <-  as(pols_t, "Spatial")
+tr <- terra::rast(r)
+
 bench <- microbenchmark::microbenchmark(
-  rasterize = r <- raster::rasterize(pols_r, r, field = "value", fun="sum"),
+ # rasterize = r <- raster::rasterize(pols_r, r, field = "value", fun="sum"),
+  terrarize = tr <- terra::rasterize(pols_t, tr, field = "value", fun = "sum"),
   fasterize = f <- fasterize(pols, r, field = "value", fun="sum"),
   unit = "ms"
 )
+
 print(bench, digits = 3)
 ```
 
-It's also quite a bit faster than terra, see the vignette. 
 
 
 How does `fasterize()` do on a large set of polygons? Here I download the IUCN shapefile for the ranges of all terrestrial mammals and generate
 a 1/6 degree world map of mammalian biodiversity by rasterizing all the layers.
 
+
+(this doesn't work anymore because the source data is gone, left as a record 2024-09-25). 
+
 ```{r download, eval=FALSE, cache=TRUE}
 if(!dir.exists("Mammals_Terrestrial")) {
   download.file(
@@ -134,7 +142,8 @@ plot(mammal_raster, axes=FALSE, box=FALSE)
 
 ## About
 
-**fasterize** is developed openly at [EcoHealth Alliance](https://github.com/ecohealthalliance) under the USAID PREDICT project.
+**fasterize** was developed openly at [EcoHealth Alliance](https://github.com/ecohealthalliance) under the USAID PREDICT project.   
+In 
 Please note that this project is released with a [Contributor Code of Conduct](CODE_OF_CONDUCT.md). By participating in this project you agree to abide by its terms.
 
 [![https://www.ecohealthalliance.org/](vignettes/eha-footer.png)](https://www.ecohealthalliance.org/)

fasterize.Rproj

@@ -16,5 +16,6 @@ AutoAppendNewline: Yes
 
 BuildType: Package
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageBuildArgs: --no-build-vignettes
 PackageCheckArgs: --as-cran --no-manual
 PackageRoxygenize: rd,collate,namespace

src/fasterize.cpp

@@ -95,18 +95,6 @@ Rcpp::S4 fasterize_cpp(Rcpp::DataFrame &sf,
     rasterdata.slot("fromdisk") = false;
     rasterdata.slot("haveminmax") = true;
 
-    // new sf only stores ()$input and ()$wkt so we have no basis to grab
-    // a PROJ.4 string from that, just assume they are the same
-    // - this wrongly would *assign* the sf projection to the raster if it
-    // was not NA before MDSumner 2020-03-02
-    // Rcpp::CharacterVector sfproj4 =
-    //   Rcpp::as<Rcpp::StringVector>(
-    //     Rcpp::as<Rcpp::List>(polygons.attr("crs"))["proj4string"]
-    //   );
-    // if(sfproj4[0] != NA_STRING) {
-    //   Rcpp::S4 rcrs(raster1.slot("crs"));
-    //   rcrs.slot("projargs") = sfproj4;
-    // }
 
     return raster1;
 
@@ -146,14 +134,6 @@ Rcpp::S4 fasterize_cpp(Rcpp::DataFrame &sf,
     rasterdata.slot("haveminmax") = true;
     rasterdata.slot("names") = "layer";
 
-    // Rcpp::CharacterVector sfproj4 =
-    //   Rcpp::as<Rcpp::StringVector>(
-    //     Rcpp::as<Rcpp::List>(polygons.attr("crs"))["proj4string"]
-    //   );
-    // if(sfproj4[0] != NA_STRING) {
-    //   Rcpp::S4 rcrs(raster1.slot("crs"));
-    //   rcrs.slot("projargs") = sfproj4;
-    // }
 
     return raster1;
   }

tests/testthat/test-01-inputcheck.R

@@ -1,16 +1,27 @@
 
 context("input checks")
-suppressPackageStartupMessages(library(sf))
+
 library(geos)
-p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
-hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
-p1 <- list(p1, hole)
-p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
-p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
-pols <- st_sf(value = c(1,2,3),
-              geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
-r1 <- raster(pols, res=1)
+# p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
+# hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
+# p1 <- list(p1, hole)
+# p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
+# p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
+# pols <- st_sf(value = c(1,2,3),
+#               geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
 
+##wk::wk_coords(pols) ## ...
+pols_df <- structure(list(feature_id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), 
+                          part_id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), 
+                          ring_id = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L,  3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L), 
+                          x = c(-180, -140, 10, -140,  -180, -150, -100, -110, -150, -10, 140, 160, 140, -10, -125, 0, 40, 15, -125), 
+                          y = c(-20, 55, 0, -60, -20, -20, -10, 20, -20, 0, 60, 0, -55, 0, 0, 60, 5, -45, 0)), row.names = c(NA, 19L), class = "data.frame")
+pols_df$xy <- wk::xy(pols_df$x, pols_df$y)
+pols <- wk::wk_polygon(pols_df, feature_id = pols_df$feature_id, ring_id = pols_df$ring_id, part_id = pols_df$part_id)
+#plot(pols)
+ex <- as.numeric(wk::wk_bbox(pols_df))[c(1, 3, 2, 4)]
+r1 <- raster::raster(raster::extent(ex), res = 1)
+lines <- wk::wk_linestring(pols)
 ## we now accept any wk-handled class
 # test_that("fasterize needs class sf", {
 #   pols_err <- pols
@@ -19,25 +30,23 @@ r1 <- raster(pols, res=1)
 # })
 
 test_that("fasterize likes wkt/wkb/geos", {
-  expect_s4_class(fasterize(wk::as_wkt(pols$geometry), r1), "BasicRaster")
-  expect_s4_class(fasterize(wk::as_wkb(pols$geometry), r1), "BasicRaster")
-  expect_s4_class(fasterize(geos::as_geos_geometry(pols$geometry), r1), "BasicRaster")
-i <- seq_along(pols$geometry)
-  expect_s4_class(fasterize(data.frame(a = i, g = wk::as_wkt(pols$geometry)), r1), "BasicRaster")
-  expect_s4_class(fasterize(data.frame(a = i, wk::as_wkb(pols$geometry)), r1), "BasicRaster")
-  expect_s4_class(fasterize(data.frame(a = i, geos::as_geos_geometry(pols$geometry)), r1), "BasicRaster")
+  expect_s4_class(fasterize(wk::as_wkt(pols), r1), "BasicRaster")
+  expect_s4_class(fasterize(wk::as_wkb(pols), r1), "BasicRaster")
+  expect_s4_class(fasterize(geos::as_geos_geometry(pols), r1), "BasicRaster")
+i <- seq_along(pols)
+  expect_s4_class(fasterize(data.frame(a = i, g = wk::as_wkt(pols)), r1), "BasicRaster")
+  expect_s4_class(fasterize(data.frame(a = i, wk::as_wkb(pols)), r1), "BasicRaster")
+  expect_s4_class(fasterize(data.frame(a = i, geos::as_geos_geometry(pols)), r1), "BasicRaster")
 
 
 })
 test_that("fasterize needs polygons", {
-  lines <- st_sf(value = c(1,2,3),
-                 geometry = st_sfc(lapply(list(p1, p2, p3),
-                                          function(x) st_linestring(x[[1]]))))
+
   expect_error(fasterize(lines, r1),
                "no polygon geometries to fasterize")
 
-    lines_wkb <- data.frame(value = c(1,2,3),
-                 geometry = wk::as_wkb(sf::st_cast(pols$geometry, "MULTILINESTRING")))
+    lines_wkb <- data.frame(value = c(1),
+                 geometry = wk::as_wkb(lines))
 
       expect_error(fasterize(lines_wkb, r1),
                "no polygon geometries to fasterize")
@@ -46,8 +55,8 @@ test_that("fasterize needs polygons", {
 
 })
 
-test_that("field value name is in sf object", {
-  expect_error(fasterize(pols, r1, field="hello"), class="Rcpp::index_out_of_bounds")
+test_that("field value name is in handleable object", {
+  expect_error(fasterize(pols, r1, field="hello"))
 })
 
 test_that("rotated rasters not supported", {
@@ -56,3 +65,71 @@ test_that("rotated rasters not supported", {
   expect_error(fasterize(pols, r1_err),
                "No current support for rotated rasters.")
 })
+
+
+vals <- 1:3
+funs <- c("sum", "first", "last", "min", "max", "count", "any")
+outval <- c(sum(vals), vals[1], vals[3], min(vals), max(vals),
+            length(vals), any(as.logical(vals)))
+
+pols <- data.frame(value = vals, g = pols)
+for (i in seq_along(funs)) {
+
+  test_that(paste(funs[1], "function works"), {
+    rastout <- fasterize(pols, r1, field = "value", fun = funs[i])
+    expect_equal(unname(rastout[60,172]), outval[i])
+  })
+
+}
+
+test_that("disallowed aggregation function is rejected", {
+  invalid_fn_name <- "yo"
+  expect_error(
+    fasterize(pols, r1, field = "value", fun =  invalid_fn_name),
+    paste0("'fun' has an invalid value: ", invalid_fn_name)
+  )
+})
+
+pols$by_1 = c("a", "a", "b")
+pols$by_2 = c(1, 1, 2)
+pols$by_3 = factor(c("a", "a", "b"))
+
+
+
+test_that("'by' argument works", {
+  expect_error(
+    rb <-fasterize(pols, r1, field="value", fun="sum", by ="by_1"), NA)
+  expect_equal(names(rb), unique(pols$by_1))
+  expect_equal(ncol(rb@data@values), length(unique(pols$by_1)))
+})
+
+test_that("'by' layers are equivalent to layers generated separately", {
+  rba <- fasterize(pols, r1, field="value", fun="sum", by ="value")
+  for(i in 1:nrow(pols)) {
+    expect_equal(raster::as.raster(rba[[i]]),
+                 raster::as.raster(fasterize(pols[i,], r1, field="value", fun="sum")))
+  }
+})
+
+test_that("'by' can handle non-character fields", {
+  expect_error(
+    rb_n <- fasterize(pols, r1, field="value", fun="sum", by ="by_2"), NA)
+  expect_error(
+    rb_f <- fasterize(pols, r1, field="value", fun="sum", by ="by_3"), NA)
+  expect_equal(rb_n@data@names, unique(as.character(pols$by_2)))
+  expect_equal(names(rb_f), unique(as.character(pols$by_3)))
+})
+
+test_that("non-NA background values allowed with by", {
+  r <- r1
+  bg <- 20
+  expect_error(
+    f0 <- fasterize(pols, r, field = "value", fun="last", background = bg,
+                    by = "by_1"), NA)
+  expect_equal(unname(f0[[1]][1,1]), bg)
+  expect_equal(f0@data@max, max(bg, max(pols$value)))
+  expect_equal(f0@data@min, min(bg, min(pols$value)))
+})
+
+
+