2

I am trying to extract the values of a multiband raster to a polygon in R.

ACVheatingdegreedayS <- stack("location/multibandraster.tif")
ACVheatingdegreedayS <- ACVheatingdegreedayS %>% 
  projectRaster(ACVheatingdegreedayS, crs="EPSG:3035")

( zCVhdd1 <- CVheatingdegreedayS[[1]] )

I want to avoid that last line of code where I break the multiband into a single band. There are 70 bands, so creating separate rasters would be very time consuming.

Anyone know of a way to extract all bands of a multiband raster to a polygon?

7
  • 1
    II would recommenced the exactextractr::exact_extract function. It benchmarks, by far, the fastest of all the options. There is also raster::extract, terra::extract and a pure tidy approach using tabularaster::cellnumbers but, it is not as straight forward to use as the others. Apr 14, 2021 at 14:27
  • Where's the polygon here? Do we care about the projectRaster? Is CVheatingdegreedayS a typo for ACVheatingdegreedayS? Why can't you loop over 1 to 70? Have you not done a basic R course where you find out about loops?
    – Spacedman
    Apr 14, 2021 at 14:38
  • thanks! the extact_extractr package was what I was looking for, not so much advice on how to make a loop.
    – Tris
    Apr 15, 2021 at 6:56
  • Question still isn't clear. Where's the weights, and where's the values you are finding the weighted mean of?
    – Spacedman
    Apr 15, 2021 at 8:37
  • I would like the average of all pixels values covering each polygon weighted by how much they cover the pixel. I am interested in creating a variable in the polygons for each raster band @Spacedman
    – Tris
    Apr 15, 2021 at 14:29

1 Answer 1

2

Here are some general methods for multi-band raster extraction and summary along with benchmarks. The main difference in multi-band verses single-band is how the data is summarized. If you want a mean weighted by the intersection of cell sizes the easiest route is exactextractr::exact_extract. This function actually has a weighted.mean option, as the output, that is calculated in the C++ routine (detailed in functions help). The "raw" results from this function includes the fractional intersection for each intersecting cell so, you could calculate something yourself.

I would note that the different methods of summary are not really adding to the benchmark results, times are reflective of how long it takes to extract the data from the raster stack.

Add libraries and create some data

library(sf)
library(raster)
library(terra)
library(exactextractr)
library(tabularaster)
library(dplyr)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
  nc <- st_cast(nc, "POLYGON")

i=500; j=500
r <- do.call(raster::stack, replicate(20, 
             raster::raster(matrix(runif(i*j), i, j)))) 
    extent(r) <- extent(nc)
      proj4string(r) <- st_crs(nc)$proj4string

plot(r[[1]])
  plot(st_geometry(nc), add=TRUE)

We can now run through and benchmark each method from raster, exactextract, terra and tabularaster. For a gut check, the results of each method are held in the eresult list object.

eresults <- list()

# raster::extract 
system.time({
  v.raster <- raster::extract(r, nc)
  eresults[["raster"]] <-   
    do.call("rbind", lapply(v.raster, function(x) apply(x, MARGIN=2, mean))) 
})

# terra::extract 
system.time({
  v.terra <- terra::extract( rast(r), vect(nc))
  eresults[["terra"]] <-    
  do.call("rbind", lapply(unique(v.terra[,"ID"]), function(i) 
        apply(v.terra[v.terra$ID==i,], MARGIN=2, mean))) 
})

# exactextractr 
system.time({
  v.exact <- exactextractr::exact_extract(r, nc) 
  eresults[["exactextract"]] <-    
    do.call("rbind", lapply(v.exact, function(x) apply(x, MARGIN=2, mean)))[,-21] 
})

# tabularaster using raster 
system.time({
  index <- cellnumbers(r[[1]], nc)
    index %>% 
      group_by(object_) %>% 
        count()     
  result <- index %>% 
     mutate(pixelvalue = raster::extract(r, cell_)) %>% 
       as.data.frame() 
    eresults[["tabularaster"]] <-    
       aggregate(result[,3:ncol(result)], by=list(result$object_), mean)       
})
5
  • The Q asks about "Extracting all bands to a polygon" which seems unclear. Do they want to polygonize each band of a raster? Do they want to crop to a polygon?
    – Spacedman
    Apr 14, 2021 at 16:04
  • @Spacedman I am going on "extract the values of a multiband raster to a polygon" but, you are correct, it is very unclear exactly what the OP is after here. Apr 14, 2021 at 16:09
  • @JeffreyEvans and Spacedman Thank you for all the help. I'm fairly new to R so please excuse my poor communication of what I'm trying to accomplish. Spacedman got it here: "I have polygons P and a multiband raster R and I want the mean of cells of each band R[[i]] inside each polygon P[[j]] including fractions of cells where the polygons cross cells". I'll update you on how this goes.
    – Tris
    Apr 19, 2021 at 0:52
  • what is the aim of [,-21] in exactextractr? Is it necessary? Thanks!
    – FraNut
    Feb 1 at 17:19
  • 1
    @FraNut The index is just dropping the column containing the fractional cell intersections. The index would be unique to each problem (number of bands in your data) and not necessary if you want to keep this data. If you wanted to automate it you could do something like [,-(nlayers(x)+1)] which would denote the last column. I also believe that in newer versions of exactextractr there is an argument to omit this column at the function call level. Feb 1 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.