Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm sure this is quite basic, but I just couldn't find anything appropriate: I am looking for a function that applies another function to variables of a SpatialPointsDataFrame, but only for locations that are the nearest neighbors of a specified location.

Alternatively, how can I subset a SpatialPointsDataFrame by distance to a given location?

For example, in the meuse data set (gstat package) I would like to calculate mean, variance and quantiles of zinc concentration for the 20 closest locations that are at most 0.1km away from location x=179997, y=331662.

I found the function applynbd() in the spatstat package, but it applies the given function to all points in the data frame insteat of a single one. Moreover, only a number of nearest neighbors OR a distance (max. radius) can be specified instead of both.

share|improve this question
1  
why don't you try to do it yourself? the class SpatialPointsDataFrame makes it easy to extract the coordinates (coordinates(meuse)). after you just need to calculate the distance from each point to the (xref,yref) and select the ones that are at most 0.1km away. –  A.R Nov 29 '12 at 15:56

1 Answer 1

up vote 2 down vote accepted

The following code will get you a list of distances from a point to a set of points by using the spDistsN1 function from the sp package. Let's assume that the point you are looking for is not part of the meuse SpatialPointsDataFrame:

require(sp)
data(meuse)
coordinates(meuse) = c("x","y")
dist_vector = spDistsN1(meuse, c(179997, 331662))

...or if it is one of the existing points in the SPDF, let's say the tenth point:

dist_vector = spDistsN1(meuse, meuse[10,])
# or even better, add it as a column to meuse
meuse[["distance"]] = spDistsN1(meuse, meuse[idx,])

and then you can select the 20 closest points:

meuse_dat = as.data.frame(meuse)
meuse_dat = meuse_dat[order(meuse_dat[["distance"]]),]
meuse_dat[["close"]] = FALSE
meuse_dat[["close]][1:20] = TRUE
meuse_subset = meuse_dat[close && distance < 1000,]

and calculate some stats:

require(plyr)
summarise(meuse_subset, mn = mean(zinc), vr = var(zinc))
share|improve this answer
    
Thank you, @Paul Hiemstra, I think that'll get me started. I just thought this task was so common that there had to be a dedicated function in some package that is optimized in terms of performance. I could imagine your apprach may not be efficient when applied to many, say, tens of thousands of points via a for loop or similar? Any suggestions on efficiency? Thanks again! –  Nima Dec 3 '12 at 23:23
    
If the points you want to know the distance of are always part of the original dataset, you can precompute the distance matrix once to make it more efficient. For any subsequent steps you then do not need to calculate anything, but just refer to the distance matrix. An easy way to get the distance matrix is to use spDists. –  Paul Hiemstra Dec 4 '12 at 8:06
    
To get raw performance, you could use Rcpp to write C++ code to calculate the distances, this should be very fast, at least orders faster than a pure R solution. –  Paul Hiemstra Dec 4 '12 at 8:07

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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