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

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

• why use the density rather than actual counts? Mar 29, 2011 at 21:37
• @mindless Typically when a count is associated with a polygon it is an aggregate statistic for the area covered by that polygon. Ideally we would distribute the population throughout the polygon into precisely the right locations, but usually we cannot because we don't know exactly where the people live. Using the density assumes all locations within the polygon are equally likely. It turns out (mathematically) that you could assign the entire count to the centroid (not "center") of the polygon and you would get the same result. (That's why your question is such a good one. :-) Mar 29, 2011 at 21:55
• I see - I'm actually using a raster that tries to do exactly what you mentioned - distribute the population precisely over the entire area. GPW, GRUMP, UNEP, LandScan are examples of this. Since I do have the actual (or predicted rather) population at some resolution, I don't use the density. Now, assuming I did what you did, where there is only one value of the polygon, doesn't that result in simply the centroid being the typical geometric center? Final twist, if I get a weighted centroid that's outside of its associated polygon, does this entail data issues, or is this possible? Mar 29, 2011 at 22:06
• @mindless (1) If the count in the polygon is 1, converting to a density effectively palces it at the centroid. (2) There is no problem with centroids lying outside their polygons. NB: "centroid" in this discussion always means the physical center-of-mass definition, not any of the other GIS alternatives (such as those used to place labelpoints). (3) (at the beginning): when you have a grid of population counts, you effectively are working with a density already (because all cell areas are the same). Mar 29, 2011 at 22:15
• @mindless So far our remarks have focused on finding the population center. For almost all other related calculations, such as estimates of dispersion around that center, placing all polygon counts at the polygon centroids can bias the dispersion estimates, depending on how those estimates are made. This is where converting to a density grid is the best approach. Mar 29, 2011 at 22:17

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.