R-与SpatialPolygonsDataFrame对象相交的SpatialLinesDataFrame列表的嵌套循环 [英] R - nested loop for list of SpatialLinesDataFrame intersected with SpatialPolygonsDataFrame objects

问题描述

我需要根据对象与多功能SpatialPolygonsDataFrame(多边形")对象内各个要素的关系，在一系列SpatialLinesDataFrame(此处为线")对象上完成一系列步骤.简而言之，每个线列表元素都起源于单个面要素内，并且可能会也可能不会穿过一个或多个其他面要素.我想更新每个线元素以将原多边形连接到该线元素相交的每个单独多边形的第一接触点.因此，每个线元素可能会变成多个新的线要素(n =相交多边形的数量).

I have a series of steps I need to complete on a list of SpatialLinesDataFrame ('lines' herein) objects based on their relationships with individual features within a multi-feature SpatialPolygonsDataFrame ('polygons') object. In short, each line list element originates inside a single polygon feature, and may or may not pass through one or more other polygon features. I want to update each line element to connect origin polygons to the first point of contact for each individual polygon intersected by the line element. So, each line element may become multiple new line features (n=number of intersected polygons).

我想有效地做到这一点，因为我的线列表和面要素很多.我在下面提供了示例数据和逐步说明.我是R语言的新手，而不是程序员，所以我不知道我提出的任何建议是否有效.

I would like to do this efficiently as my lines lists and polygon features are numerous. I have provided example data and STEP-wise description of what I am trying to do below. I am new to R and not a programmer, so I don't know if any of what I propose is valid.

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包的解决方案.尽管您可以将shapefile读入sf对象并从头开始创建它们，但我将使用您创建的对象并将它们转换为sf，因此，如果没有其他理由使用sp对象，我建议您这样做.

Here's a solution using the sf package. I'll work with the objects you create and convert them to sf, although you can read shapefiles into sf objects and create them from scratch so if there's no other reason to use sp objects I'd recommend that.

使用这些软件包:

library(sf)
library(dplyr)

转换多边形.我从parks_sub删除了一些列，以便它可以打印得更整洁-如果需要，请不要删除它们:

Convert polygons. I'm dropping a load of columns from parks_sub just so it can print neater - if you need them don't drop them:

p = st_as_sf(parks_sub)
p = p[,c("OBJECTID","PARK_NAME")]

转换行.我将只使用您的第一行对象，在您的列表上编写一个循环以完成整个设置:

Convert lines. I'm only going to work with your first line object, write a loop over your list to do a whole set:

l1 = st_as_sf(lineldf[[1]])

下一步是计算线和面之间的所有交点.您必须:a)将多边形转换为线串，否则多边形和直线的交点为直线，并且b)使用st_cast:

Next step is to compute all the intersection points between your line and your polygons. You have to: a) convert the polygons to linestrings otherwise the intersection of a polygon and a line is a line, and b) convert the "MULTIPOINT" intersections when a line crosses a polygon more than once into a set of POINT objects using st_cast:

pts = st_cast(st_cast(st_intersection(l1,
             st_cast(p,"MULTILINESTRING")
             ),"MULTIPOINT"),"POINT")

接下来，我们需要该行的第一点.对于您在示例中创建的数据，多边形的线端实际上是第二点，因此我将在此处进行提取.

Next we need the first point of the line. For the data you create in the example, the line end in the polygon is actually the second point, so I'll extract that here.

first_point = st_cast(l1$geom,"POINT")[2]

如果在您的真实数据中第一点，则在其中放置[1].如果这取决于那那是另一个小问题.

If in your real data its the first point then put [1] in there. If it depends then that's another little problem.

现在计算从第一个点到所有相交点的距离:

Now compute the distances from that first point to all the intersection points:

pts$d_first = st_distance(first_point, pts)[1,]

所以我们现在想要的是由相同多边形ID定义的每组点中最接近的交点.

So what we want now is the nearest intersection point in each group of points defined by the same polygon ID.

near_pts = pts %>% group_by(OBJECTID)  %>% filter(d_first==min(d_first))

然后，从第一个点到最近的点构建所需的线:

Then the desired lines are constructed from the first point to those nearest points:

nlines = st_cast(st_union(near_pts, first_point),"LINESTRING")

以减小的宽度绘制多边形和线条以显示重叠:

Plot the polygons and the lines in decreasing width to show the overlap:

plot(p$geom)
plot(nlines$geom, lwd=c(10,7,4), col=c("black","grey","white"), add=TRUE)

请注意，三行包括从第一个点到多边形所处边界的短线(白色)-如果您不希望这样做，则可以滤除距数据框最近距离的点在构造线条之前-但假定第一个点在多边形内...

Note the three lines include a short one (in white) from the first point to the boundary of the polygon it is in - if you don't want this you can filter out the point with the nearest distance from the data frame before constructing the lines - but that assumes the first point is inside a polygon...

nlines保留线相交的多边形的属性以及线的ID:

nlines retains the attributes of the polygons the line intersects, as well as the ID of the line:

> 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...

现在将所有内容包装到一个函数中，并循环遍历您的行并完成工作！?

so now wrap all that into a function and loop that over your lines and job done!?

这篇关于R-与SpatialPolygonsDataFrame对象相交的SpatialLinesDataFrame列表的嵌套循环的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！