R/snapPointsToLines_v2.R

In KevinSee/QRFcapacity: Carrying Capacity with QRF

#' @title Snap Points to Lines
#'
#' @description This function snaps a set of points to a set of lines based on the minimum distance of each point to any of the lines. This function does not work with geographic coordinates. This function replaces the maptools::snapPointsToLines because it stopped working correctly
#'
#' @author Kevin See
#'
#' @param points An object of the class SpatialPoints or SpatialPointsDataFrame.
#' @param lines An object of the class SpatialLines or SpatialLinesDataFrame.
#' @param maxDist Numeric value for establishing a maximum distance to avoid snapping points that are farther apart; its default value is NA.
#' @param withAttrs Boolean value for preserving (TRUE) or getting rid (FALSE) of the original point attributes. Default: TRUE. This parameter is optional.
#' @param idField A string specifying the field which contains each line's id. This id will be transferred to the snapped points data set to distinguish the line which each point was snapped to.
#'
#' @import rgeos sp maptools
#' @return SpatialPointsDataFrame object as defined by the R package 'sp'. This object contains the snapped points, therefore all of them lie on the lines.
#' @export

snapPointsToLines_v2 <- function( points, 
                                  lines, 
                                  maxDist = 0.1, 
                                  withAttrs = TRUE, 
                                  idField = NA) {
  
  if (rgeosStatus()) {
    if (!requireNamespace("rgeos", quietly = TRUE))
      stop("package rgeos required for snapPointsToLines")
  } else
    stop("rgeos not installed")
  
  if (class(points) == "SpatialPoints" && missing(withAttrs))
    withAttrs = FALSE
  
  if (class(points) == "SpatialPoints" && withAttrs==TRUE)
    stop("A SpatialPoints object has no attributes! Please set withAttrs as FALSE.")
  
  d = rgeos::gDistance(points, lines, byid=TRUE) 
  
  if(!is.na(maxDist)){
    distToLine <- apply(d, 2, min, na.rm = TRUE)  
    validPoints <- distToLine <= maxDist  # indicates which points are within maxDist of a line
    distToPoint <- apply(d, 1, min, na.rm = TRUE)
    validLines <- distToPoint <= maxDist 
    
    # Drop elements beyond maxdist
    points <- points[validPoints, ]
    lines = lines[validLines, ]
    d  = d[ validLines,  validPoints, drop = FALSE]
    distToLine <- distToLine[validPoints]
    
    # If no points are within maxDist return an empty SpatialPointsDataFrame object 
    if(!any(validPoints)){
      if(is.na(idField)){
        idCol = character(0)
      } else {
        idCol = lines@data[,idField][0]
      }
      newCols = data.frame(nearest_line_id  = idCol, snap_dist = numeric(0))
      if(withAttrs) df <- cbind(points@data, newCols) else df <- newCols
      res <- SpatialPointsDataFrame(points, data=df, 
                                    proj4string=CRS(proj4string(points)), match.ID = FALSE)
      return(res)
    }
    
  } else { # If no maxDist arg still calculate distToLine so it can be returned
    distToLine = apply(d, 2, min, na.rm = 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 (If no id field has been specified, take the sp-lines id)
  if (!is.na(idField)) { 
    nearest_line_id = lines@data[,idField][nearest_line_index] 
  }  else {
    nearest_line_id = sapply(slot(lines, "lines"), function(i) slot(i, "ID"))[nearest_line_index] 
  }
  # Create data frame and sp points
  if (withAttrs) df = cbind(points@data, 
                            data.frame(nearest_line_id, snap_dist = distToLine)) 
  else df = data.frame(nearest_line_id, 
                       snap_dist = distToLine, 
                       row.names=names(nearest_line_index))
  
  SpatialPointsDataFrame(coords=t(mNewCoords), 
                         data=df, 
                         proj4string=CRS(proj4string(points)))
}