4

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()

figure example

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

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)

  • 1
    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
3

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

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.