# Stats on individual grid points nearest neighbours

I have shp of points where each point is the centroid of the individual cells of a 1*1km grid (~ 500 000 points). For each of these points, I have their coordinates and elevation.

What I want: for each of these points, I want to compute the mean, std dev, min, max and median elevation of the nearest 8 surrounding points.

What I did so far: I wrote an R script to identify each of the nearest points for each of the points of the grid. But considering the size of the dbf and the general slowness of R, it takes forever.

My tools: QGIS, ArcGIS, SAS, R or other free software if needed.

Does any one have a suggestion to do this (1- identify the direct neighbors, 2- compute summary stats from the elevation of these neighbors) in a more efficient way, using GIS?

-
If you have access to any form of raster GIS, then please go back to the original raster DEM and compute focal (aka neighborhood) statistics: the calculations will be instantaneous. Also, although R is slow, using the raster library it will still complete this task faster than anything using your point shapefile as input. – whuber Nov 20 '12 at 19:05
Thanks a lot Whuber for these suggestions. In a first glance, it seams quite straight forward to do this using ArcGIS Focal Statistics tool or R raster package. I wander if there is a way to do it in QGIS (the neighour analysis tool in QGIS is of no help here)? ... – Helene Genet Nov 20 '12 at 19:20
See here for QGIS nearest neighbor. The rest of the analysis (mean, std dev, etc.) should be done elsewhere. – Zach May 23 '13 at 14:31

There is the spdep R package, which has a K-nearest neighbor algorithm implemented. I am not sure if its more efficient than what you have, but it might be worth a try.

Here is an example script to show how you could calculate the values you wanted and put them into a data frame.

library(spdep)

example(columbus)

# Extract coordinates from shape
coords <- coordinates(columbus)

# Use the area of the example data as a replacement for elevation
# (select column with elevation value from attribute table here)
area.vector <- columbus\$AREA

# Calculate K nearest neighbors using spdep package
col.knn <- knearneigh(coords, k=8)

# Extract neighbor index matrix
neighbors <- col.knn\$nn

# Prepare empty dataframe for calculations
output<-data.frame(area.vector, mean=NA, sd=NA, min=NA, max=NA, median=NA)

# Loop over neighbors matrix, calculate values
for(i in 1:nrow(columbus)){
values <- area.vector[neighbors[i,]]
output[i,2:6] <- c(mean(values), sd(values), min(values), max(values), median(values))
}

# View head of output to see if it worked