如何通过循环对多份野火烟雾多边形Shapefile执行空间交集分析？

需求说明 #

需识别2005-2019年期间，每日野火烟雾多边形（Shapefile）覆盖范围内的诊所（含经纬度、ID信息的CSV文件）。已实现单个烟雾多边形的处理逻辑，现需批量处理约5475份每日文件，最终输出包含覆盖日期与对应诊所ID的数据框。

单文件处理代码 #

hms191030 <- readOGR("/Users/jin/Documents/HMS/hms_smoke20191030.shp") 
clinics <- read_csv('/Users/jin/Documents/ESRD/Clinics_cont.csv') 
clinicsW <- subset(clinics, PHYS_STATE_NM=="WASHINGTON" | PHYS_STATE_NM=="OREGON" | PHYS_STATE_NM=="CALIFORNIA")
clinicsW_sf <- st_as_sf(clinicsW, coords = c("LONG", "LAT"),  crs = 4326) 

hms191030_sf <- st_as_sf(hms191030, crs = 4326) 
sf::sf_use_s2(FALSE) 
intersection <- st_intersection(x = hms191030_sf, y = clinicsW_sf)

批量处理尝试与错误分析 #

尝试批量处理2019年1月1日-5日的文件时，执行相交函数阶段报错：

尝试代码 #

# 1. 读取诊所数据
clinics <- read_csv('/Users/songhyeonjin/Documents/ESRD/Clinics_cont.csv')
clinicsW <- subset(clinics, PHYS_STATE_NM=="WASHINGTON" | PHYS_STATE_NM=="OREGON" | PHYS_STATE_NM=="CALIFORNIA")
clinicsW_sf <- st_as_sf(clinicsW, coords = c("LONG", "LAT"),  crs = 4326) 

# 2. 读取多个HMS shp文件并生成列表
setwd("/Users/songhyeonjin/Documents/HMS/Shapefile_yrly/2019sub")
shps <- dir(getwd(), recursive = TRUE,"*.shp")
for (shp in shps) assign(shp, readOGR(shp))
my.list <- lapply(paste0('hms_smoke',20190101:20190105,'.shp'), get)
for (h in my.list) {
  st_as_sf(h, coords = c("long", "lat"),  crs = 4326) 
}

# 3. 对列表应用相交函数
sf::sf_use_s2(FALSE)
my_fun <- function(b) {
  intersection <- st_intersection(x = b, y = clinicsW_sf)
  return(intersection)
}
lapply(my.list, my_fun)

错误信息 #

Error in UseMethod("st_intersection") : 
no applicable method for 'st_intersection' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialVector', 'SpatialPolygonsNULL')"

错误原因 #

my.list中的对象仍为readOGR生成的SpatialPolygonsDataFrame类型，而st_intersection仅支持sf格式对象；且循环中转换sf后未重新赋值回列表，导致列表未更新为sf格式。

修正后的可运行代码 #

方案1：直接读取为sf格式（推荐） #

# 示例诊所数据
ID <- c(101, 102, 103, 104, 105)
LAT <- c(33.78595, 34.26310, 44.64489, 46.19070, 47.20550)
LONG <- c(-118.1900, -119.2293, -123.1116, -123.8365, -123.7543)
STATE <- c("California", "California", "Oregon", "Oregon", "Washington")
clinicsW <- data.frame(ID, LAT, LONG, STATE)
clinicsW_sf <- st_as_sf(clinicsW, coords = c("LONG", "LAT"), crs = 4326)

# 读取批量shp文件为sf对象
setwd("/Users/songhyeonjin/Documents/HMS/Shapefile_yrly/2019sub")
shps <- dir(getwd(), recursive = TRUE,"*.shp")
for (shp in shps) assign(shp, st_read(shp))
my.list <- lapply(paste0('hms_smoke',20190101:20190105,'.shp'), get)

# 应用空间相交函数
sf::sf_use_s2(FALSE)
my_fun <- function(b) {
  intersection <- st_intersection(x = b, y = clinicsW_sf)
  return(intersection)
}
result_list <- lapply(my.list, my_fun)

方案2：转换现有列表为sf格式 #

# 1. 读取诊所数据
clinics <- read_csv('/Users/jin/Documents/ESRD/Clinics_cont.csv')
clinicsW <- subset(clinics, PHYS_STATE_NM=="WASHINGTON" | PHYS_STATE_NM=="OREGON" | PHYS_STATE_NM=="CALIFORNIA")
clinicsW_sf <- st_as_sf(clinicsW, coords = c("LONG", "LAT"),  crs = 4326) 

# 2. 读取并转换为sf列表
setwd("/Users/jin/Documents/HMS/Shapefile_yrly/2019sub")
shps <- dir(getwd(), recursive = TRUE,"*.shp")
my.list <- lapply(shps, function(shp) {
  spatial_obj <- readOGR(shp)
  st_as_sf(spatial_obj, crs = 4326) # 多边形无需指定coords参数
})

# 3. 应用相交函数
sf::sf_use_s2(FALSE)
newlist <- lapply(my.list, function(b) st_intersection(x=b, y=clinicsW_sf))

最终结果整理 #

将所有日期的匹配结果合并为一个包含日期的数据集：

# 提取文件名中的日期信息
dates <- gsub("hms_smoke(\\d{8})\\.shp", "\\1", shps)
# 合并结果并添加日期列
final_df <- do.call(rbind, Map(function(df, date) {
  if(nrow(df) > 0) df$date <- date
  df
}, newlist, dates))

# 保留核心列（诊所ID、覆盖日期）
final_df <- final_df[, c("ID", "date")]

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