1

I am using the extract function from the raster package in R.

The extract function is taking incredibly long to assign the values to the polygon. In ArcMap, this process is much faster (but of course more time consuming to set up).

CVsteepslope is the raster (100m resolution) and grids (3,000,000 objects) is my polygon.

Here is the code:

CVsteepslope <- CVsteepslope %>% 
  projectRaster(CVsteepslope, crs="EPSG:3035")
# here I am reducing the raster to the polygon extent in an attempt to speed things up
CVsteepslope <- crop(CVsteepslope, extent(gridsdiss), snap="out")
  CVsteepslope <- rasterize(gridsdiss, CVsteepslope)
    CVsteepslope <- mask(x=CVsteepslope, mask=CVp) 
grids$steepslope <- as.numeric(extract(CVsteepslope, grids, method='bilinear', small=TRUE, fun=mean))

Side note: I am using a virtual environment with 128 gb memory, R is only using a third of this.

Tris
  • 77
  • 1
  • 7
  • 1
    You could try the "exactextractr" package from CRAN: "exactextractr: Fast Extraction from Raster Datasets using Polygons". We can't really test anything without your data (or something like it) to see if its faster than raster::extract – Spacedman Apr 14 '21 at 09:28
  • 1
    Did you bother to read the answer to your previous post? You were aimed to exactextractr, by this same person, and responded that that is what you were looking for. I also posted methods and benchmarks that clearly show that you would not want to use raster::extract. https://gis.stackexchange.com/questions/393822/extract-all-bands-as-weighted-mean-values-of-a-multiband-raster-to-polygons-r/393830#393830 – Jeffrey Evans Apr 16 '21 at 14:12
  • Hey @JeffreyEvans I did read the answer to my previous post, apologies for the duplicate. I should have deleted this one, it doesn't specify that I am interested in extracting all of the values from a multiband to a polygon. On the post you linked, I marked yours as the answer. – Tris Apr 19 '21 at 00:42
  • Also, thank you for taking the time to answer my questions :) @JeffreyEvans – Tris Apr 19 '21 at 00:48

2 Answers2

3

There are tons of question on how to speed up extraction in r, for example this. Furthermore you did not give any example date to test you code or to generate better code. However, as @Spacedman said in the comment the best option so far is exact_extract() from exactextractr package (more detail here). Some example below

library(exactextractr)  
library(raster)
library(maptools)
library(sf)
library(terra)
## Example SpatialPolygonsDataFrame
data(wrld_simpl) #polygon of world countries
wrld_simpl_sf <- st_as_sf(wrld_simpl)

Example RasterLayer


r <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, xmx=180, 
ymn=-90, ymx=90)
r[] <- 1:length(r)


#plot, so you can see it
plot(r)

plot(wrld_simpl, add=TRUE)


extract -----------------------------------------------------------------


rr <- rast(r)
v <- vect(wrld_simpl)


system.time(e1 <- raster::extract(r,wrld_simpl,small=T,method="bilinear",fun=mean))#76 
sec
system.time(e2 <- terra::extract(rr,v,method="bilinear",fun=mean))#5 sec
system.time(e3 <- exact_extract(r,wrld_simpl_sf,fun="mean"))#0.61sec

For now, the function takes only polygons object (not points, not lines), but is super fast and accurate (don't need option small=TRUE). It also returns the fraction of cell coverage by polygons

Elia
  • 288
  • 2
  • 7
  • To use terra::extract you need to do something like rr <- rast(r); w <- vect(wrld_simpl); e2 <- extract(rr, w, fun=mean). This is because extract is a generic function. What gets executed depends on the argument(s) supplied. If you do terra::extract(RasterLayer) the raster::extract method will be used (thankfully). – Robert Hijmans Apr 16 '21 at 04:57
  • Sorry, I've forgotten 2 lines of code. I edited my answer, to take into account the class of the input in terra::extract – Elia Apr 16 '21 at 10:06
2

extract function from raster package is quite slow when you are extracting data from a big raster or you have many polygons. A good alternative is using the extract function from package terra. This package is an updated version of raster. It works better with native formats SpatVector and SpatRaster. Therefore, the way to proceed is converting your files to these formats and then extract.

library(terra)

convert to native format


poly<-vect(grids)
rs<-rast(CVsteepslope)


extract


example<-terra::extract(rs,poly)

More info about terra package: https://cran.r-project.org/web/packages/terra/terra.pdf

HMSP
  • 1,630
  • 11
  • 18