Feature request: allow extract() to do multiple functions

[dfriend21]

I've been wanting to find a way to extract multiple summary stats for polygons that doesn't require calling extract() multiple times, since that's not especially efficient. I've found that I can accomplish this by creating a function that calculates the statistics and then passing that to extract():

library(terra)
v <- vect(system.file("ex/lux.shp", package = "terra"))
r <- rast(system.file("ex/elev.tif", package = "terra"))

f <- function(x){
  return(cbind(mean = mean(x, na.rm = TRUE), 
               sd = sd(x, na.rm = TRUE)))
}
e <- extract(r, v, f)
head(e)
#>      ID elevation elevation.1
#> [1,]  1  467.1052    34.58480
#> [2,]  2  333.8629    68.03058
#> [3,]  3  377.3712    77.14169
#> [4,]  4  373.6000    82.79129
#> [5,]  5  418.6490    48.35488
#> [6,]  6  314.9969    49.02558

Based on some quick investigation into the code, it looks like if you use common functions like mean, median, sum, etc. then it
uses C++ functions that are presumably more optimized (correct me if I'm wrong). So as it is right now there doesn't seem to be any way to calculate two of the optimized C++ functions at the same time, since the function I created will be using the base R functions. It'd be nice if there was a way to do this.

Just a suggestion, feel free to close if you don't want to implement this.

[kadyb]

BTW: Do you know {exactextractr}? This is the fastest package for zonal statistics (but requires {sf} too for vector data).

library("sf")
library("terra")
library("exactextractr")

v <- read_sf(system.file("ex/lux.shp", package = "terra"))
r <- rast(system.file("ex/elev.tif", package = "terra"))

e <- exact_extract(r, v, fun = c("mean", "stdev"))
head(e)
#       mean    stdev
# 1 467.3792 33.98623
# 2 334.6855 68.61759
# 3 377.2069 77.09263
# 4 372.2498 81.72250
# 5 418.7867 48.63316
# 6 314.7698 49.40102

[dfriend21]

Yes, I'm familiar with it. However, I've been struggling to make it work effectively in my situation. I'm working with rasters via HTTP (using vsicurl/) - I want to be able to calculate summary statistics for polygons without having to download the entire raster (some of them are 30+ GB), which is why I really want to calculate the statistics simultaneously rather than with separate extract calls. I tried using exactextractr - in a lot of the test cases I've tried, it works just fine. But I've had a few instances where it hangs. I'm guessing it has something to do with the max_cells_in_memory parameter. The situation where it hangs for me is one where I'm trying to extract data for two polygons that aren't close to each other - so I'm wondering if it's trying to download all of the data covered by the extent, even though the individual polygons aren't that big. That being said, I tried a different test case with two small polygons that were very far away from each other, and it worked just fine. Perhaps if I fiddled around with it a bit more, I could get it to work how I want it to - I haven't spent a ton of time messing with it. But terra::extract() runs quite quickly in the same situation, and since it's been more consistent for me, that's what I'm using right now.

[kadyb]

Ok, I see. I don't know how to fix it, but maybe it is worth creating the issue in {exactextractr} and Dan will suggest a solution.

Another workaround, I think, is to use {data.table} (or {collapse}) and these statistics should be calculated in parallel for groups, which should also be faster than the current {terra}:

library("terra")
library("data.table")

v <- vect(system.file("ex/lux.shp", package = "terra"))
r <- rast(system.file("ex/elev.tif", package = "terra"))

e <- extract(r, v)
e <- setDT(e, key = "ID")
e <- e[, .(mean = mean(elevation, na.rm = TRUE),
           sd = sd(elevation, na.rm = TRUE)),
       by = ID]

[dfriend21]

Thanks for the suggestion. I hadn't thought of just extracting the values and then doing the calculations on that (probably should have...). That's a good alternative.

[rhijmans]

Thank you for the suggestion. You can now do:

library(terra)
#terra 1.7.60
v <- vect(system.file("ex/lux.shp", package = "terra"))
r <- rast(system.file("ex/elev.tif", package = "terra"))
e <- extract(r, v, fun=c("min", "mean", "max"), na.rm=TRUE)
head(e)
#  ID min_elevation mean_elevation max_elevation
#1  1           339       467.1052           547
#2  2           195       333.8629           514
#3  3           256       377.3712           517
#4  4           213       373.6000           520
#5  5           293       418.6490           511
#6  6           164       314.9969           403

[dfriend21]

Great, thanks!