R - как написать вложенный цикл for с различными пространственными данными. - PullRequest
1 голос
/ 01 апреля 2019

Мне нужно выполнить итерацию по многофункциональному элементу 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,])
}}
...