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