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