[GIS] Calculating average raster to point distances within census tracts

arcgis-10.0arcgis-desktopdistancerastervector

I have some experience using ArcGIS, but mostly with the field calculator in the attribute table/spatial joins, so I apologize if I say something that doesn't make sense.

To provide a bit of background, I have a vector layer containing census tract data. I have another layer with geocoded supermarket locations. I'm trying to determine the number of supermarkets per census tract, which I can then use for various regression models. When I first tackled this issue I simply used raw counts for each census tract. Those counts, however, assume that residents from one census tract have absolutely no access to supermarkets in other census tracts. I've been struggling with a way to account for this in my regression model, so I decided to ditch the regression statistics for now and find a way to account for this issue in my count data.

After doing some searching, I found a method described by Liam Downey (2003) that would seem appropriate for what I'm trying to do. I'm not sure if we're allowed to directly quote academic articles, so I'll paraphrase his method in steps.

Convert census tract maps to raster maps (he used 25x25m grid cells).
Determine the number of facilities within a 1 km radius of the center of each grid cell.
Calculate the distance (m) from the center of each grid cell to the nearest facility.
Sum the values in (2) and (3) for each tract and divide by the number of cells in each tract.

You're then left with two new variables: the average number of facilities in each tract, and the average distance to the nearest facility for each tract.

This method sounds interesting, but beyond converting the census tract vector map to a raster map using polygon to raster (and truth be told, I'm not even sure I did THAT correctly), I'm not really sure how to proceed. What tools should I be using for this?

Best Answer

Both variables are zonal means.

The average distance to the nearest facility is the zonal mean of the Euclidean distance grid (based on the facilities).
The average number of facilities is the zonal mean of a one-kilometer radius focal sum of the facilities grid. (This is merely a grid whose cell values count the number of facilities within each cell. Typically it will be an indicator grid whose values are just zeros and ones.)

These both are fast, computationally efficient calculations to perform. Moreover, great flexibility is afforded (if desired) by using a weighted focal sum in (2), which is not any more difficult to compute but can weight facilities according to their distances.

Related Solutions

[GIS] Calculation of average point distances from one another within a polygon grid

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")

Results

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.

[GIS] How to all polylines within a polygon be selected and then operated on to determine the average value from an attribute field

I doubt the null output could be due to this

Instead of using the SpatialJoin Geoprocessing too, you could try the below method.

Right Click your Polygon Layer > Joins and Relates - Select Join.

give the parameters as per the picture below. The counts and the sum of each attributes will be updated for each polygon grid

enter image description here

Best Answer

Related Solutions

[GIS] Calculation of average point distances from one another within a polygon grid

[GIS] How to all polylines within a polygon be selected and then operated on to determine the average value from an attribute field

Related Question