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

Question

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"

# Load data points
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))

# Load data lines
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?

Answer 1

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 

Answer 2

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)