如何使用回圈對多個多邊形形狀檔案使用交集函式？

我想識別 2005 年至 2019 年每天可用的野火煙霧多邊形 (shapefile) 下的診所（包含經度和緯度資訊以及診所 ID 的 csv 檔案）。

我可以使用下面的 R 代碼識別具有單個野火煙霧多邊形的診所。

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, coords = c("long", "lat"), crs = 4326) 
clinicsW_sf <- st_as_sf(clinicsW , coords = c("LONG", "LAT"), crs = 4326) 

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

由于我必須對多個每日野火煙霧多邊形（365 天 * 15 年）重復此程序，因此我嘗試將回圈函式與 shapefile 的子集（從 2019 年 1 月 1 日到 2019 年 1 月 5 日）一起使用。

# 1. Open clinics data
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. Open multiple HMS shapefiles and make a list
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. Apply intersection function to the list of SpatialPolygonsDataFrame
sf::sf_use_s2(FALSE)
my_fun <- function(b) {
  intersection <- st_intersection(x = b, y = clinicsW_sf)
  return(intersection)
}
lapply(my.list, my_fun)

然后對于行程 #3，我收到錯誤訊息

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

我期望得到的是一個資料框，其中顯示了 2005 年至 2019 年野火煙霧下的診所 ID。

這是我第一次在堆疊溢位中提問。在此先感謝您，非常感謝您的幫助。

[**這是帶有示例資料的更新代碼]

# Example data for the clinics location
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)

# Wildfire smoke shapefiles can be downlaoded here by selecting 'SMOKE' and 'Shapefile' and date range (2019 Jan 1 - 2019 Jan 5)
<https://www.ospo.noaa.gov/Products/land/hms.html#data>

# Edited code by John
setwd("/Users/songhyeonjin/Documents/HMS/Shapefile_yrly/2019sub")
shps <- dir(getwd(), recursive = TRUE,"*.shp")
# Just read in the shapefiles with sf to start. 
# Then you don't need the conversion.
for (shp in shps) assign(shp, st_read(shp))
my.list <-lapply( paste0('hms_smoke',20190101:20190105,'.shp'),get)

# 3. Apply intersection function to the list of sf objects
sf::sf_use_s2(FALSE)
my_fun <- function(b) {
  intersection <- st_intersection(x = b, y = clinicsW_sf)
  return(intersection)
}
lapply(my.list, my_fun)

所以這是根據約翰的回答作業的更新代碼。希望這對其他人也有幫助！

# 1. Open clinics data
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. Open multiple HMS shapefiles and make a list
setwd("/Users/jin/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. Apply intersection function to the list of SpatialPolygonsDataFrame
sf::sf_use_s2(FALSE)
newlist<-lapply(my.list, function(b) st_intersection(x=b, y=clinicsW_sf))

uj5u.com熱心網友回復：

我沒有你的資料，這讓這更難了。您從sp物件開始并將其中一些物件轉換為sf物件，而不是全部轉換。我認為問題出在您嘗試轉換的第二步。我認為您實際上并沒有轉換任何東西。你可能想要這樣的東西。

setwd("/Users/songhyeonjin/Documents/HMS/Shapefile_yrly/2019sub")
shps <- dir(getwd(), recursive = TRUE,"*.shp")
# Just read in the shapefiles with sf to start. 
# Then you don't need the conversion.
for (shp in shps) assign(shp, st_read(shp))
my.list <-lapply( paste0('hms_smoke',20190101:20190105,'.shp'),get)

# 3. Apply intersection function to the list of sf objects
sf::sf_use_s2(FALSE)
my_fun <- function(b) {
  intersection <- st_intersection(x = b, y = clinicsW_sf)
  return(intersection)
}
lapply(my_list, my_fun)

這是 QGIS 中檔案的螢屏截圖。由于診所位于目錄的底部，因此任何相交的煙霧形狀檔案都會出現在診所點的頂部。它們都不重疊，例如沒有交集。

不是說st_intersections不行。沒有交叉點，將回傳 0 行。

此時，所有代碼都可以作業。

標籤：r循环路口形状文件

上一篇：如何計算加權回圈？

下一篇：動態集的陣列和回圈實作？