我有一个分类土地利用栅格 (habitat_raster.tif) 和一个约 200 个城市的 shapefile (municipality_shapefile.shp)。我想获得每个城市中每种栖息地类型的面积。市镇由代码唯一标识。
栅格有 2 个类与基于水位的 ID 相关联（class_low，class_high）。请参阅下面的属性表示例：
| ID | 低级 | class_high |
|---|---|---|
| 1 | 开放水域 | 开放水域 |
| 2 | 森林 | 水淹森林 |
| 3 | 草本 | 水淹草本植物 |
到目前为止，请参见以下代码：
#load in municipality as sf object
tmp_municipality_sf <-
st_read('spatial_analyses/data/municipality_shapefile.shp') %>%
# same projection string as habitat raster
st_transform(crs = '+proj=aea +lat_0=-6 +lon_0=-63 +lat_1=-13 +lat_2=3 +x_0=6377563.4 +y_0=6356256.9 +datum=WGS84 +units=m +no_defs')
# load raster data
habitat_raster <-
# CRS: Albers Equal Area, resolution (100 x 100):
terra::rast('data_cleaning/spatial/habitat_raster.tif')
# load raster attribute table
habitat_classes <-
read_csv(habitat_classes.csv)
# set habitat classes as levels
levels(habitat_raster) <- habitat_classes
# crop raster to municipality shapefile
habitat_crop <-
habitat_raster %>%
terra::crop(tmp_municipality_sf)
我试过使用 terra::extract()，但它为我提供了一个 ID 为 class_low 的数据框，并且它与市政代码无关。
municipality_extract <-
tmp_municipality_sf %>%
group_by(code) %>%
terra::vect() %>%
terra::extract(habitat_crop,
.,
mask = TRUE)
我还尝试将这些值放入市政府的空间数据框中：
municipaltiy_extract_df <-
tmp_municipality_sf %>%
mutate(habitats = tmp_municipality_sf %>%
terra::vect() %>%
terra::extract(habitat_crop,
.) %>%
# pulls a vector of values for each municipality
pull())
我不确定我是否可以改用 terra::zonal() 并将每个自治市形状用作“区域”。 但我想要的是一个数据框，其中包含市政代码和列中每个栖息地类别 (km^2) 的面积。我认为最好的方法是从每个栖息地类别中获取像素数，然后计算面积？
| mun_code | ID_1 | ID_2 | ID_3 |
|---|---|---|---|
| 100 | 10000 | 5000 | 2000 |
| 101 | 9500 | 4000 | 1500 |
示例数据
library(terra)
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"))
按区域计算每个类别的单元格数量
e <- extract(r, v)
e <- tapply(e[,2], e[,1], table)
e <- do.call(rbind, e)
d <- data.frame(NAME_2 = v$NAME_2, e)
d
# 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
总结每个区域的每个类别覆盖的面积
# compute the area of each cell
a <- cellSize(rast(r[[1]]), unit="km")
x <- c(r, a)
# extract and aggregate
f <- extract(x, v)
f <- aggregate(f[,3], f[,1:2], sum)
names(f)[3] <- "area"
f <- cbind(zone=v$NAME_2[f$ID], f[,-1])
head(f)
# zone name area
#1 Diekirch forest 1.110530
#2 Echternach forest 7.784033
#3 Remich forest 27.951955
#4 Grevenmacher forest 10.591752
#5 Diekirch water 68.311151
#6 Redange water 60.668712
重塑为宽格式
f <- reshape(f, idvar="zone", timevar="name", direction="wide")
f <- round(f, 1)
f[is.na(f)] <- 0
f <- f[order(f$zone), ]
f
zone area.forest area.water area.urban area.crops area.grass
#12 Capellen 0.0 13.9 170.1 0.0 0.0
#16 Clervaux 0.0 0.0 15.5 253.8 40.9
#1 Diekirch 1.1 68.3 111.1 36.6 1.7
#2 Echternach 7.8 42.3 129.6 0.6 0.0
#13 Esch-sur-Alzette 0.0 106.2 128.0 8.4 0.0
#4 Grevenmacher 10.6 123.2 76.3 1.1 0.0
#14 Luxembourg 0.0 102.7 124.4 8.9 0.0
#15 Mersch 0.0 92.9 133.0 7.8 0.0
#6 Redange 0.0 60.7 94.6 94.5 9.4
#3 Remich 28.0 82.2 13.4 0.0 0.0
#7 Vianden 0.0 17.2 21.6 31.6 1.7
#8 Wiltz 0.0 0.6 89.3 168.0 4.4
你也可以用
zon alrv <- rasterize(v, r, 1:nrow(v))
x <- c(rv, r)
levels(x) <- NULL
names(x) <- c("zone", "lulc")
a <- cellSize(r, unit="km")
zones <- unique(x, as.raster=TRUE)
z <- zonal(a, zones, sum)
g <- na.omit(cats(zones)[[1]])
m <- merge(g, z, by="label")[,3:5]
m$zone <- v$NAME_2[m$zone]
m$lulc <- levels(r)[[1]][m$lulc,2]
head(m)
# zone lulc area
#1 Clervaux urban 15.499144
#2 Clervaux crops 253.829542
#3 Clervaux grass 40.900555
#4 Esch-sur-Alzette water 106.247915
#5 Esch-sur-Alzette urban 128.048387
#6 Esch-sur-Alzette crops 8.392079