Search code examples
rraster

Calculate zonal statistics in R as in GIS


I have multiple rasters in a folder. I need to extract mean of each of these rasters over a polygon shape file (has more 2500 polygons).

I came across two functions zonal and extract. It says extract can be used for points, lines and polygons too. Is it the only difference ? (Yes/No expected)

How can I extract mean from these multiple rasters and specify different column names as per their filenames for these extracted mean values ?

Edit::

I found a code somewhere and implemented it. But it is taking forever and no progress at all.

grids <- list.files("my_path", pattern = "*.tif$")

#check the number of files in the raster list (grids)
length <- length(grids)

#read-in the polygon shapefile
poly <- readShapePoly("my_path/supplimentY.shp")

#create a raster stack
s <- stack(paste0("my_path/", grids))

#extract raster cell count (sum) within each polygon area (poly)
for (i in 1:length(grids)){
  ex <- extract(s, poly, fun='mean', na.rm=TRUE, df=TRUE, weights = TRUE)
# the code doesnot progress from here onwards. 
# i checked it by adding this line:: print(i)

}

#write to a data frame
dfr <- data.frame(ex)

Solution

  • You do not need the loop (you repeat the same operation at each iteration!).

    It should be like this:

    library(terra)   
    ff <- list.files("my_path", pattern = "\\.tif$", full=TRUE)
    s <- rast(ff)
    
    poly <- vect("my_path/supplimentY.shp")
    ex <- extract(s, poly, fun='mean', na.rm=TRUE)
    #or  ex <- extract(s, poly, fun='mean', na.rm=TRUE, exact=TRUE)
    

    You can also use terra::zonal

    With the (now obsolete) "raster" package you can do

    library(raster)   
    ff <- list.files("my_path", pattern = "\\.tif$", full=TRUE)
    s <- stack(ff)
    poly <- shapefile("my_path/supplimentY.shp")
    ex <- extract(s, poly, fun='mean', na.rm=TRUE, df=TRUE, weights = TRUE)