R - вложенный цикл для списка SpatialLinesDataFrame, пересекаемого с объектами SpatialPolygonsDataFrame

Question

У меня есть ряд шагов, которые мне нужно выполнить для списка 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)
      }}}}}

Answer 1

Вот решение с использованием пакета sf.Я поработаю с созданными вами объектами и преобразую их в sf, хотя вы можете читать шейп-файлы в 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]])

Следующим шагом будет вычисление всех точек пересечения между вашей линией и полигонами.Вы должны: а) преобразовать многоугольники в линейные линии, в противном случае пересечение многоугольника и линии является линией, и б) преобразовать пересечения "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,]

Итак, что нам нужно, так это ближайшая точка пересечения в каждомгруппа точек, определенных одним и тем же идентификатором многоугольника.

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 сохраняет атрибуты многоугольников, которые пересекает линия, а также идентификатор линии:

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

, поэтому теперь оберните все это в функцию и зациклитеВаши строки и работа выполнены!?