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