0

I have a large SpatialLinesDataFrame "o" (~500,000 roads) and ~3000 survey points "r", I'd like to calculate the distance to the nearest line for each survey point and extract relevant columns.

My solution at present is to split the SpatialLinesDataFrame into chunks and calculate distance using gDistance in the "rgeos" package. This will take over 18 hours though on my machine.

Is there a faster way to do this, perhaps by converting lines to points first?

chunk <- 10000
n<-dim(o@data)[1]
m<-rep(1:ceiling(n/chunk),each=chunk)[1:n]

m.list<-vector("list", length=dim(r)[1])
for (i in 1:dim(r)[1]){
df.list<-vector("list", length=length(unique(m)))
for (j in 1:length(unique(m))){
p<-o[m==j,]
df.list[[j]]<-data.frame(ID=i,Dist=min(gDistance(r[i,],p,byid=TRUE)),p@data[which.min(gDistance(r[i,],p,byid=TRUE)),])}
d<-data.frame(do.call("rbind", df.list),stringsAsFactors=F)
m.list[[i]]<-d[which.min(d$Dist),]
}
  • Converting lines to vertex points won't work - a survey point could be a long way from two line vertex points but the line segment between the vertex points could go straight through it. Try with the sf spatial classes and st_distance function. – Spacedman Jan 30 at 17:00
1

I'm not sure that this is the fastest way to achieve what you are after, but I am fairly sure it is faster than your current implementation. As @spacedman recommends in his comment using the sf package is usually faster than using sp. Obviously I cannot verify my claim as your example is not reproducible. Here I present a reproducible example using some built-in data from the mapview package. These data are small compared to the data you refer to.

To simply append the attributes of the nearest line to each point we can do:

library(sf)
library(mapview)

trls = trails
trls$id = 1:nrow(trls) # make sure to have unique id to trace selected features later
brew = st_transform(breweries, st_crs(trls)) # make sure crs is the same

brew_w_nearest_trail = st_join(brew, trls, join = st_nearest_feature)

In case you also need the actual distances we can calculate the pairwise distances of the above identified nearest line features for each point:

dists = st_distance(brew, trls[brew_w_nearest_trail$id, ], by_element = TRUE)

To verify:

head(dists)
Units: [m]
[1] 2397.932533  481.878180    7.261746  453.968063   98.623719    6.525047

So the third feature pair should be about 7m apart. A visual check with mapview confirms:

mapview(brew[3, ]) + trls[brew_w_nearest_trail$id[3], ]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.