R - 与 SpatialPolygonsDataFrame 对象相交的 SpatialLinesDataFrame 列表的嵌套循环
R - nested loop for list of SpatialLinesDataFrame intersected with SpatialPolygonsDataFrame objects
我需要根据 SpatialLinesDataFrame（此处为 'lines'）对象的列表完成一系列步骤，这些步骤基于它们与多特征 SpatialPolygonsDataFrame 中各个特征的关系('polygons') 对象。简而言之，每个线列表元素都起源于单个多边形要素内部，并且可能穿过也可能不穿过一个或多个其他多边形要素。我想更新每个线元素以将原点多边形连接到与线元素相交的每个单独多边形的第一个接触点。因此，每个线元素可能成为多个新的线要素（n=相交多边形的数量）。
我想高效地执行此操作，因为我的线列表和多边形要素很多。我在下面提供了示例数据和我要执行的操作的逐步说明。我是 R 的新手，不是程序员，所以我不知道我提出的任何建议是否有效。
library(sp)
library(rgdal)
library(raster)
###example data prep START
#example 'RDCO Regional Parks' data can be downloaded here:
https://data-rdco.opendata.arcgis.com/datasets group_ids=1950175c56c24073bb5cef3900e19460
parks <- readOGR("/path/to/example/data/RDCO_Regional_Parks/RDCO_Regional_Parks.shp")
plot(parks)
#subset watersheds for example
parks_sub <- parks[parks@data$Shapearea > 400000,]
plot(parks_sub, col='green', axes = T)
#create SpatialLines from scratch
pts_line1 <- cbind(c(308000, 333000), c(5522000, 5530000))
line1 <- spLines(pts_line1, crs = crs(parks_sub))
plot(line1, axes=T, add=T) #origin polygon = polyl[[4]] = OBJECTID 181
pts_line2 <- cbind(c(308000, 325000), c(5524000, 5537000))
line2 <- spLines(pts_line2, crs = crs(parks_sub))
plot(line2, axes=T, add=T) #origin polygon = polyl[[8]] = OBJECTID 1838
linel <- list()
linel[[1]] <- line1
linel[[2]] <- line2
#convert to SpatialLinesDataFrame objects
lineldf <- lapply(1:length(linel), function(i) SpatialLinesDataFrame(linel[[i]], data.frame(id=rep(i, length(linel[[i]]))), match.ID = FALSE))
#match id field value with origin polygon
lineldf[[1]]@data$id <- 181
lineldf[[2]]@data$id <- 1838
###example data prep END
#initiate nested for loop
for (i in 1:length(lineldf)) {
for (j in 1:length(parks_sub[j,])) {
#STEP 1:for each line list feature (NB: with ID matching origin polygon ID)
#identify whether it intersects with a polygon list feature
if (tryCatch(!is.null(intersect(lineldf[[i]], parks_sub[j,])), error=function(e) return(FALSE)) == 'FALSE'){
next
}
#if 'FALSE', go on to check intersect with next polygon in list
#if 'TRUE', go to STEP 2
#STEP 2: add intersected polygon OBJECTID value to SLDF new column in attribute table
#i.e., deal with single intersected polygon at a time
else {
lineldf[[i]]@data["id.2"] = parks_sub[j,]@data$OBJECTID
#STEP 3: erase portion of line overlapped by intersected SPDF
line_erase <- erase(lineldf[[i]],parks_sub[j,])
#STEP 4: erase line feature(s) that no longer intersect with the origin polygon
#DO NOT KNOW HOW TO SELECT feature (i.e., line segment) within 'line_erase' object
if (tryCatch(!is.null(intersect(line_erase[???], parks_sub[j,])), error=function(e) return(FALSE)) == 'FALSE'){
line_erase[???] <- NULL}
else {
#STEP 5: erase line features that only intersect with origin polygon
if (line_erase[???]@data$id.2 = parks_sub[j,]@data$OBJECTID){
line_erase[???] <- NULL
}
else {
#STEP 6: write valid line files to folder
writeOGR(line_erase,
dsn = paste0("path/to/save/folder", i, ".shp"),
layer = "newline",
driver = 'ESRI Shapefile',
overwrite_layer = T)
}}}}}
这是一个使用 sf 包的解决方案。我将使用您创建的对象并将它们转换为 sf，尽管您可以将 shapefile 读入 sf 对象并从头开始创建它们，因此如果没有其他理由使用 sp 对象我会推荐。
使用这些包：
library(sf)
library(dplyr)
转换多边形。我从 parks_sub 中删除了一些列，以便打印更整洁 - 如果您需要它们，请不要删除它们：
p = st_as_sf(parks_sub)
p = p[,c("OBJECTID","PARK_NAME")]
转换行。我只打算使用你的第一行对象，在你的列表上写一个循环来做一整套：
l1 = st_as_sf(lineldf[[1]])
下一步是计算直线和多边形之间的所有交点。您必须：a) 将多边形转换为线串，否则多边形和一条线的交点是一条线，以及 b) 当一条线多次穿过多边形时将 "MULTIPOINT" 交点转换为一组 POINT 对象使用 st_cast:
pts = st_cast(st_cast(st_intersection(l1,
st_cast(p,"MULTILINESTRING")
),"MULTIPOINT"),"POINT")
接下来我们需要直线的第一个点。对于您在示例中创建的数据，多边形中的线端实际上是第二个点，因此我将在此处提取它。
first_point = st_cast(l1$geom,"POINT")[2]
如果在您的真实数据中它是第一个点，则将 [1] 放在那里。如果这取决于那是另一个小问题。
现在计算从第一个点到所有交点的距离：
pts$d_first = st_distance(first_point, pts)[1,]
所以我们现在要的是同一多边形ID定义的每组点中最近的交点。
near_pts = pts %>% group_by(OBJECTID) %>% filter(d_first==min(d_first))
然后从第一个点到最近的点构建所需的线：
nlines = st_cast(st_union(near_pts, first_point),"LINESTRING")
绘制宽度递减的多边形和线以显示重叠：
plot(p$geom)
plot(nlines$geom, lwd=c(10,7,4), col=c("black","grey","white"), add=TRUE)
请注意，这三行包括从第一个点到它所在的多边形边界的一条短线（白色）- 如果您不想要它，您可以过滤掉距离最近的点构建线之前的数据框 - 但假设第一个点在多边形内...
nlines保留线相交的多边形的属性，以及线的ID：
> nlines
Simple feature collection with 3 features and 4 fields
geometry type: LINESTRING
dimension: XY
bbox: xmin: 310276 ymin: 5522728 xmax: 333000 ymax: 5530000
epsg (SRID): 26911
proj4string: +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
id OBJECTID PARK_NAME d_first geometry
1 181 2254 Mission Creek 6794.694 m LINESTRING (326528.6 552792...
2 181 1831 Glen Canyon 23859.161 m LINESTRING (310276 5522728,...
3 181 1838 Black Mountain 1260.329 m LINESTRING (331799.6 552961...
所以现在将所有这些都包装到一个函数中并在你的代码中循环并完成工作！？
我需要根据 SpatialLinesDataFrame（此处为 'lines'）对象的列表完成一系列步骤，这些步骤基于它们与多特征 SpatialPolygonsDataFrame 中各个特征的关系('polygons') 对象。简而言之，每个线列表元素都起源于单个多边形要素内部，并且可能穿过也可能不穿过一个或多个其他多边形要素。我想更新每个线元素以将原点多边形连接到与线元素相交的每个单独多边形的第一个接触点。因此，每个线元素可能成为多个新的线要素（n=相交多边形的数量）。
我想高效地执行此操作，因为我的线列表和多边形要素很多。我在下面提供了示例数据和我要执行的操作的逐步说明。我是 R 的新手，不是程序员，所以我不知道我提出的任何建议是否有效。
library(sp)
library(rgdal)
library(raster)
###example data prep START
#example 'RDCO Regional Parks' data can be downloaded here:
https://data-rdco.opendata.arcgis.com/datasets group_ids=1950175c56c24073bb5cef3900e19460
parks <- readOGR("/path/to/example/data/RDCO_Regional_Parks/RDCO_Regional_Parks.shp")
plot(parks)
#subset watersheds for example
parks_sub <- parks[parks@data$Shapearea > 400000,]
plot(parks_sub, col='green', axes = T)
#create SpatialLines from scratch
pts_line1 <- cbind(c(308000, 333000), c(5522000, 5530000))
line1 <- spLines(pts_line1, crs = crs(parks_sub))
plot(line1, axes=T, add=T) #origin polygon = polyl[[4]] = OBJECTID 181
pts_line2 <- cbind(c(308000, 325000), c(5524000, 5537000))
line2 <- spLines(pts_line2, crs = crs(parks_sub))
plot(line2, axes=T, add=T) #origin polygon = polyl[[8]] = OBJECTID 1838
linel <- list()
linel[[1]] <- line1
linel[[2]] <- line2
#convert to SpatialLinesDataFrame objects
lineldf <- lapply(1:length(linel), function(i) SpatialLinesDataFrame(linel[[i]], data.frame(id=rep(i, length(linel[[i]]))), match.ID = FALSE))
#match id field value with origin polygon
lineldf[[1]]@data$id <- 181
lineldf[[2]]@data$id <- 1838
###example data prep END
#initiate nested for loop
for (i in 1:length(lineldf)) {
for (j in 1:length(parks_sub[j,])) {
#STEP 1:for each line list feature (NB: with ID matching origin polygon ID)
#identify whether it intersects with a polygon list feature
if (tryCatch(!is.null(intersect(lineldf[[i]], parks_sub[j,])), error=function(e) return(FALSE)) == 'FALSE'){
next
}
#if 'FALSE', go on to check intersect with next polygon in list
#if 'TRUE', go to STEP 2
#STEP 2: add intersected polygon OBJECTID value to SLDF new column in attribute table
#i.e., deal with single intersected polygon at a time
else {
lineldf[[i]]@data["id.2"] = parks_sub[j,]@data$OBJECTID
#STEP 3: erase portion of line overlapped by intersected SPDF
line_erase <- erase(lineldf[[i]],parks_sub[j,])
#STEP 4: erase line feature(s) that no longer intersect with the origin polygon
#DO NOT KNOW HOW TO SELECT feature (i.e., line segment) within 'line_erase' object
if (tryCatch(!is.null(intersect(line_erase[???], parks_sub[j,])), error=function(e) return(FALSE)) == 'FALSE'){
line_erase[???] <- NULL}
else {
#STEP 5: erase line features that only intersect with origin polygon
if (line_erase[???]@data$id.2 = parks_sub[j,]@data$OBJECTID){
line_erase[???] <- NULL
}
else {
#STEP 6: write valid line files to folder
writeOGR(line_erase,
dsn = paste0("path/to/save/folder", i, ".shp"),
layer = "newline",
driver = 'ESRI Shapefile',
overwrite_layer = T)
}}}}}
这是一个使用 sf 包的解决方案。我将使用您创建的对象并将它们转换为 sf，尽管您可以将 shapefile 读入 sf 对象并从头开始创建它们，因此如果没有其他理由使用 sp 对象我会推荐。
使用这些包：
library(sf)
library(dplyr)
转换多边形。我从 parks_sub 中删除了一些列，以便打印更整洁 - 如果您需要它们，请不要删除它们：
p = st_as_sf(parks_sub)
p = p[,c("OBJECTID","PARK_NAME")]
转换行。我只打算使用你的第一行对象，在你的列表上写一个循环来做一整套：
l1 = st_as_sf(lineldf[[1]])
下一步是计算直线和多边形之间的所有交点。您必须：a) 将多边形转换为线串，否则多边形和一条线的交点是一条线，以及 b) 当一条线多次穿过多边形时将 "MULTIPOINT" 交点转换为一组 POINT 对象使用 st_cast:
pts = st_cast(st_cast(st_intersection(l1,
st_cast(p,"MULTILINESTRING")
),"MULTIPOINT"),"POINT")
接下来我们需要直线的第一个点。对于您在示例中创建的数据，多边形中的线端实际上是第二个点，因此我将在此处提取它。
first_point = st_cast(l1$geom,"POINT")[2]
如果在您的真实数据中它是第一个点，则将 [1] 放在那里。如果这取决于那是另一个小问题。
现在计算从第一个点到所有交点的距离：
pts$d_first = st_distance(first_point, pts)[1,]
所以我们现在要的是同一多边形ID定义的每组点中最近的交点。
near_pts = pts %>% group_by(OBJECTID) %>% filter(d_first==min(d_first))
然后从第一个点到最近的点构建所需的线：
nlines = st_cast(st_union(near_pts, first_point),"LINESTRING")
绘制宽度递减的多边形和线以显示重叠：
plot(p$geom)
plot(nlines$geom, lwd=c(10,7,4), col=c("black","grey","white"), add=TRUE)
请注意，这三行包括从第一个点到它所在的多边形边界的一条短线（白色）- 如果您不想要它，您可以过滤掉距离最近的点构建线之前的数据框 - 但假设第一个点在多边形内...
nlines保留线相交的多边形的属性，以及线的ID：
> nlines
Simple feature collection with 3 features and 4 fields
geometry type: LINESTRING
dimension: XY
bbox: xmin: 310276 ymin: 5522728 xmax: 333000 ymax: 5530000
epsg (SRID): 26911
proj4string: +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
id OBJECTID PARK_NAME d_first geometry
1 181 2254 Mission Creek 6794.694 m LINESTRING (326528.6 552792...
2 181 1831 Glen Canyon 23859.161 m LINESTRING (310276 5522728,...
3 181 1838 Black Mountain 1260.329 m LINESTRING (331799.6 552961...
所以现在将所有这些都包装到一个函数中并在你的代码中循环并完成工作！？