Here we will conduct zonal statistics in R using the highly efficient R package, exactextractr. See more information and an excellent vignette here: https://isciences.gitlab.io/exactextractr/

We will calculate the mean terrestrial vertebrate richness per ADM 1 boundary in Peru for our example.

Attach packages

library(terra) # rasters 
library(sf) # polygons 
library(rnaturalearth) # spatial boundary data
library(rnaturalearthhires) # may have to use devtools::install_github("ropensci/rnaturalearthhires")
library(exactextractr) # zonal statistics
library(kableExtra) # html tables
library(RColorBrewer) # palettes
library(tidyverse) # everything else

Define raster

# example raster
# using species richness
rast <- rast("verts_richness.tif")/100

# define crs
my_crs <- crs(rast)

plot(rast)

Define polygons

# get some polygons from natural earth
# we'll be working w/ subnational boundaries
polys <- ne_states(
  country = "Peru",
  returnclass = "sf") %>% 
  dplyr::select(name_en) %>% 
  st_transform(my_crs) # ensure crs is aligned

ggplot() + 
  geom_sf(data = polys,
          aes(fill = name_en),
          size = 0.1) + 
  theme_void()

Crop raster to region of interest

rast_roi <- rast %>% 
  crop(vect(polys)) %>% # needs to be vect to work with rast 
  mask(vect(polys))

plot(rast_roi)

Calculate zonal statistics with exactextractr

# now summarize within boundaries with exact_extract
zonal <- exact_extract(
  x = rast_roi, # rast
  y = polys, # vect
  fun = "mean",
  append_cols = "name_en") %>%  # want to keep names 
  rename("mean_richness" = mean)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%
kable(zonal) %>% kable_styling()
name_en mean_richness
12 Tacna Region 30.03053
281 Madre de Dios Region 746.80774
1352 Loreto Region 698.93250
1356 Amazonas Region 412.35754
1358 Cajamarca Region 132.77730
1360 Tumbes region 109.76875
1362 Piura Region 81.85780
1365 Ucayali Region 732.71399
1370 Puno Region 207.97925
1804 Moquegua Region 33.13277
1805 Arequipa Region 38.58153
1806 Ica Region 27.71490
1807 Lima Province 64.62955
1808 Lima Province 43.44174
1809 Callao Province 39.81061
1810 Áncash Region 67.92921
1811 La Libertad Region 89.44556
1812 Lambayeque Region 59.38707
3567 Cusco Región 326.31274
3568 Ayacucho Region 88.72918
3569 Apurímac Region 81.15025
3570 Huancavelica Region 87.56619
3571 San Martín Region 527.91577
3572 Huanuco Region 357.24387
3573 Pasco Region 399.99280
3574 Junín Region 310.83301 
# now link back with polygons by name
zonal_polys <- polys %>% 
  left_join(zonal, by = "name_en")

Voila!

ggplot() + 
  geom_sf(
    data = zonal_polys,
    aes(fill = mean_richness),
    size = 0.1) + 
  scale_fill_distiller(palette = "YlGnBu", 
                       direction = 1) + 
  theme_void()