3

I am trying to extract zonal statistics for different classes in R. I have a raster with two classes (0,1) and need the area (or percentage) of each class under the area of each polygon. I tried exactextractr::exact_extract but I don't manage to extract a frequency of each class. I can get the sum of class 1, but that wont tell me the proportions.

6

You can get the frequency of a single class by passing a custom summary function to exact_extract. For example, to get the fraction of pixels that have a value of 1, you could run:

exact_extract(rast, polys, function(value, fraction) {
     sum(fraction[value == 1]) / sum(fraction)
  })

If you have an arbitrary number of classes, here's a solution that will provide the frequencies of each class in each polygon. It does not require knowing the classes in advance, nor does it require loading all intersected pixels into memory.

library(dplyr)
library(tidyr)

freqs <- exact_extract(r, g, function(value, coverage_fraction) {
  data.frame(value = value,
             frac = coverage_fraction / sum(coverage_fraction)) %>%
    group_by(value) %>%
    summarize(freq = sum(frac), .groups = 'drop') %>%
    pivot_wider(names_from = 'value',
                names_prefix = 'freq_',
                values_from = 'freq')
}) %>% 
  mutate(across(starts_with('freq'), replace_na, 0))

Basically, we provide a function to exact_extract that returns a one-row data frame for each polygon, with a column containing the frequency of each class found in that polygon. Doing this with a callback (i.e., specifying the fun argument) is important, because otherwise R must store every pixel that intersects every polygon in memory at the same time. With the callback, these pixels are reduced to a frequency table as each polygon is processed. Internally, exact_extract uses dplyr::bind_rows to merge the data frames for each polygon, which handles the fact that not all classes are present in each data frame. However, it fills in NA for the frequency of missing classes, so we use a final mutate call to replace these with zero.

8
  • 1
    If you have NA values in your raster, you'd need to add na.rm = TRUE to the sum function in the numerator and possibly the denominator, depending on whether you want the fraction of defined pixels that have a value of 1 or the fraction of all pixels that have a value of 1. – dbaston Jul 21 '20 at 18:59
  • 2
    exact_extract(rast, polys, function(value, coverage_fraction) {sum(coverage_fraction[value == 1],na.rm = T) / sum(coverage_fraction,na.rm = T)}) worked. Thanks so much! – mace Jul 21 '20 at 19:01
  • 1
    Why is it not sum(fraction[value == 1]) / length(fraction)? If your vector is [0,1] then the sum would only be elements that are positive and not the length of the vector, which would always result in 1 and not the proportion of elements that are one. – Jeffrey Evans Jul 21 '20 at 19:59
  • 1
    I think we're after the fraction of the polygon area that is associated with a value of 1. fraction is the fraction of the cell that is covered by the polygon, so it ranges from zero to 1. length(fraction) would give you the number of cells that are touched by the polygon, while sum(fraction) would give you the total covered area of cells that are touched by the polygon. – dbaston Jul 21 '20 at 20:40
  • 1
    Where using fractional pixel area is quite clever, it is also something that should be made clear. It does, in fact, effect the results and in a quick-and-dirty evaluation of the differences on simple circular polygons the Root Mean Square Error was 0.058. That is high enough to be considered a systematic bias in an analysis. I would imagine that with complex irregular polygons of varying size the RMSE between absolute intersection verses fractional corrected proportion would be much higher. I would prefer setting a threshold representing amount of intersection and just drop pixels. – Jeffrey Evans Jul 21 '20 at 22:15
3

You can write a simple function using prop.table and table to return proportions of multiple classes. The catch is that you have to know what all of the classes are before hand so you can fix the number of expected classes.

Here is an example of what is going on.

Here we set our "known" classes and then set up a loop that randomly samples a vector of 1:10 (some values may be missing in a given iteration). We can take the know classes and create an empty factor in x and then calculate our class proportions. If a value is missing then the resulting freq is 0.

classes <- 1:10

p <- list()
  for(i in 1:10) {
    x <- sample(1:10, 10, replace=TRUE)
    p[[i]] <- as.data.frame(prop.table(table(factor(x, levels = classes))))
  }
p

Now we can expand this idea to zonal statistics using exact_extract.

Add libraries and create some example data

library(raster)
library(sp)  
library(sf)                                                                          
library(rgeos)
library(exactextractr)

r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))

  r[] <- rpois(ncell(r), lambda=1)

x <- gBuffer(sampleRandom(r, 10, na.rm = TRUE, sp = TRUE),  
             byid = TRUE, width = 1000) 
x@data <- data.frame(x@data, ID=paste0("poly", 1:nrow(x)))

plot(r)
  plot(x, add=TRUE)

Now we extract the data and use lapply to apply a function to the resulting list object. We create the known classes by using unique on the raster object. Because you have to read the raster into memory, this could be a real processing bottleneck though.

( e <- exact_extract(r, as(x, "sf")) )   

classes <- sort(unique(r[]))
cp <- lapply(e, FUN=function(x) { as.data.frame(prop.table(table(factor(x[,1], 
                          levels = classes))))} ) 
names(cp) <- x$ID   
cp 
                      

You can perform some fancy data wrangling to get a data.frame that relates back to your polygons using a simple for loop with transpose. I set up an empty data.frame first so I can populate it using a simple assignment.

props <- data.frame(matrix(vector(), length(cp), length(classes)+1,
                    dimnames=list(c(), c("ID", paste0("class_",classes)))))
  props$ID <- names(cp)  
for(i in 1:length(cp)){ props[i,][2:ncol(props)] <- t(cp[[i]][,2]) }
props
1
  • Thanks for the very elaborate answer! – mace Jul 23 '20 at 16:25

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.