Add references to exactextractr

[Lvulis]

The exactextractr package enables super fast zonal extraction of raster values from polygons. It also reports the coverage fraction of each cell by the input polygon which some could find useful. It works with sf polygons and terra/raster raster data. From the benchmarking reported on the package GitHub, it seems that the improvement with exactextractr::exact_extract over terra::extract is not as significant as the improvement over raster::extract, but could still be useful if someone is working with very large datasets.

Maybe a reference to the package could be added in Chapter 6.

[Robinlovelace]

👍 to that, it's well-maintained, fast and clearly useful. You up for taking a look at this @Lvulis? It's more on Jakub's area, what do you reckon? Thanks for the idea!

[kadyb]

It is also worth mentioning that {exactextractr} is more precise than {raster}, {terra} or {stars}. See: 1, 2

[Nowosad]

@kadyb including when terra::extract() has exact = TRUE?

[kadyb]

There are some differences in results, but are they significant?

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

f = system.file("ex/elev.tif", package = "terra")
r = rast(f)

geom = spatSample(r, size = 3, as.points = TRUE, na.rm = TRUE)
geom = buffer(geom, width = 2000)

extract(r, geom, fun = "mean", exact = FALSE)
#>   ID elevation
#> 1  1  379.6957
#> 2  2  451.9130
#> 3  3  457.9565

extract(r, geom, fun = "mean", exact = TRUE)
#>      ID elevation
#> [1,]  1  380.0747
#> [2,]  2  452.3375
#> [3,]  3  460.0331

exact_extract(r, st_as_sf(geom), fun = "mean")
#> [1] 380.0745 452.3372 460.0329

[Lvulis]

Those are O(10^-6) differences and should be negligible.
I can take a run at writing something integrated with the text and submit it here? Don't have the book source code downloaded currently.

[kadyb]

The timings for precise aggregation are also interesting. It seems that the {terra} algorithm is much less efficient.

library("sf")
library("terra")
library("exactextractr") # devel version from GitHub

n = 1000
# in-memory raster
r = rast(nrows = n, ncols = n, vals = rnorm(n^2),
         xmin = 0, xmax = n,  ymin = 0, ymax = n)

s = 1 # numer of points
geom = matrix(c(sample(n, size = s), sample(n, size = s)), ncol = 2)
geom = st_as_sf(as.data.frame(geom), coords = c(1, 2))
geom = st_buffer(geom, dist = 300)
geom2 = vect(geom)

system.time(x <- exact_extract(r, geom, fun = "mean", progress = FALSE))
#> x == 0.00179644
#>  user  system elapsed
#> 0.041   0.008   0.048

system.time(y1 <- extract(r, geom2, fun = "mean", exact = FALSE))
#> y1 == 0.001754078
#>  user  system elapsed
#> 0.971   0.000   0.972

system.time(y2 <- extract(r, geom2, fun = "mean", exact = TRUE))
#> y2 == 0.00179644
#>   user  system elapsed
#> 21.520   0.175  21.697

[Nowosad]

@kadyb and what are the timings without aggregation?

[kadyb]

The aggregate statistic has a marginal impact (at least in this case):

## `geom` and `geom2` consist of 1 geometry

system.time(x <- exact_extract(r, geom, progress = FALSE))
#>  user  system elapsed
#> 0.042   0.001   0.042

system.time(y1 <- extract(r, geom2, exact = FALSE))
#>  user  system elapsed
#> 0.614   0.001   0.615

system.time(y2 <- extract(r, geom2, exact = TRUE))
#>   user  system elapsed
#> 20.979   0.451  21.430

[Lvulis]

@kadyb Your comment motivated me to double check the benchmarking example that's reported on the package GitHub as the discrepancy didn't seem so large. I just noticed that in the example, the close to similar performance between {terra} and {exactextractr} is with terra::extract(exact=F). Running a subset of the example with exact=T yields a ~300x runtime improvement, similar to the 500x improvement you reported in that synthetic data. It seems the package does give a very significant performance boost, so well worth including.

library(raster)
library(sf)
library(exactextractr)
library(microbenchmark)

brazil <- st_as_sf(getData('GADM', country='BRA', level=1))
brazil_spat <- as(brazil, 'SpatVector')

prec_rast <- getData('worldclim', var='prec', res=10)
prec_terra <- terra::rast(prec_rast) 
microbenchmark(
  terra::extract(prec_terra, brazil_spat, mean, na.rm = TRUE, exact = T),
  exact_extract(prec_terra, brazil, 'mean', progress = FALSE),
  times = 10)

Unit: milliseconds
                                                                   expr
 terra::extract(prec_terra, brazil_spat, mean, na.rm = TRUE, exact = T)
            exact_extract(prec_terra, brazil, "mean", progress = FALSE)
       min         lq       mean     median         uq         max neval
 75268.602 75995.8173 90945.6605 78579.3624 81447.6076 203122.7589    10
   212.249   234.2676   263.9664   245.0914   259.0936    455.1932    10

[kadyb]

Did you test the stable version from CRAN or development from GitHub? Recently there have been some fixes (isciences/exactextractr@5b88ecf) that make {exactextractr} even faster for SpatRaster object.

terra::extract(exact = FALSE) processing time is just similar to exactextractr::exact_extract() but the biggest differences is for terra::extract(exact = TRUE). Here is one more benchmark from me on satellite data: https://github.com/kadyb/raster-benchmark but extract(exact = FALSE) only.

[Lvulis]

The previous benchmark was using the stable version from CRAN. With the development version, for the same example above I have almost the same result. it seems to be a very minor improvement, there was one iteration that was very slow in the previous example and dragged up the mean, but otherwise similar speed.

> microbenchmark(
+   terra::extract(prec_terra, brazil_spat, mean, na.rm = TRUE, exact = T),
+   exact_extract(prec_terra, brazil, 'mean', progress = FALSE),
+   times = 10)
Unit: milliseconds
                                                                   expr
 terra::extract(prec_terra, brazil_spat, mean, na.rm = TRUE, exact = T)
            exact_extract(prec_terra, brazil, "mean", progress = FALSE)
        min         lq      mean     median         uq
 77587.2348 77971.3209 78466.483 78320.0923 79112.2975
   207.5322   231.2555   246.865   243.4413   249.6305
        max neval
 79635.6012    10
   340.1172    10

[kadyb]

You can check it on larger data set, because these differences may not be noticeable on small dataset. For this example I got speedup from 17.857 s to 5.572 s on my PC.

[Lvulis]

EDIT - I compared the two stable release of {exactextractr} (v0.8.2) versus the GH version (0.9) and did not really find improvement. Also fixed variable names below. See old comment that I did compare against {terra} and found {exactextractr} to be faster.

OLD:
Good idea. I did not do the comparison with different {exactextractr} versions but interestingly I found {exactextractr} to be significantly faster than {terra} with this large dataset, even with exact=F in the extract call. I am a little hesitant to run the extract=T and do replicate benchmarks for how long it'll take.

download.file("https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_2.5m_prec.zip",
              "prec_2.5m.zip")
unzip("prec_2.5m.zip", exdir = "prec_2.5m")
f = list.files("prec_2.5m", pattern = "\\.tif$", full.names = TRUE)

prec = rast(f)
prec = prec * 1

brazil <- st_as_sf(getData('GADM', country='BRA', level=1))
brazil_spat <- as(brazil, 'SpatVector')

microbenchmark(
  terra::extract(prec, brazil_spat, mean, na.rm = TRUE, exact = F),
  exact_extract(prec, brazil, 'mean', progress = FALSE),
  times = 10)

Unit: seconds
                                                         expr
 terra::extract(prec, brazil_spat, mean, na.rm = TRUE, exact = F)
        exact_extract(prec, brazil, "mean", progress = FALSE)
      min       lq     mean   median       uq       max
 9.044934 9.168709 9.561833 9.446062 9.607415 11.124775
 2.037946 2.045342 2.148452 2.115836 2.224045  2.365771
 neval
    10
    10
# exact_extract from GH development version 0.9.0


> microbenchmark(
+   terra::extract(prec, brazil_spat, mean, na.rm = TRUE, exact = F),
+   exact_extract(prec, brazil, 'mean', progress = FALSE),
+   times = 10)
Error in x$.self$finalize() : attempt to apply non-function
Unit: seconds
                                                             expr
 terra::extract(prec, brazil_spat, mean, na.rm = TRUE, exact = F)
            exact_extract(prec, brazil, "mean", progress = FALSE)
      min       lq     mean   median       uq      max
 9.288340 9.578843 9.615942 9.630635 9.733694 9.810699
 2.072167 2.090151 2.375785 2.168442 2.418485 3.256852
 neval
    10
    10
# exact_extract from stable version 0.8.2