3

I recently find a problem of snapPointsToLines. It can't keep the attributes of the spatial point dataframe. The example is as below:

# Generate a spatial line dataframe
l1 = cbind(c(1,2,3),c(3,2,2))
l1a = cbind(l1[,1]+.05,l1[,2]+.05)
l2 = cbind(c(1,2,3),c(1,1.5,1))
Sl1 = Line(l1)
Sl1a = Line(l1a)
Sl2 = Line(l2)
S1 = Lines(list(Sl1, Sl1a), ID="a")
S2 = Lines(list(Sl2), ID="b")
Sl = SpatialLines(list(S1,S2))
df = data.frame(z = c(1,2), row.names=sapply(slot(Sl, "lines"), function(x) slot(x, "ID")))
Sldf = SpatialLinesDataFrame(Sl, data = df)

# Generate a spatial point dataframe
xc = c(1.2,1.5,2.5)
yc = c(1.5,2.2,1.6)
Spoints = SpatialPoints(cbind(xc, yc))
Spdf <- SpatialPointsDataFrame(Spoints, data = data.frame(value = 1:length(Spoints)))

#use the function SpatialPointsDataFrame 
res <- snapPointsToLines(Spdf, Sldf)

res only has "nearest_line_id" and "snap_dist". It doesn't have "value" field from Spdf, which I need. 

#use the function SpatialPointsDataFrame with "withAttrs = TRUE" parameter
res <- snapPointsToLines(Spdf, Sldf, withAttrs = TRUE)

It reports error: 

"Error in snapPointsToLines(Spdf, Sldf, withAttrs = TRUE) : 
  A SpatialPoints object has no attributes! Please set withAttrs as FALSE."

But Spdf is the spatialpointdataframe with attribute. I don't know what problem it is. When I used this function several weeks ago, it didn't have this problem.

Freeego
  • 135
  • 6

1 Answers1

2

I think the problem may be due to the function itself. When you look at the codes of this function, we can see the codes at the beginning part as below.

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.")

Sometimes a SpatialPointsDataFrame could be identified as SpatialPoints. So the function will treat your SpatialPointsDataFrame as SpatialPoints and will not keep the attributes in the function. You can make a little modification in the the codes of the function as below.

snapPointsToLines1 <-  function (points, lines, maxDist = NA, withAttrs = TRUE, idField = NA) 
{
  if (rgeosStatus()) {
    if (!requireNamespace("rgeos", quietly = TRUE)) 
      stop("package rgeos required for snapPointsToLines")
  }
  else stop("rgeos not installed")
  if (is(points, "SpatialPointsDataFrame")==FALSE && missing(withAttrs)) 
    withAttrs = FALSE
  if (is(points, "SpatialPointsDataFrame")==FALSE && withAttrs == TRUE) 
    stop("A SpatialPointsDataFrame object is needed! 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
    distToPoint <- apply(d, 1, min, na.rm = TRUE)
    validLines <- distToPoint <= maxDist
    points <- points[validPoints, ]
    lines = lines[validLines, ]
    d = d[validLines, validPoints, drop = FALSE]
    distToLine <- distToLine[validPoints]
    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 {
    distToLine = apply(d, 2, min, na.rm = TRUE)
  }
  nearest_line_index = apply(d, 2, which.min)
  coordsLines = coordinates(lines)
  coordsPoints = coordinates(points)
  mNewCoords = vapply(1:length(points), function(x) nearestPointOnLine(coordsLines[[nearest_line_index[x]]][[1]], 
                                                                       coordsPoints[x, ]), FUN.VALUE = c(0, 0))
  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]
  }
  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)))
}

Then using this new function snapPointsToLines1, you can get the attributes what you want.

Alvin_z
  • 36
  • 2