How do I calculate a weighted centroid with respect to a raster, e.g. population weighted centroid

Question

In ArcGIS 10.0, how do I calculate the weighted centroid of a collection of polygons, where I want the weight to be a raster layer. In particular I am using a raster of population counts and would like to generate the population weighted centroid as a result.

Answer 1

There's a neat trick: you can create grids of row and column indexes by performing FlowAccumulation calculations on constant grids, as in FlowAccumulation(1) (column indexes, starting at 0, increasing to the right) and FlowAccumulation(64) (row indexes, starting at 0, increasing upwards).

For your purposes that's good enough. When you really need $$XMap and $$YMap, rescale the row/column indexes by the cellsize and shift by the origin (plus, perhaps, another half a cellsize in both directions to obtain coordinates of cell centers).

No scripting needed (and much faster in execution, too)!

Now to answer your question: convert the polygons to grid format after first computing the population per unit area. Call this grid [density], because it's a population density. Compute [row index] and [column index] grids as shown above. To obtain the x-coordinates of the centroids, divide the zonal sums of [density]*[column index] by the zonal sums of [density] (there's one of each per polygonal zone). Do a similar operation with the row indexes to obtain the y-coordinates of the centroids. If desired, convert these centroids (which are averaged row/column indexes) to coordinates by scaling by the cellsize and adding the coordinates of the origin (plus one-half the cellsize).

Answer 2

This post on ESRI forums provides the logic on using $$XMap, $$YMap, some map algebra, and zonal summaries to come up with the x and y coordinates for each weighted centroid.

The following code is my attempt at replicating the now removed $$XMap and $$YMap grid expressions in python in ArcGIS 10.0:

import arcpy
import numpy

# reference to the weight raster
weightRas = arcpy.Raster("weighted_raster")

# dimensions for the arrays
numColumns = weightRas.width
numRows = weightRas.height

# create arrays that will hold coordinate values
xCoords = numpy.zeros( (numRows,numColumns), dtype='float32' )
yCoords = numpy.zeros( (numRows,numColumns), dtype='float32' )

# each cell center is offset by mean width/height
xOffset = weightRas.meanCellWidth
yOffset = weightRas.meanCellHeight

# set direction of the offsets
if( (weightRas.extent.upperRight.X - weightRas.extent.upperLeft.X) < 0 ) :
    xOffset = -xOffset
if( (weightRas.extent.lowerLeft.Y - weightRas.extent.upperLeft.Y) < 0 ) :
    yOffset = -yOffset

# the start of x coords for the raster is the center of the first pixel in the
# upper left corner, which is the upper left extent's x coord plus half the offset
x = weightRas.extent.upperLeft.X + (xOffset / 2)
y = weightRas.extent.upperLeft.Y + (yOffset / 2)

# now fill the arrays
for i in range(0,numColumns):
    xCoords[:,i] = x
    x = x + xOffset

for i in range(0,numRows):
    yCoords[i] = y
    y = y + yOffset

# finish by converting back to raster, specifying the lower left corner of the
# raster as well as cell size for x and y, which we pull from our weight raster
xMap = arcpy.NumPyArrayToRaster( xCoords, arcpy.Point(weightRas.extent.lowerLeft.X, weightRas.extent.lowerLeft.Y), weightRas, weightRas )
yMap = arcpy.NumPyArrayToRaster( yCoords, arcpy.Point(weightRas.extent.lowerLeft.X, weightRas.extent.lowerLeft.Y), weightRas, weightRas )

At this point I have rasters of both x and y coordinates and can continue on using the logic in the above link.