1nearestPointOnSegment = function(s, p){
2    # Adapted from http://pastebin.com/n9rUuGRh
3    ap = c(p[1] - s[1,1], p[2] - s[1,2])
4    ab = c(s[2,1] - s[1,1], s[2,2] - s[1,2])
5    t = sum(ap*ab) / sum(ab*ab)
6    t = ifelse(t<0,0,ifelse(t>1,1,t))
7    t = ifelse(is.na(t), 0, t) # when start and end of segment are identical t is NA
8    x = s[1,1] + ab[1] * t
9    y = s[1,2] + ab[2] * t
10    result = c(x, y, sqrt((x-p[1])^2 + (y-p[2])^2))  # Return nearest point and distance
11    names(result) = c("X","Y","distance")
12    result
13}
14
15nearestPointOnLine = function(coordsLine, coordsPoint){
16    nearest_points = vapply(2:nrow(coordsLine),
17        function(x)
18            nearestPointOnSegment(coordsLine[(x-1):x,], coordsPoint),
19                FUN.VALUE=c(0,0,0))
20
21    # Return coordinates of the nearest point on this line
22    nearest_points[1:2, which.min(nearest_points[3,])]
23}
24
25snapPointsToLines <- function( points, lines, maxDist=NA, withAttrs=TRUE, idField=NA) {
26
27    if (rgeosStatus()) {
28    	if (!requireNamespace("rgeos", quietly = TRUE))
29			stop("package rgeos required for snapPointsToLines")
30    } else
31        stop("rgeos not installed")
32
33    if (is(points, "SpatialPoints") && missing(withAttrs))
34        withAttrs = FALSE
35
36    if (is(points, "SpatialPoints") && withAttrs==TRUE)
37        stop("A SpatialPoints object has no attributes! Please set withAttrs as FALSE.")
38
39    d = rgeos::gDistance(points, lines, byid=TRUE)
40
41    if(!is.na(maxDist)){
42      distToLine <- apply(d, 2, min, na.rm = TRUE)
43      validPoints <- distToLine <= maxDist  # indicates which points are within maxDist of a line
44      distToPoint <- apply(d, 1, min, na.rm = TRUE)
45      validLines <- distToPoint <= maxDist
46
47      # Drop elements beyond maxdist
48      points <- points[validPoints, ]
49      lines = lines[validLines, ]
50      d  = d[ validLines,  validPoints, drop = FALSE]
51      distToLine <- distToLine[validPoints]
52
53      # If no points are within maxDist return an empty SpatialPointsDataFrame object
54      if(!any(validPoints)){
55        if(is.na(idField)){
56          idCol = character(0)
57        } else {
58          idCol = lines@data[,idField][0]
59        }
60        newCols = data.frame(nearest_line_id  = idCol, snap_dist = numeric(0))
61        if(withAttrs) df <- cbind(points@data, newCols) else df <- newCols
62        res <- SpatialPointsDataFrame(points, data=df,
63                               proj4string=CRS(proj4string(points)), match.ID = FALSE)
64        return(res)
65      }
66
67    } else { # If no maxDist arg still calculate distToLine so it can be returned
68      distToLine = apply(d, 2, min, na.rm = TRUE)
69    }
70
71    nearest_line_index = apply(d, 2, which.min) # Position of each nearest line in lines object
72
73    coordsLines = coordinates(lines)
74    coordsPoints = coordinates(points)
75
76    # Get coordinates of nearest points lying on nearest lines
77    mNewCoords = vapply(1:length(points),
78        function(x)
79            nearestPointOnLine(coordsLines[[nearest_line_index[x]]][[1]],
80                coordsPoints[x,]), FUN.VALUE=c(0,0))
81
82    # Recover lines' Ids (If no id field has been specified, take the sp-lines id)
83    if (!is.na(idField)) {
84      nearest_line_id = lines@data[,idField][nearest_line_index]
85    }  else {
86      nearest_line_id = sapply(slot(lines, "lines"), function(i) slot(i, "ID"))[nearest_line_index]
87    }
88    # Create data frame and sp points
89    if (withAttrs) df = cbind(points@data, data.frame(nearest_line_id, snap_dist = distToLine))
90    else df = data.frame(nearest_line_id, snap_dist = distToLine, row.names=names(nearest_line_index))
91
92    SpatialPointsDataFrame(coords=t(mNewCoords), data=df,
93        proj4string=CRS(proj4string(points)))
94}
95