Мне нужно выполнить итерацию по многофункциональному элементу SpatialPolygonsDataFrame
(в данном документе SPDF) и стереть, где каждый многоугольник пересекается с пространственными линиями, содержащимися в списке однофункциональных SpatialLinesDataFrames
(SLDF), и сохранить обновленные «стертые» SLDF вновый список.Если линия пересекается с двумя различными полигональными объектами, я хочу, чтобы два обновленных «стертых» SLDF были созданы и добавлены в новый список SLDF.В приведенном ниже примере SLDF пересекаются ровно с одним SPDF, за исключением того, что один из SLDF пересекается с двумя различными элементами полигона SPDF.Следовательно, обновленный список должен содержать дополнительный элемент SLDF.
Однако, когда я запускаю вложенный цикл for, результирующий «стертый» список SLDF содержит то же количество элементов, что и исходный список SLDF.Я думаю, что есть проблема с моей структурой цикла, но я не могу понять это.
library(rgdal)
library(raster)
library(rgeos)
library(sp)
library(gdistance)
#Reproducible example data prep:
#'RDCO Regional Parks' data can be downloaded here: https://data-
rdco.opendata.arcgis.com/datasets?
group_ids=1950175c56c24073bb5cef3900e19460
parks <- readOGR("/Users/rachelfield/Documents/UBC/Data/Regional
Districts/RDCO/RDCO_Regional_Parks/RDCO_Regional_Parks.shp")
#DEM data downloaded here: https://pub.data.gov.bc.ca/datasets/175624/82e/
(files for '082e14_w.dem')
dem <- raster("/path/to/example/data/082e14_w.dem")
demproj <- "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
#reproject parks to match dem
p <- spTransform(parks, demproj)
#subset of parks data to reduce for example
e <- extent(dem)
p_crop <- crop(p, e)
p_sub <- p_crop[p_crop@data$Shapearea > 100000,]
p_sub2 <- p_sub[p_sub@data$CommonName != "Mission Creek Greenway Regional
Park",]
#fyi I delete polygon [7,] because violates rules of my actual data
(polygons do not touch)
#create polygon centroids and convert to SpatialPointsDataFrame using data
from each 'origin' polygon data.frame
p_cent <- gCentroid(p_sub, byid = TRUE)
p_centdf <- SpatialPointsDataFrame(p_cent, data = data.frame(p_sub),
match.ID = FALSE)
#manually create approx location of lowest elevation cell
lowest <- SpatialPoints(coords = cbind(-119.47,49.86), proj4string =
CRS("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"))
#find LCPs from 'origin' polygon centroids to lowest elevation cell
tr <- transition(dem, transitionFunction=function(x) 1/min(x),
directions=8)
trCost <- geoCorrection(tr, type="c")
#run shortestPath (LCP) analysis for all polygon centroids
lcp_list <- list()
for(i in 1:nrow(p_centdf)){
origin <- p_centdf[i,]
#find LCP from each centroid
lcp <- shortestPath(trCost, origin, goal = lowest,
output="SpatialLines")
#convert LCP SpatialLines to SpatialLinesDataFrame object type, and
preserve ID from original centroid
lcp_list[[i]] <- SpatialLinesDataFrame(lcp, data =
data.frame(p_centdf[i,]), match.ID = FALSE)
}
#my nested for-loop attempt to get resulting SLDF list
line_erasel <- list()
#iterate thru all SPDF features
for (i in seq_along(p_sub2)) {
#iterate thru all SLDF list elements
for (j in seq_along(lcp_list)) {
#if a SLDF intersects with a SPDF feature, execute erase function
if (tryCatch(!is.null(raster::intersect(lcp_list[[j]], p_sub2[i,])),
error=function(e) return(FALSE)) == 'TRUE'){
#erase part of line overlapped by intersected polygon and add to new list
line_erasel[[i]] <- erase(lcp_list[[j]],p_sub2[i,])
}}