使用R terra包提取各市镇生境类型面积的技术问题

你提到的需求完全可以用terra包实现，下面提供两种可靠的方法，均通过像素计数转换为面积：

方法一：extract + 分组统计 #

之前用extract没关联到市镇code，是因为没把提取结果的多边形ID和矢量的code对应起来。按以下步骤操作：

library(terra)
library(dplyr)
library(tidyr)

# 加载数据（你已完成，这里仅示意）
habitat <- rast("habitat_raster.tif")
munis <- vect("municipality_shapefile.shp")

# 提取栅格值，df=TRUE返回带多边形ID的数据框
extract_df <- extract(habitat, munis, df = TRUE, na.rm = TRUE)

# 给市镇矢量添加行ID，用于关联提取结果
munis$row_id <- seq_len(nrow(munis))

# 关联提取结果和市镇code
combined_df <- merge(extract_df, munis[, c("row_id", "code")], by.x = "ID", by.y = "row_id")

# 统计每个市镇各生境的像素数
pixel_counts <- combined_df %>%
  group_by(code, habitat) %>%
  summarise(count = n(), .groups = "drop")

# 计算单个像素的面积（单位：km²，假设栅格分辨率为米）
pixel_size_m <- res(habitat)[1]
pixel_area_km2 <- (pixel_size_m / 1000) ^ 2

# 转换为面积并整理成宽格式（每个市镇一行，两类生境各一列）
final_area <- pixel_counts %>%
  mutate(area_km2 = count * pixel_area_km2) %>%
  pivot_wider(
    names_from = habitat,
    values_from = area_km2,
    names_prefix = "area_",
    values_fill = 0  # 没有对应生境的市镇填充0
  )

# 重命名列匹配生境类型（根据你的栅格ID调整，比如class_low对应1，class_high对应2）
colnames(final_area) <- c("code", "area_low_water", "area_high_water")

方法二：zonal统计（更高效） #

terra::zonal()专门用于按分区统计栅格数据，需要先把市镇矢量转为和生境栅格同分辨率的分区栅格：

library(terra)

# 加载数据
habitat <- rast("habitat_raster.tif")
munis <- vect("municipality_shapefile.shp")

# 将市镇矢量转为栅格，每个像素赋值为对应code
muni_rast <- rasterize(munis, habitat, field = "code")

# 按分区统计各生境的像素数
zonal_stats <- zonal(habitat, muni_rast, fun = "table")

# 转换为数据框并计算面积
pixel_size_m <- res(habitat)[1]
pixel_area_km2 <- (pixel_size_m / 1000) ^ 2

final_area <- as.data.frame(zonal_stats)
colnames(final_area) <- c("code", "area_low_water", "area_high_water")
final_area[, -1] <- final_area[, -1] * pixel_area_km2

注意事项 #

• 确保你的市镇矢量code字段是唯一且非NA的，避免统计出错
• 栅格的类别值要和class_low、class_high对应，比如如果class_low的栅格ID是1，class_high是2，那么宽格式的列顺序要匹配
• 如果栅格有NoData值，记得在extract时加na.rm=TRUE，或者在zonal前用mask()过滤无效区域

内容的提问来源于stack exchange，提问作者megsruppUNBC