-
-
Save barryrowlingson/89914f5ed94894f1ba30 to your computer and use it in GitHub Desktop.
computing nearest neighbours with RANN
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## lets use this nearest neighbour package: | |
require(RANN) | |
## create a village as a cluster of houses around x,y: | |
village <- function(x, y, n, id, sd){ | |
## return a data frame of x, y, and id values: | |
data.frame(x=x+rnorm(n,0,sd), | |
y=y+rnorm(n,0,sd), | |
id=id | |
) | |
} | |
## create some villages with some number of houses and spatial sd | |
## with centres in unit square: | |
villages <- function(nvillages, mu, sd){ | |
## return a data frame of x, y, and village id values: | |
do.call(rbind, | |
lapply(1:nvillages,function(i){ | |
village(runif(1),runif(1),rpois(1,mu)+1,i,sd) | |
}) | |
) | |
} | |
## create a grid over the x and y coordinates: | |
vgrid <- function(villages, ng=100){ | |
grid = expand.grid( | |
x = seq(min(villages$x), max(villages$x),len=ng), | |
y = seq(min(villages$y), max(villages$y),len=ng) | |
) | |
grid | |
} | |
## make some data: | |
set.seed(3) | |
## 3 villages, with about 5 houses, spread about 0.06 | |
v = villages(3, 5, .06) | |
vg = vgrid(v, 100) | |
## now compute first (k=1) nearest neighbour to each grid point | |
## limit search to 0.1 units: | |
nns = nn2(v[,c("x","y")], vg, k=1, radius=.1, searchtype="radius") | |
## index of nearest village locations is this component. | |
## when a grid point has no neighbour in range this is zero: | |
iv = nns$nn.idx | |
## we are going to use this to index into the village data frame, | |
## and indexing by zero is going to fail (comment this line out to | |
## see how badly), so set to NA | |
iv[iv==0] <- NA | |
## plot the grid coloured by colour of corresponding village id: | |
plot(vg$x, vg$y, col=v$id[iv], pch=19, asp=1) | |
## add village points with a white outline marker in the village colour. | |
## you should just see white circles with the inner colour matching the | |
## region colour: | |
points(v$x, v$y, bg=v$id, asp=1, pch=21, col="white") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment