r – Calculating Distance of Points Along a Path Using R

distancelinepointrsp

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?

Best Answer

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

Related Solutions

PostGIS Distance – Calculating Total Great Circle Distance Along a Path of Points

It is fairly straightforward to write a single query for this, you just need to break it down into logical blocks. To start with, you need your point geometry sorted by track and time:

SELECT track_id, gpstime, ST_AsText(the_geom) FROM points ORDER BY track_id, gpstime;

Here I assume you have your points in a geometry column. A geography column in this case would achieve the same thing, and actually save a bit of faff in the next step.

(Note here I've used ST_AsText() to make things readable. You should remove this in the final query).

With my simple database of 2 tracks and 6 points, that query gives me:

 track_id |            gpstime            |         st_astext
----------+-------------------------------+----------------------------
        1 | 2012-01-27 14:36:26.354665+00 | POINT(-71.064544 42.28787)
        1 | 2012-01-27 14:36:53.608661+00 | POINT(-71.063544 42.28987)
        1 | 2012-01-27 14:37:10.841856+00 | POINT(-71.066544 42.3)
        2 | 2012-01-27 14:37:46.400954+00 | POINT(-0.8 53.2)
        2 | 2012-01-27 14:37:55.752122+00 | POINT(-0.81 53.21)
        2 | 2012-01-27 14:38:07.773077+00 | POINT(-0.82 53.22)

The next step is to aggregate those points into linestrings. This is done with the ST_MakeLine() function:

SELECT gps.track_id, ST_AsText(ST_MakeLine(gps.the_geom)) AS track_line
    FROM (SELECT track_id, gpstime, the_geom FROM points ORDER BY track_id, gpstime) as gps
    GROUP BY gps.track_id
    ORDER BY gps.track_id;

Which yields:

 track_id |                             track_line
----------+---------------------------------------------------------------------
        1 | LINESTRING(-71.064544 42.28787,-71.063544 42.28987,-71.066544 42.3)
        2 | LINESTRING(-0.8 53.2,-0.81 53.21,-0.82 53.22)

Finally you need to get the length of that line over the spheroid of your choice. If you're using GPS then inevitably it'll be the WGS84 spheroid. The function to use is ST_Length_Spheroid().

So putting that all together, my final query ends up looking like this:

SELECT gps.track_id, tracks.name, ST_Length_Spheroid(ST_MakeLine(gps.the_geom), 'SPHEROID["WGS 84",6378137,298.257223563]') AS track_len
    FROM (SELECT track_id, gpstime, the_geom FROM points ORDER BY track_id, gpstime) as gps, tracks
    WHERE gps.track_id = tracks.track_id
    GROUP BY gps.track_id, tracks.name
    ORDER BY gps.track_id;

And I get these results:

 track_id | name  |    track_len
----------+-------+------------------
        1 | Larry | 1389.08015675763
        2 | Curly | 2596.09131911171

If you want to limit it to a specific track, change your WHERE clause to something like this:

 WHERE gps.track_id = tracks.track_id AND tracks.name = 'Larry'

If you're using geography columns instead of geometry, you can eliminate the spheroid part and just use ST_Length(). Note that the second parameter of ST_Length() when using geography columns should be set to TRUE if you want more accurate but slower calculations. The only drawback is that geography columns are limited to the WGS84 spheroid.

[GIS] How to sort points from a SpatialPoints object over a SpatialLines object according to line direction using R

I could answer my question traveling through every line segment in the SpatialLines object and finding the intersected points from the SpatialPoints object. I made a function that is all explained by comments.

# Load packages
library('rgeos')
library('data.table')

# Vector norm
norm_vec <- function(x) sqrt(sum(x^2))    

# OrderPointsAlongLine function
OrderPointsAlongLine <- function(spatialLinesObject, spatialPointsObject) {

  lineCoords = as.data.frame(coordinates(spatialLinesObject))
  lineCoordsLength = nrow(lineCoords)

  # Make a buffer 
  spatialPointsObjectBuffer = gBuffer(spgeom = spatialPointsObject, width = 1, byid = TRUE)

  # Aux list and data frame
  cellsAndDistances = list()
  cd = data.frame('cell' = NA, 'dist' = NA)

  # Go sequentially through every line segment in the SpatialLine object looking for which SpatialPoints are intersected
  for(i in 1:(lineCoordsLength - 1)) {

    # Unommented: show iteration number
    # cat(paste("\n iteration", i, sep = ""))

    # Line segment points
    p0 = lineCoords[i,]
    p1 = lineCoords[i+1,]

    # Build line segment
    SL = SpatialLines(LinesList = list(Lines(Line(coords = rbind(p0,p1)), ID = paste(i, sep = ""))), proj4string = CRS(epsg.32721))

    # Find intersected points with line segment
    intersect = which(gIntersects(spgeom1 = spatialPointsObjectBuffer, spgeom2 = SL, byid = TRUE))

    # Points intersected
    p2 = spatialPointsObject[intersect,]

    p2Length = nrow(p2@coords)

    # If there is almost one point intersected
    if(p2Length != 0) {

      # If there are more than one intersected points
      for(j in 1:p2Length) {

        # Commented: show sub-iteration number
        # cat(paste("\n sub-iteration", j, sep = ""))

        cd[j,] = data.frame('cell' = intersect[j], 'dist' = norm_vec(coordinates(p2)[j,] - p0))

      }
    }

    # Order by distance
    data.table::setorder(cd, dist) 
    cellsAndDistances[[i]] = cd

  }

  # Bind list into one data farme
  cellsAndDistancesBind = do.call("rbind", cellsAndDistances)

  # Remove duplicated
  cellsAndDistancesBind <- cellsAndDistancesBind[-which(duplicated(cellsAndDistancesBind$cell)),]

  # Return
  return(cellsAndDistancesBind)

}

# Call function
orderedPoints <- OrderPointsAlongLine(spatialLinesObject = spatialLinesObject, spatialPointsObject = spatialPointsObject)

# Plots
plot(spatialLinesObject, xlab = "Longitude", ylab = "Latitude", main = "", col = "green")
plot(spatialPointsObject[orderedPoints$cell,], add = TRUE, pch = 19, cex = 0.5)
text(coordinates(spatialPointsObject[orderedPoints$cell,]), labels = seq_along(spatialPointsObject[orderedPoints$cell,]), cex = 0.75, pos = 1)

Improvements are welcome!

Best Answer

Related Solutions

PostGIS Distance – Calculating Total Great Circle Distance Along a Path of Points

[GIS] How to sort points from a SpatialPoints object over a SpatialLines object according to line direction using R

Related Question