在R中计算加权多边形质心

Question

我需要根据单独的人口网格数据集计算一组空间区域的质心。感谢您指导如何通过以下示例实现这一目标。

提前致谢。

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

dat <- raster(volcano)   # simulated population data
polys <- readShapePoly(system.file("etc/shapes/columbus.shp",package="spdep")[1])

# 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)
plot(polys, add = TRUE)

Answer 1

三个步骤：

首先，找到每个多边形中的所有单元格，返回一个包含单元格编号和值的 2 列矩阵列表：

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

其次，变成类似矩阵的列表，但添加 x 和 y 坐标：

cell_value_xy = llply(cell_value, function(x)cbind(x,xyFromCell(dat,x[,"cell"])))
head(cell_value_xy[[1]])
     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

第三，计算加权平均坐标。这忽略了任何边缘效应并假设所有网格单元的大小相同：

centr = laply(cell_value_xy, function(m){c(weighted.mean(m[,3],m[,2]), weighted.mean(m[,4],m[,2]))})
head(centr)
            1        2
[1,] 8.816277 14.35309
[2,] 8.327463 14.02354
[3,] 8.993655 13.82518
[4,] 8.467312 13.71929
[5,] 9.011808 13.28719
[6,] 9.745000 13.47444

现在centr 是一个 2 列矩阵。在您的示例中，它非常接近 coordinates(polys)，因此我会制作一个带有一些极端权重的人为示例，以确保其按预期工作。

Answer 2

另一种选择。

我喜欢它的紧凑性，但它可能只有在您相当熟悉整个 raster 函数系列的情况下才有意义：

## Convert polygons to a raster layer
z <- rasterize(polys, dat)

## Compute weighted x and y coordinates within each rasterized region
xx <- zonal(init(dat, v="x")*dat, z) / zonal(dat,z)
yy <- zonal(init(dat, v="y")*dat, z) / zonal(dat,z)

## Combine results in a matrix
res <- cbind(xx[,2],yy[,2])
head(res)
#          [,1]     [,2]
# [1,] 8.816277 14.35309
# [2,] 8.327463 14.02354
# [3,] 8.993655 13.82518
# [4,] 8.467312 13.71929
# [5,] 9.011808 13.28719
# [6,] 9.745000 13.47444

Answer 3

Spacedman 和 Josh 的答案真的很棒，但我想分享另外两个相对快速和简单的替代方案。

library(data.table)
library(spatialEco)
library(raster)
library(rgdal)

使用data.table 方法：

# get centroids of raster data
  data_points <- rasterToPoints(dat, spatial=TRUE)

# intersect with polygons
  grid_centroids <- point.in.poly(data_points, polys)

# calculate weighted centroids
  grid_centroids <- as.data.frame(grid_centroids)
  w.centroids <- setDT(grid_centroids)[, lapply(.SD, weighted.mean, w=layer), by=POLYID, .SDcols=c('x','y')]

使用wt.centroid{spatialEco}：

  # get a list of the ids from each polygon
    poly_ids <- unique(grid_centroids@data$POLYID)

  # use lapply to calculate the weighted centroids of each individual polygon
    w.centroids.list <- lapply(poly_ids, function(i){wt.centroid( subset(grid_centroids, grid_centroids@data$POLYID ==i)
                                                                  , 'layer', sp = TRUE)} )

Answer 4

下面是我自己不太优雅的解决方案。给出与 Spacedman 和 Josh 完全相同的结果。

# raster to pixels
p = rasterToPoints(dat) %>% as.data.frame()
coordinates(p) = ~ x + y
crs(p) = crs(polys)

# overlay pixels on polygons
ol = over(p, polys) %>% mutate(pop = p$layer) %>% cbind(coordinates(p)) %>% 
  filter(COLUMBUS_ %in% polys$COLUMBUS_) %>%     # i.e. a unique identifier
  dplyr::select(x, y, pop, COLUMBUS_) %>% as_data_frame()

# weighted means of x/y values, by pop
pwcs = split(ol, ol$COLUMBUS_) %>% lapply(function(g){
  data.frame(x = weighted.mean(g$x, g$pop), y = weighted.mean(g$y, g$pop))
}) %>% bind_rows() %>% as_data_frame()