R计算网格单元土地利用比例及binomial GLM实现方法咨询

• sf：矢量空间数据处理核心包，负责图层读取、坐标系转换、网格生成、空间叠加、面积计算全流程操作
• tidyr/dplyr：表格数据整理工具，负责面积统计、占比计算、格式转换，生成符合建模要求的输入表
• stats：R原生包，无需额外安装，直接用于拟合二项分布GLM模型
• 可选拓展包：pROC用于计算模型AUC判别精度，DHARMa用于GLM残差分布检验，ggplot2用于结果可视化

1. 数据导入与预处理 #

首先读入已利用、未利用两个polygon图层，统一投影坐标系（必须使用等积投影/适合研究区的投影坐标系，禁止直接用经纬度计算面积，否则结果会有数量级误差），合并图层并给二分类因变量赋值：1代表监测到猫头鹰鸣叫声的已利用单元，0代表未监测到鸣叫声的未利用单元。

library(sf)
library(dplyr)
library(tidyr)

# 读入两个矢量图层，替换为本地文件实际路径
used_poly <- st_read("owl_used_habitat.gpkg")
unused_poly <- st_read("owl_unused_habitat.gpkg")

# 添加二分类标签
used_poly$use_label <- 1
unused_poly$use_label <- 0

# 合并图层，将EPSG编号替换为研究区对应的投影坐标系编码（如对应UTM分带的EPSG码）
target_crs <- st_crs("EPSG:XXXX")
all_poly <- bind_rows(used_poly, unused_poly) %>%
  st_transform(crs = target_crs)

2. 生成研究区网格 #

根据研究的猫头鹰物种家域范围确定网格尺寸（如雕鸮类家域多在1km²以上，可设置500m-1000m边长的网格），基于研究区总边界生成规则网格，再通过空间叠加确定每个网格的归属类别：对落在两类区域边界的网格，取重叠面积最大的类别作为该网格的标签，过滤掉重叠度过低的边缘网格。

# 提取研究区整体边界
study_extent <- st_union(all_poly)

# 生成规则网格，cellsize单位与投影坐标系一致，设为500即代表500m边长
grid_layer <- st_make_grid(study_extent, cellsize = 500, square = TRUE) %>%
  st_sf() %>%
  mutate(grid_id = row_number())

# 计算网格和两类利用区的重叠面积，匹配网格标签
grid_inter <- st_intersection(grid_layer, all_poly) %>%
  mutate(intersect_area = st_area(.)) %>%
  # 每个网格仅保留重叠面积最大的类别
  group_by(grid_id) %>%
  filter(intersect_area == max(intersect_area)) %>%
  ungroup() %>%
  select(grid_id, use_label)

# 生成最终带标签的网格层，过滤无标签的无效网格
grid_labeled <- grid_layer %>%
  left_join(st_drop_geometry(grid_inter), by = "grid_id") %>%
  filter(!is.na(use_label)) %>%
  mutate(grid_total_area = st_area(.)) # 预存每个网格的总面积，后续计算占比使用

3. 计算每个网格的土地利用类型占比 #

将带标签的网格和土地利用图斑做空间叠加，统计每个网格内各土地利用类型的总面积，除以网格总面积得到各类型占比，最后转成宽表格式（每列对应一种土地利用类型的占比）作为模型输入。

# 提取所有图斑的土地利用类型字段
lu_poly <- all_poly %>% select(land_use_type)

# 网格和土地利用图斑叠加，计算单块图斑落在网格内的面积
lu_inter <- st_intersection(grid_labeled, lu_poly) %>%
  mutate(lu_part_area = st_area(.)) %>%
  st_drop_geometry()

# 统计各类型占比，转换为建模用宽表
model_input <- lu_inter %>%
  group_by(grid_id, land_use_type) %>%
  summarise(lu_total = sum(lu_part_area), .groups = "drop") %>%
  left_join(st_drop_geometry(grid_labeled %>% select(grid_id, use_label, grid_total_area)), by = "grid_id") %>%
  mutate(lu_ratio = as.numeric(lu_total / grid_total_area)) %>%
  select(grid_id, use_label, land_use_type, lu_ratio) %>%
  # 转宽表，网格内不存在的土地类型占比填0
  pivot_wider(names_from = land_use_type, values_from = lu_ratio, values_fill = 0) %>%
  select(-grid_id) # 删除ID列，仅保留建模变量

注意：所有土地利用类型的占比和为1，建模时需要留一类作为参照组，否则会出现完全共线性导致模型报错；如果存在总占比低于1%的稀有土地类型，建议提前合并，避免模型自由度不足。

4. 拟合二项GLM与模型检验 #

直接调用R原生glm()函数，指定分布族为二项分布、连接函数为logit即可，拟合完成后可开展常规模型检验。

# 拟合模型，公式中.代表用除use_label外的所有列作为自变量
owl_model <- glm(use_label ~ ., data = model_input, family = binomial(link = "logit"))

# 查看模型系数结果
summary(owl_model)

# 可选：模型检验
# 1. 残差检验（需提前加载DHARMa包）
# library(DHARMa)
# sim_res <- simulateResiduals(owl_model)
# plot(sim_res)

# 2. 判别精度计算（需提前加载pROC包）
# library(pROC)
# pred <- predict(owl_model, type = "response")
# auc(model_input$use_label, pred)

建模前建议先检查自变量相关性，将相关系数高于0.7的变量剔除其一，避免多重共线性影响结果可靠性；如果数据存在明显空间自相关，可后续改用混合效应模型加入空间随机效应优化。

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