5

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?

2

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

  • You do not need maptools to coerce to a ppp object. You can just pass the coordinates to the ppp function: x.ppp <- ppp(coordinates(x), win) although, ppp does expect a window so you will have to create one using something like "convexhull.xy", which also accepts the matrix in the coordinates slot. – Jeffrey Evans Nov 11 '14 at 19:07
2

Starting from Robin's solution, here is an alternative that uses less packages and spDists's ability to compute pairwise distances between two point sets:

library(sp)
set.seed(2014)
x <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2))
y <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2)) 
plot(x, col = "red")
points(y, col = "green")
snap = apply(spDists(x, y), 1, which.min)
points(y[snap,], pch = 3)
  • glad to have you on board! Your expertise in R spatial classes will be greatly appreciated. I had not thought of wrapping spDists in apply. Very efficient. – Jeffrey Evans Feb 24 '15 at 0:12
0

You can use spDistsN1 and which.min to assign the coordinates of the nearest points in another spatial points object. Things will remain ordered so, you could just assign the data from the original data back to the adjusted points.

Add package and create offset points

require(sp)
  data(meuse)
    pts <- meuse[1:10,]
      pts2 <- data.frame( x=jitter(pts[,1],factor=10), y=jitter(pts[,2],factor=10), 
                          pts[,3:ncol(pts)])  
                            coordinates(pts) <- ~x+y  
                              coordinates(pts2) <- ~x+y  

plot(pts, pch=20, col="black")
  plot(pts2, pch=20, col="red", cex=0.75, add=TRUE) 

Loop to assign coordinates, based on smallest distance, from another sp point object

new.coords <- matrix(ncol=2)
  for(i in 1:nrow(pts2)) {
    d <- spDistsN1(pts, pts2[i,])
      new.coords <- rbind(new.coords,coordinates(pts)[which.min(d),])
    }  
adj.pts <- SpatialPoints(new.coords[-1,]    )

plot(pts, pch=20, col="green")
  plot(pts2, pch=20, col="red", cex=0.75, add=TRUE)
    plot(adj.pts, pch=20, col="black", add=TRUE) 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.