Search code examples
rcalculatorpolygonrasterterra

How to Zonal Statistic NDVI (tif) from land cover(shp) with multiple features in R

I have many NDVI tifs and a land cover shp with multiple features.

enter image description here

I want to know average NDVI in forest or grass and other land cover. I have tried terra::extract, but it extract by the CODE1 feature,not land TYPE

u <-vect("D:/muusVutm/type_49N.shp")
head(u)
f <- rast("D:/ndvi2001.tif")
f %>% terra::extract(u,fun='mean',na.rm=T)

enter image description here

In ArcGIS Pro, there is a tool called Zonal Statistics as Table, it could choose which feature of shp can be zonal statistics. But I have too many tifs(may be thousands), it couldn't be use this tool one by one.

Thanks for any help!

Solution

  • Here are three ways to accomplish that.

    Example data. We have raster r and polygons v and want the mean values for variable "NAME_1"; but the values of "NAME_1" are not unique (the same name is associated with multiple polygons).

    library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
v <- vect(system.file("ex/lux.shp", package="terra"))
v$NAME_1
# [1] "Diekirch"     "Diekirch"     "Diekirch"     "Diekirch"     "Diekirch"     "Grevenmacher"
# [7] "Grevenmacher" "Grevenmacher" "Luxembourg"   "Luxembourg"   "Luxembourg"   "Luxembourg"

    We could first aggregate the polygons and then use extract

    va <- aggregate(v, "NAME_1")
x <- extract(r, va, "mean", na.rm=TRUE, ID=FALSE)
cbind(NAME_1=va$NAME_1, x)
#        NAME_1 elevation
#1     Diekirch  403.1779
#2 Grevenmacher  283.8853
#3   Luxembourg  316.1935

    Or extract and then aggregate

    e <- extract(r, v)
aggregate(e[,-1,drop=FALSE], data.frame(NAME_1=v$NAME_1[e[,1]]), mean, na.rm=TRUE)
#        NAME_1 elevation
#1     Diekirch  403.1779
#2 Grevenmacher  283.8853
#3   Luxembourg  316.1935

    You can also use zonal

    z <- rasterize(v, r, "NAME_1")
zonal(r, z, mean, na.rm=TRUE)
#        NAME_1 elevation
#1     Diekirch  403.1779
#2 Grevenmacher  283.8853
#3   Luxembourg  316.1935

    In all cases, you can summarize the values of multiple rasters that have the same geometry in one step. For example

    rr <- c(r, r/2, r/10)
zonal(rr, z, mean, na.rm=TRUE)
#        NAME_1 elevation elevation elevation
#1     Diekirch  403.1779  201.5889  40.31779
#2 Grevenmacher  283.8853  141.9426  28.38853
#3   Luxembourg  316.1935  158.0968  31.61935