This has been simplified in terra version 1.7.21. That is currently the development version. You should be able to install that version with
install.packages('terra', repos='https://rspatial.r-universe.dev')
Example data
library(terra)
#terra 1.7.21
v <- vect(system.file("ex/lux.shp", package="terra"))
r <- rast(system.file("ex/elev.tif", package="terra"))
r <- round((r-50)/100)
levels(r) <- data.frame(id=1:5, name=c("forest", "water", "urban", "crops", "grass"))
You can use extract
to get cell counts (use exact=TRUE
to consider cell fractions).
e <- extract(r, v, fun="table", na.rm=TRUE, exact=FALSE)
# add the regions label
data.frame(NAME_2=v$NAME_2[e[,1]], e[,-1])
# NAME_2 forest water urban crops grass
#1 Clervaux 0 0 28 459 74
#2 Diekirch 2 123 200 66 3
#3 Redange 0 109 170 170 17
#4 Vianden 0 31 39 57 3
#5 Wiltz 0 1 161 303 8
#6 Echternach 14 76 233 1 0
#7 Remich 50 147 24 0 0
#8 Grevenmacher 19 221 137 2 0
#9 Capellen 0 25 305 0 0
#10 Esch-sur-Alzette 0 190 229 15 0
#11 Luxembourg 0 184 223 16 0
#12 Mersch 0 167 239 14 0
To compute areas, you can use expanse
or zonal
.
First with expanse
:
# rasterize the zones
zone <- rasterize(v, r, "NAME_2")
expanse(r, unit="km", byValue=TRUE, zones=zone, wide=TRUE)[,-1]
# zone forest water urban crops grass
#1 Diekirch 1.110530 68.3111509 111.08929 36.621714 1.662192
#2 Echternach 7.784033 42.2720014 129.61202 0.555644 0.000000
#3 Grevenmacher 10.591752 123.1872271 76.33934 1.114409 0.000000
#4 Remich 27.951955 82.1601534 13.41321 0.000000 0.000000
#5 Capellen 0.000000 13.9496106 170.14699 0.000000 0.000000
#8 Esch-sur-Alzette 0.000000 106.2479147 128.04839 8.392079 0.000000
#10 Luxembourg 0.000000 102.7248564 124.44157 8.921887 0.000000
#11 Mersch 0.000000 92.9318291 133.02760 7.792452 0.000000
#12 Redange 0.000000 60.6687117 94.55675 94.459184 9.444716
#14 Vianden 0.000000 17.1997610 21.63354 31.602716 1.662192
#15 Wiltz 0.000000 0.5546965 89.26629 168.018372 4.437191
#17 Clervaux 0.000000 0.0000000 15.49914 253.829542 40.900555
Now with zonal
x <- cellSize(r, unit="km")
zonal(x, c(r, zone), fun="sum", wide=TRUE)
# group water urban crops grass forest
#1 Capellen 13.9496106 170.14699 0.000000 0.000000 0.000000
#3 Clervaux 0.0000000 15.49914 253.829542 40.900555 0.000000
#6 Diekirch 68.3111509 111.08929 36.621714 1.662192 1.110530
#11 Echternach 42.2720014 129.61202 0.555644 0.000000 7.784033
#15 Esch-sur-Alzette 106.2479147 128.04839 8.392079 0.000000 0.000000
#18 Grevenmacher 123.1872271 76.33934 1.114409 0.000000 10.591752
#22 Luxembourg 102.7248564 124.44157 8.921887 0.000000 0.000000
#25 Mersch 92.9318291 133.02760 7.792452 0.000000 0.000000
#28 Redange 60.6687117 94.55675 94.459184 9.444716 0.000000
#32 Remich 82.1601534 13.41321 0.000000 0.000000 27.951955
#35 Vianden 17.1997610 21.63354 31.602716 1.662192 0.000000
#39 Wiltz 0.5546965 89.26629 168.018372 4.437191 0.000000
To get the ("exact") area with extract
, you could do
e <- extract(c(x,r), v, exact=TRUE, na.rm=TRUE) |> na.omit()
e$area <- e$area * e$fraction
a <- aggregate(e[, "area", drop=FALSE], e[, c("ID", "name")], sum)
a <- data.frame(NAME_2 = v$NAME_2[a[,1]], a)[,-2]
reshape(a, idvar="NAME_2", timevar="name", direction="wide")
# NAME_2 area.forest area.water area.urban area.crops area.grass
#1 Diekirch 1.110530 67.6791652 109.78461 37.3792518 1.835912
#2 Echternach 8.224001 42.6656508 130.98765 0.3979812 NA
#3 Remich 26.893107 81.6930623 13.35597 NA NA
#4 Grevenmacher 9.685963 122.5968312 74.52150 1.1144093 NA
#6 Redange NA 60.6745706 93.69050 93.9922323 9.444456
#7 Vianden NA 17.0706616 22.18677 30.8643367 1.454355
#8 Wiltz NA 0.7361318 89.07461 167.5479682 4.768619
#12 Capellen NA 14.7992764 169.55486 NA NA
#13 Esch-sur-Alzette NA 105.7853834 127.80048 8.2993843 NA
#14 Luxembourg NA 103.0913118 125.09963 8.9218868 NA
#15 Mersch NA 92.6554656 132.88209 7.7924516 NA
#16 Clervaux NA NA 13.88577 251.4770800 40.600391
raster::extract
function, it is notoriously slow and the package is slowly being replaced with terra. Instead, try theterra::extract
function or even betterexactextractr::exact_extract
from the exactextractr package which can manage memory in a way that can process large problems. Just note that your polygon vector data must be in an sf class.method="bilinear"
there? land cover classes shouldn't be interpolated (averaged) or you end up with classes like "3.25". Also, do you needsmall=TRUE
? You'll get approximately the right answer without it for a big increase in speed.