# Calculation of average point distances from one another within a polygon grid

I am working with processed LiDAR data of forested areas in ArcMap 10, with various area based statistics extracted for each 30x30m polygon grid-cell in a shapefile. There are approximately 20,000 of these cells. I also have a layer of points relating to individual tree locations for the whole site.

What I need to do is calculate the average distance between the tree-points and their standard deviations within each grid-cell polygon. And then write those two values to the corresponding grid-cell.

I have used the GME toolbox (http://www.spatialecology.com/gme/) to derive this data for areas corresponding to my field site locations; however this is simply not feasible over such a large area.

Does anyone have any ideas on how I can solve this issue?

-

Here is an R solution, intended to function as pseudocode for implementation on any appropriate platform (C++, Python, etc) and to be a working prototype. It begins with a function to compute the mean and SD distances of an array of points:

#
# Compute distance statistics for points (x,y).
#
dist.stats <- function(x,y) {
#
# Compute all distances between distinct points.
#
xy <- cbind(x,y)
f <- function(w) apply(xy, 1, function(v) sum((v - w)^2))
distances <- apply(xy, 1, f)
distances <- sqrt(distances[lower.tri(distances)])
#
# Return the distance statistics.
#
c(mean(distances), sd(distances), length(x))
}

This function needs to be applied cell by cell. Here, I presume the cells are regularly arranged in a grid parallel to the coordinate axes. This enables the points to be grouped by means of arithmetic operations. (If they have already been grouped by polygon (by virtue of a spatial join), the code would be simpler: two lines would split the x and y coordinates by polygon id and a third line would apply the block statistics to each group.)

blockstats <- function(f,x,y,cellsize=1,n.cols=1,n.rows=1,origin.x=0,origin.y=0) {
#
# Split points by column.
#
cols <- factor(floor((x-origin.x)/cellsize) + 1, 1:n.cols)
x.cols <- split(x, cols)
y.cols <- split(y, cols)
rows.cols <- split(rows, cols)
#
# Split columns by rows.
#
g <- function(z,g) split(z, factor(g, 1:n.rows))
x.groups <- mapply(g, x.cols, rows.cols)
y.groups <- mapply(g, y.cols, rows.cols)
#
# Apply summary function `f` to each group.
#
mapply(dist.stats, x.groups, y.groups)
}

To illustrate their use, let's create a small sample dataset:

cellsize <- 30
n.rows <- 10
n.cols <- 20
n <- n.rows * n.cols * 5
set.seed(17)
x <- rnorm(n, mean=1/2, sd=1/4) * cellsize * n.cols
y <- rnorm(n, mean=1/2, sd=1/6) * cellsize * n.rows

With these points in hand, the block statistics are computed, converted to raster format, and plotted:

system.time(results <- blockstats(dist.stats, x, y, cellsize, n.cols, n.rows))
r.mean <- matrix(results[1,], nrow=n.rows, ncol=n.cols)
r.sd <- matrix(results[2,], nrow=n.rows, ncol=n.cols)
r.n <- matrix(results[3,], nrow=n.rows, ncol=n.cols)

require("raster")
raster.plot <- function(r,s) {
plot(raster(r, xmn=0, xmx=n.cols*cellsize, ymn=0, ymx=n.rows*cellsize), main=s)
points(x, y, cex = min(1, 10/sqrt(n)))
}
par(mfrow=c(2,2))
raster.plot(r.mean, "Mean block distance")
raster.plot(r.sd, "SD block distance")
raster.plot(r.n, "# points")

Increasing n.rows from 10 to 100 and n.cols from 20 to 200 simulates the situation in the question: about 100,000 points covering 20,000 cells on a 30 m cellsize grid. The timing is 20 seconds. Dedicated compiled code ought to go several orders of magnitude faster, but even this slow speed may be adequate for the problem.

-
Thank you for the detailed answer! You are correct the polygon cells are regularly arraigned in a grid parallel to the British National grid coordinate system. Each tree point has 3 attributes: the XY coordinates and the intersecting polygon cell ID. You’ll have to forgive me, but I am still in the process of learning how to use R and as such some of what you wrote I did not understand. How might the code be modified to group the tree points by polygon ID? – MattS May 9 '12 at 15:34
Change blockstats to include an array of polygon ids among its arguments and split the coordinates on that array, rather than computing row and column indexes and splitting on those. Finish up with mapply exactly as before. – whuber May 9 '12 at 21:02
Thanks again! I've run the script and got the results I needed. – MattS May 14 '12 at 10:22