Search code examples
rr-rasterterra

How to calculate zonal statistics class and polygon wise using terra R package?

I want to calculate class and polygon-wise zonal statistics using terra R package. I could able to calculate polygon-wise mean values using the following code:

library(terra)

#Read a raster file
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

#Read a polygon file
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)

plot(v, add = T)

#Reclassify the elevation raster
m_r <- c(-Inf, 200, 1,  
         200, 300, 2,
         300, 400, 3,
         400, 500, 4,
         500, Inf, 5)

rclmat <- matrix(m_r, ncol=3, byrow=TRUE)
rc <- terra::classify(r, rclmat, include.lowest=TRUE)
plot(rc)

#Polygon wise mean calculation
mean_elev <- terra::extract(r, v, mean, na.rm=TRUE, bind = T, 
                               exact = T)
as.data.frame(mean_elev)

which returns me

  ID_1       NAME_1 ID_2           NAME_2 AREA    POP elevation
1     1     Diekirch    1         Clervaux  312  18081  467.3792
2     1     Diekirch    2         Diekirch  218  32543  334.6856
3     1     Diekirch    3          Redange  259  18664  377.2070
4     1     Diekirch    4          Vianden   76   5163  372.2499
5     1     Diekirch    5            Wiltz  263  16735  418.7867
6     2 Grevenmacher    6       Echternach  188  18899  314.7698
7     2 Grevenmacher    7           Remich  129  22366  240.2105
8     2 Grevenmacher   12     Grevenmacher  210  29828  283.2307
9     3   Luxembourg    8         Capellen  185  48187  329.8955
10    3   Luxembourg    9 Esch-sur-Alzette  251 176820  310.3833
11    3   Luxembourg   10       Luxembourg  237 182607  314.0103
12    3   Luxembourg   11           Mersch  233  32112  313.5930

My expected output is to have elevation statistics according to class as well as polygon-wise like

enter image description here

Solution

  • Couldn't think of a more direct method, not sure how it will scale.

    library(terra)

# Load SpatRaster
r <- rast(system.file("ex/elev.tif", package = "terra"))

# Load SpatVector
v <- vect(system.file("ex/lux.shp", package = "terra"))

# Reclassify the elevation Spatraster
m_r <- c(-Inf, 200, 1,  
         200, 300, 2,
         300, 400, 3,
         400, 500, 4,
         500, Inf, 5)

rclmat <- matrix(m_r, ncol = 3, byrow = TRUE)
rc <- classify(r, rclmat, include.lowest = TRUE)

# Create empty list to hold results
results <- list()

# Return mean elevation by class for each polygon in v
results <- lapply(1:nrow(v), function(i) {
  
  # Get feature in v
  p <- v[i, ]
  
  # Mask r and rc using p
  r1 <- mask(r, p)
  rc1 <- mask(rc, p)
  
  # Calculate mean elevation by class in rc1
  mean_ele <- zonal(r1, rc1, fun = mean, na.rm = TRUE)
  
  # Create vector, return NAME_2 value, pad with NAs
  res <- c(p[[4]], rep(NA, 5))
  
  # Add mean elevation values to res
  res[mean_ele[,1] + 1] <- mean_ele[,2]
  
  return(res)
  
})

# Convert the list to a dataframe
final <- data.frame(do.call(rbind, results))

# Rename columns
colnames(final) <- c("NAME_2", paste0("elevation_class", 1:5))

final
#              NAME_2 elevation_class1 elevation_class2 elevation_class3 elevation_class4 elevation_class5
# 1          Clervaux               NA               NA         377.8378         464.6545         516.1067
# 2          Diekirch            198.6         260.5344         346.9874         442.2857            511.5
# 3           Redange               NA         282.6525         349.3186         455.7701         507.7647
# 4           Vianden              200         263.9444         343.6875         447.6154         512.6667
# 5             Wiltz               NA              292         364.8659         447.3684         509.1818
# 6        Echternach         181.5556         267.1579         340.1055              404               NA
# 7            Remich         173.4808         251.1879         321.7037               NA               NA
# 8      Grevenmacher         175.9615         266.1818         329.5625            401.5               NA
# 9          Capellen               NA         290.2857         334.8779               NA               NA
# 10 Esch-sur-Alzette               NA         280.5859         329.4641           412.25               NA
# 11       Luxembourg               NA         277.0495         337.9529          413.375               NA
# 12           Mersch               NA         263.3077         343.6127         406.7692               NA