geotheory - 11 months ago 28
R Question

# Calculating weighted polygon centroids in R

I need to calculate the centroids of a set of spatial zones based on a separate population grid dataset. Grateful for a steer on how to achieve this for the example below.

``````require(raster)
require(spdep)
require(maptools)

dat <- raster(volcano)   # simulated population data

# set consistent coordinate ref. systems and bounding boxes
proj4string(dat) <- proj4string(polys) <- CRS("+proj=longlat +datum=NAD27")
extent(dat) <- extent(polys)

# illustration plot
plot(dat, asp = TRUE)
``````

Three steps:

First, find all the cells in each polygon, return a list of 2-column matrices with the cell number and the value:

``````require(plyr) # for llply, laply in a bit...
cell_value = extract(dat, polys,cellnumbers=TRUE)
cell value
[1,]   31   108
[2,]   32   108
[3,]   33   110
[4,]   92   110
[5,]   93   110
[6,]   94   111
``````

Second, turn into a list of similar matrices but add the x and y coords:

``````cell_value_xy = llply(cell_value, function(x)cbind(x,xyFromCell(dat,x[,"cell"])))
cell value        x        y
[1,]   31   108 8.581164 14.71973
[2,]   32   108 8.669893 14.71973
[3,]   33   110 8.758623 14.71973
[4,]   92   110 8.581164 14.67428
[5,]   93   110 8.669893 14.67428
[6,]   94   111 8.758623 14.67428
``````

Third, compute the weighted mean coordinate. This neglects any edge effects and assumes all grid cells are the same size:

``````centr = laply(cell_value_xy, function(m){c(weighted.mean(m[,3],m[,2]), weighted.mean(m[,4],m[,2]))})
Now `centr` is a 2-column matrix. In your example its very close to `coordinates(polys)` so I'd make a contrived example with some extreme weights to make sure its working as expected.