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

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

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)

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!

CodePudding user response：

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