R/GIS: Find orthogonal distance between a location and nearest line
If I understand correctly, you can do this simply enough with gDistance
in the rgeos
package.
Read in the lines as SpatialLines/DataFrame
and points as SpatialPoints/DataFrame
and then loop over each point calculating the distance each time:
require(rgeos)
## untested code
shortest.dists <- numeric(nrow(sp.pts))
for (i in seq_len(nrow(sp.pts)) {
shortest.dists[i] <- gDistance(sp.pts[i,], sp.lns)
}
Here sp.pts
is the Spatial points object, and sp.lns
is the Spatial lines object.
You must loop so that you only compare a single coordinate in sp.pts
with the entirety of all lines geometries in sp.lns
, otherwise you get the distance from an aggregate value across all points.
Since your data are in latitude/longitude you should transform both the lines and points to a suitable projection since the gDistance
function assumes Cartesian distance.
MORE DISCUSSION AND EXAMPLE (edit)
It would be neat to get the nearest point on the line/s rather than just the distance, but this opens another option which is whether you need the nearest coordinate along a line, or an actual intersection with a line segment that is closer than any existing vertex. If your vertices are dense enough that the difference doesn't matter, then use spDistsN1
in the sp
package. You'd have to extract all the coordinates from every line in the set (not hard, but a bit ugly) and then loop over each point of interest calculating the distance to the line vertices - then you can find which is the shortest and select that coordinate from the set of vertices, so you can have the distance and the coordinate easily. There's no need to project either since the function can use ellipsoidal distances with longlat = TRUE
argument.
library(maptools)
## simple global data set, which we coerce to Lines
data(wrld_simpl)
wrld_lines <- as(wrld_simpl, "SpatialLinesDataFrame")
## get every coordinate as a simple matrix (scary but quick)
wrld_coords <- do.call("rbind", lapply(wrld_lines@lines, function(x1) do.call("rbind", lapply(x1@Lines, function(x2) x2@coords[-nrow(x2@coords), ]))))
Check it out interactively, you'll have to modify this to save the coords or minimum distances. This will plot up the lines and wait for you to click anywhere in the plot, then it will draw a line from your click to the nearest vertex on a line.
## no out of bounds clicking . . .
par(mar = c(0, 0, 0, 0), xaxs = "i", yaxs = "i")
plot(wrld_lines, asp = "")
n <- 5
for (i in seq_len(n)) {
xy <- matrix(unlist(locator(1)), ncol = 2)
all.dists <- spDistsN1(wrld_coords, xy, longlat = TRUE)
min.index <- which.min(all.dists)
points(xy, pch = "X")
lines(rbind(xy, wrld_coords[min.index, , drop = FALSE]), col = "green", lwd = 2)
}
The geosphere
package has the dist2line function that does this for lon/lat data. It can use Spatial* objects or matrices.
line <- rbind(c(-180,-20), c(-150,-10), c(-140,55), c(10, 0), c(-140,-60))
pnts <- rbind(c(-170,0), c(-75,0), c(-70,-10), c(-80,20), c(-100,-50),
c(-100,-60), c(-100,-40), c(-100,-20), c(-100,-10), c(-100,0))
d <- dist2Line(pnts, line)
d
Illustration of the results
plot( makeLine(line), type='l')
points(line)
points(pnts, col='blue', pch=20)
points(d[,2], d[,3], col='red', pch='x')
for (i in 1:nrow(d)) lines(gcIntermediate(pnts[i,], d[i,2:3], 10), lwd=2)
Looks like this can be done in the sf
package using the st_distance
function.
You pass your two sf
objects to the function. Same issue as with the other solutions in that you need to iterate over your points so that the function calculates the distance between every point to every point on the roadways. Then take the minimum of the resulting vector for the shortest distance.
# Solution for one point
min(st_distance(roads_sf, points_sf[1, ]))
# Iterate over all points using sapply
sapply(1:nrow(points_sf), function(x) min(st_distance(roads_sf, points_sf[x, ])))