Snapping points to lines is a common task when dealing with points that have to be located over lines for some specific reason.

One of the most common applications is to **locate trajectory data on a road network**. As such, it is a georeferencing problem that aims to get enhanced locations of trajectories getting rid of their deviations by using road networks as reference.

This problem, also called **point-to-curve matching**, is one of the simplest approaches for **map-matching** algorithms, which can also include topology, probability and fuzzy logic, among others. A complete review of map-matching algorithms is given by (Quddus et al., 2007).

In terms of free and open source software for Geomatics (**FOSS4G**), libraries such as Java Topology Suite (JTS) and its port to C called GEOS provide classes and functions that can be used to snap points to lines. Some examples on how to use them are provided by P. Ramsey (using SQL and PostGIS functions), Marcello Benigno (using GRASS) and GeoTools User Guide.

However, **the R package rgeos**, which is the R interface to GEOS, still doesn’t provide those functions and one has to build his own ones. In the rest of this post I will show you my solution, which **snaps the trajectory points to the nearest points lying on the road network**.

First, import some road data into R (you can download the shapefile from this link, source: OpenStreetMap), don’t forget to adjust the path to the data:

1 2 3 4 5 |
require("rgdal") crs = CRS("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs") l = readOGR(dsn="/docs/geodata/munster/osm/", layer="subset_street_graph") proj4string(l) = crs |

Import trajectory data from http://ifgi.uni-muenster.de/~epebe_01/pm10/All4_1a.txt and ensure that both (road network and trajectory) have the same projected reference system (for distance calculation):

6 7 8 9 10 11 12 |
x = read.csv("http://ifgi.uni-muenster.de/~epebe_01/pm10/All4_1a.txt") x$t = as.POSIXct(strptime(do.call(paste, x[c("Date", "Time")]), "%Y/%m/%d %H:%M:%S")) m = matrix(c(x[,2],x[,1]),ncol=2) df = data.frame(x$t, x[,5],x[,6]) p = SpatialPointsDataFrame(coords=m,data=df,proj4string=CRS("+proj=longlat +datum=WGS84")) p = spTransform(p,crs) |

Define the functions **nearestPointOnSegment**, **nearestPointOnLine**, this way:

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
nearestPointOnSegment = function(s, p){ # Adapted from http://pastebin.com/n9rUuGRh ap = c(p[1] - s[1,1], p[2] - s[1,2]) ab = c(s[2,1] - s[1,1], s[2,2] - s[1,2]) t = sum(ap*ab) / sum(ab*ab) t = ifelse(t<0,0,ifelse(t>1,1,t)) x = s[1,1] + ab[1] * t y = s[1,2] + ab[2] * t c(x, y, (x-p[1])^2 + (y-p[2])^2) # Return nearest point and distance } nearestPointOnLine = function(coordsLine, coordsPoints){ nearest_points = vapply(2:nrow(coordsLine), function(x) nearestPointOnSegment(coordsLine[(x-1):x,], coordsPoints), FUN.VALUE=c(0,0,0)) # Return coordinates of the nearest point in this line nearest_points[1:2, which.min(nearest_points[3,])] } |

Finally, define the main function called **snapPointsToLines**:

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 |
snapPointsToLines <- function( points, lines, maxDist=NA, withAttrs=TRUE) { require("rgeos") if (!is.na(maxDist)) { w = gWithinDistance(points, lines, dist=maxDist, byid=TRUE) validPoints = apply(w,2,any) validLines = apply(w,1,any) points = points[validPoints,] lines = lines[validLines,] } d = gDistance(points, lines, byid=TRUE) nearest_line_index = apply(d, 2, which.min) # Position of each nearest line in lines object coordsLines = coordinates(lines) coordsPoints = coordinates(points) # Get coordinates of nearest points lying on nearest lines mNewCoords = vapply(1:length(points), function(x) nearestPointOnLine(coordsLines[[nearest_line_index[x]]][[1]], coordsPoints[x,]), FUN.VALUE=c(0,0)) # Recover lines' Ids (Ids and index differ if maxDist is given) if (!is.na(maxDist)) nearest_line_id = as.numeric(rownames(d)[nearest_line_index])+1 else nearest_line_id = nearest_line_index # Create data frame and sp points if (withAttrs) df = cbind(points@data, nearest_line_id) else df = data.frame(nearest_line_id, row.names=names(nearest_line_index)) SpatialPointsDataFrame(coords=t(mNewCoords), data=df, proj4string=CRS(proj4string(points))) } |

Now you can use the snapPointsToLines function with the whole datasets:

68 |
snappedPoints = snapPointsToLines(p, l) |

Or you can establish a maximum distance (in this case in meters) to avoid snapping points that are farther apart:

69 |
snappedPoints = snapPointsToLines(p, l, 100) |

You can also get rid of the original attributes and just get the coordinates of the snapped points and the id of their nearest line:

70 |
snappedPoints = snapPointsToLines(p, l, 100, FALSE) |

This is what I get executing the line 69 and calling the function plot (Red crosses are original trajectory points and green ones are snapped points):

71 72 73 |
plot(p,col="red") plot(snappedPoints,add=TRUE,col="green") plot(l,add=TRUE) |

A more detailed view (Red crosses are original trajectory points and blue ones are snapped points):

A basic way to improve the result is to pass a filtered road network, e.g., only main roads. Of course, it depends on the trajectory. This is the result of snapping the points to a filtered road network (using 100 meters as maximum distance):

As I’m still unfamiliar with R data structures, this solution could be improved in terms of performance. With these test datasets this code is almost twice slower than a PostGIS implementation.

**REFERENCES**

- M. Quddus, W. Ochieng and R. Noland. Current map-matching algorithms for transport applications : State-of-the art and future research directions. In Elsevier, Kidlington, ROYAUME-UNI, editor, Transportation research. Part C, Emerging technologies, vol. 15, no5, pages 312-328, 2007.

**UPDATE (2012.06.30)**

I’m pleased to announce the functions *snapPointsToLines*, *nearestPointOnSegment* and *nearestPointOnLine* have been included into the R package called **maptools** (v. 0.8-16).

The functions can be called directly after importing **maptools** into the R session, this way:

# install.packages("maptools") # For installing the package require("maptools") # Create both lines and points (see lines 1-5 and 6-12 in the post) # Finally, call the main function snappedPoints = snapPointsToLines(p, l, 100)

## Recent Comments