zonal.stats speed boost!

[pederengelstad]

An issue and a suggestion all wrapped in one!

I've had the great fortune of finding spatialEco very useful for some zonal summaries I've been performing. However, it appeared that the rasterize() function was a bit of a processing speed bottleneck. Recently, a colleague of mine pointed me toward the fasterize library. It's designed to do basically one thing: make rasterizing vectors more efficient. I took the liberty of grabbing the zonal.stats() function code and came up with a slight edit, utilizing fasterize.

zonal.stats <- function(x, y, stat, trace = TRUE, plot = TRUE) {
  if (class(y) != "RasterLayer") 
    stop("y must be a raster object")
  if (class(x) != "SpatialPolygonsDataFrame") 
    stop("x must be a SpatialPolygonsDataFrame object")
  results <- vector()
  for (j in 1:nrow(x)) {
    if (trace == TRUE) {
      cat("processing", j, "of", nrow(x), "\n")
    }
    lsub <- sf::st_as_sf(x[j, ])
    cr <- raster::crop(y, raster::extent(lsub), snap = "out")
    crop.NA <- raster::setValues(cr, NA)
    fr <- fasterize::fasterize(lsub, cr)
    r <- raster::mask(x = cr, mask = fr)
    if (plot == TRUE) {
      plot(r, main = paste("Polygon: ", j, sep = " "))
    }
    r <- raster::values(r)
    r <- stats::na.omit(r)
    if (length(r) < 1) {
      results <- append(results, NA)
    } else {
      results <- append(results, stat(r))
    }
  }
  results
} 

In my testing (400 NPS polygons, 90 m raster) this change dropped the processing time from ~ 60 min to ~ 3 min. That's the only dataset I've tested it on but I'd assume similar results would be evident in other examples.

In any case, I figured I would drop this here to see if you felt it was appropriate to incorporate.

Cheers!

Peder

[jeffreyevans]

I have been meaning to rewrite this function to make it more efficient. The current modifications that I have implemented actually use the velox extract function, which skips the need for rasterization completely and is a massive speed improvement. Unfortunately, the way velox handles its interface to Rcpp makes it difficult to call a function from the velox namespace, thus making it difficult to pass CRAN checks. I am working on the issue and should have a modified version of the function available in the next package release. The fasterize package, until recently, has only been available on GitHub so, I have not incorporated it into any of my functions. However, I will benchmark you solution against my velox implementation and go with the faster solution. Best, Jeff From: pederengelstad <notifications@github.com> Sent: Monday, April 15, 2019 3:57 PM To: jeffreyevans/spatialEco <spatialEco@noreply.github.com> Cc: Subscribed <subscribed@noreply.github.com> Subject: [jeffreyevans/spatialEco] zonal.stats speed boost! (#10) An issue and a suggestion all wrapped in one! I've had the great fortune of finding spatialEco very useful for some zonal summaries I've been performing. However, it appeared that the rasterize() function was a bit of a processing speed bottleneck. Recently, a colleague of mine pointed me toward the fasterize library. It's designed to do basically one thing: make rasterizing vectors more efficient. I took the liberty of grabbing the zonal.stats() function code and came up with a slight edit, utilizing fasterize. zonal.stats <- function(x, y, stat, trace = TRUE, plot = TRUE) { if (class(y) != "RasterLayer") stop("y must be a raster object") if (class(x) != "SpatialPolygonsDataFrame") stop("x must be a SpatialPolygonsDataFrame object") results <- vector() for (j in 1:nrow(x)) { if (trace == TRUE) { cat("processing", j, "of", nrow(x), "\n") } lsub <- sf::st_as_sf(x[j, ]) cr <- raster::crop(y, raster::extent(lsub), snap = "out") crop.NA <- raster::setValues(cr, NA) fr <- fasterize::fasterize(lsub, cr) r <- raster::mask(x = cr, mask = fr) if (plot == TRUE) { plot(r, main = paste("Polygon: ", j, sep = " ")) } r <- raster::values(r) r <- stats::na.omit(r) if (length(r) < 1) { results <- append(results, NA) } else { results <- append(results, stat(r)) } } results } In my testing (400 NPS polygons, 90 m raster) this change dropped the processing time from ~ 60 min to ~ 3 min. That's the only dataset I've tested it on but I'd assume similar results would be evident in other examples. In any case, I figured I would drop this here to see if you felt it was appropriate to incorporate. Cheers! Peder — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#10>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AJaj_2cki5NUWsmeunFNUDwBLvf-gECcks5vhPWngaJpZM4cw9Sb>.