[GIS] How to snap one set of points to another in R

rsp

I would like to snap one set of points (SpatialPoints* objects) to another in R using Euclidean distance. I'm hoping for a function like maptools::snapPointsToLines.

The attributes of the points don't need to be transferred.

Any ideas?

Best Answer

Here's a reproducible answer and a function that I think solves the problem. It all relies on nncross from the spatstat package.

Step 1: Load the packages we'll be using

library(sp)
library(spatstat)
library(maptools) # to convert to ppp

Step 2: Create two small sets of points, give one attribute data:

set.seed(2014) # ensure reproducibility 
x <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2))
y <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2)) 
# add some attribute data
x <- SpatialPointsDataFrame(x, data = data.frame(value = 1:length(x)))

Step 3: Create a function to allocate coordinates of x those of the nearest y:

gSnap <- function(x, y){
  x_ppp <- as.ppp(x)
  y_ppp <- as.ppp(y)
  nearest_y_from_x <- nncross(x_ppp, y_ppp)
  x_new_coords <- y_ppp[ nearest_y_from_x$which, ]
  x_new <- as.SpatialPoints.ppp(x_new_coords)
  x_new <- SpatialPointsDataFrame(x_new, x@data)
  x_new
}

Step 4: Test the output:

xSnapped <- gSnap(x, y)

plot(spRbind(x, y), col = "white")
points(x, col = "red")
points(y, col = "green")
points(xSnapped, pch = 3)

The output is shown below - all the red 'x' points have snapped to only 3 green 'y' points and their attribute data is maintained. Please test on your (perhaps larger and more complex) dataset and let me know if it works.

There may be a pre-existing function to do this and almost certainly a more efficient implementation.

gSnap-in-action

Related Question