# Calculate distance of points - SpatialPoints object - along a path - SpatialLines object - using R

I have a SpatialPoints object and a SpatialLines object. I want to calculate all the distances from the start point (first coordinate of the SpatialLines object) to every point in the SpatialPoints object along the line.

Here is a code example:

``````# Load packages
library('sp')

# Define projection
epsg.32721 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

x <- c(788690, 788722, 788824, 788833)
y <- c(6193235, 6193202, 6193217, 6193224)

# Create SpatialPoints object
spatialPointsObject <- SpatialPoints(coords = cbind(x,y), proj4string = CRS(epsg.32721))

x2 <- c(788665.0, 788671.0, 788671.0, 788677.0, 788687.0, 788697.0, 788699.0, 788720.0, 788717.7, 788714.0, 788721.6, 788731.0, 788731.4, 788745.0, 788755.0, 788755.7, 788764.0, 788765.6, 788774.0, 788774.6, 788784.0, 788784.4, 788796.0, 788796.9, 788803.6, 788805.0, 788815.0, 788820.9, 788829.0, 788842.0)
y2 <- c(6193252, 6193244, 6193244, 6193236, 6193235, 6193235, 6193235, 6193230, 6193224, 6193215, 6193202, 6193202, 6193202, 6193216, 6193217, 6193217, 6193217, 6193216, 6193212, 6193212, 6193212, 6193212, 6193197, 6193198, 6193205, 6193207, 6193211, 6193215, 6193221, 6193230)

# Create SpatialLines object
spatialLinesObject <- SpatialLines(LinesList = list(Lines(slinelist = list(Line(coords = cbind(x2,y2))), ID = "1")), proj4string = CRS(epsg.32721))

# Plot SpatialPoints + SpatialLines
plot(spatialLinesObject, xlab = "Longitude", ylab = "Latitude", main = "SpatialPoints + SpatialLines")
plot(spatialPointsObject, col = "blue", pch = 19, add = TRUE)
plot(SpatialPoints(as.data.frame(coordinates(spatialLinesObject))[1,]), pch = 19, col = "red", add = TRUE)
text(as.data.frame(coordinates(spatialLinesObject))[1,], labels = "Start", cex = 0.75, pos = 3)
text(coordinates(spatialPointsObject), labels = seq_along(spatialPointsObject), cex = 0.75, pos = 3)
box()
``````

Is there a more efficient function in terms of RAM memory usage to use with a big data set?

This function calculates the distance from the startpoint to an intersecting point:

``````# Load packages
library('sp')
library('rgeos')

# Function to calculate the distance between startpoint of a line and an intersecting point
# The buffer can be used to make points intersecting even if their position is not exactly on the line
lengthToPoint = function(spatialLinesObj, intersectionPoint, buffer = 0.001){
intersectionPointCoordinates = intersectionPoint@coords
bufferedIntersectionPoint = gBuffer(spgeom = intersectionPoint, width = buffer)
# return 0 if there is no intersection at all
if(!gIntersects(spatialLinesObj, bufferedIntersectionPoint)){
print("The line does not intersect with the (buffered) point!")
return(0)
}
lineCoordinates = spatialLinesObj@lines[[1]]@Lines[[1]]@coords
numOfCoordinates = length(lineCoordinates[,1])
calculatedLength = 0
# split line into segments
for(i in 2:numOfCoordinates){
# create new spatiallines for the current segment and check if point is intersecting
currentLine = SpatialLines(LinesList = list(Lines(slinelist = list(Line(coords = lineCoordinates[(i-1):i,])), ID = "1")), spatialLinesObj@proj4string)
# no intersection
if(!gIntersects(currentLine, bufferedIntersectionPoint)){
calculatedLength = calculatedLength + gLength(spgeom = currentLine)
# intersection
} else {
# create line from start of current segment to intersection point
coordinates = matrix(data = c(lineCoordinates[i-1,], intersectionPointCoordinates), nrow = 2, byrow = T)
lastLine = SpatialLines(LinesList = list(Lines(slinelist = list(Line(coords = coordinates)), ID = "1")), spatialLinesObj@proj4string)
calculatedLength = calculatedLength + gLength(spgeom = lastLine)
print(paste("Found point on line segment", (i-1), "! Length from start point: ", calculatedLength))
return(calculatedLength)
}
}
}
``````

Test and Output:

``````firstPoint = spatialPointsObject[1]
secondPoint = spatialPointsObject[2]
fourthPoint = spatialPointsObject[4]

lengthToPoint(spatialLinesObject, firstPoint)
#Output:
> lengthToPoint(spatialLinesObject, firstPoint)
[1] "Found point on the line. Length from start point:  33.0498756211209"
[1] 33.04988

lengthToPoint(spatialLinesObject, secondPoint)
#Output:
> lengthToPoint(spatialLinesObject, secondPoint)
[1] "Found point on the line. Length from start point:  95.2520693698013"
[1] 95.25207

lengthToPoint(spatialLinesObject, fourthPoint)
# Output:
[1] "The line does not intersect with the (buffered) point!"
[1] 0
# -> change buffer parameter
lengthToPoint(spatialLinesObject, fourthPoint, buffer = 0.5)
# Output:
[1] "Found point on line segment 29 ! Length from start point:  230.227748237136"
[1] 230.2277
``````

This is just a quick solution that does not involve error handling. E.g. you might want to check for same projections.

EDIT: I adjusted the function to take a buffer parameter. See example `lengthToPoint(spatialLinesObject, fourthPoint)`

• The points are derived from the function `maptools::snapPointsToLines`, however the point 3 and point 4 aren't exactly over the line (distances from line: 0.238091 m and 0.1897367 m respectively). I found some problems using `snapPointsToLines` because sometimes the points are very close but not exactly over the line. So there is any intercept. I think I will need to add another parameter to the function `lengthToPoint` like width to use the function `gBuffer` inside and make the `gIntersects` between the currentLine and the buffer object. – Guzmán Sep 5 '16 at 18:44
• You are right. I adjusted the answer to include the buffer option. – Lars Sep 5 '16 at 20:04
• If you have a large dataset and you need to go parallel: `library(parallel) # N cores no_cores <- detectCores() - 1 # Initiate cluster cl <- makeCluster(no_cores) i <- 1:nrow(spatialPointsObject) buffer <- 0.001 # Objects and libraries to export to clusters clusterExport(cl, c("i", "spatialLinesObject", "spatialPointsObject", "buffer", "lengthToPoint")) clusterEvalQ(cl, c(library(sp), library(rgeos))) tryCatch(models <- parSapply(cl, i, function(i) lengthToPoint(intersectionPoint = spatialPointsObject[i,])), error = function(e) print(e)) # Finish stopCluster(cl)` – Guzmán Sep 7 '16 at 20:39

The GEOS library has some linear referencing functions (`GEOSInterpolate`/`GEOSProject`). I've wrapped these functions in the `rgeos` package (`gInterpolate`/`gProject`), since version 0-3.20 you can use the following:

``````# get distances of closest points along line
R> d <- gProject(spatialLinesObject, spatialPointsObject, normalized=FALSE)
R> d
# [1]  33.05  95.25 218.83 230.22

# coordinates of closest points (from distances)
R> gInterpolate(spatialLinesObject, d, normalized=FALSE)
# SpatialPoints:
#        x       y
# 1 788690 6193235
# 2 788722 6193202
# 3 788824 6193217
# 4 788833 6193224
# Coordinate Reference System (CRS) arguments: +proj=utm +zone=21 +south
# +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
``````
• Super cool solution, it's work fine for me! – Isaque Daniel Jul 28 '18 at 3:18