R - 与 SpatialPolygonsDataFrame 对象相交的 SpatialLinesDataFrame 列表的嵌套循环

我有一系列需要完成的步骤SpatialLinesDataFrame（此处的“线”）基于对象与多特征中各个特征的关系SpatialPolygonsDataFrame（“多边形”）对象。简而言之，每个线列表元素源自单个面要素内部，并且可能会也可能不会穿过一个或多个其他面要素。我想更新每个线元素，以将原点多边形连接到与线元素相交的每个单独多边形的第一个接触点。因此，每个线元素可能会成为多个新的线要素（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...

所以现在将所有这些包装到一个函数中并在您的线路上循环并完成工作！？