多边形内的人口密度

Question

因此，我对R中的栅格数据包有一些疑问。我有一个栅格，其中每个栅格点的人口估计都在其中。我也有一个shapefile，其中包含区域的多边形。我想找出每个区域内人口密度最高的社区的坐标。假设每个邻域都是5 x 5网格点的同质正方形。

以下玩具示例模仿了我的问题。

library(raster)
library(maptools)

set.seed(123)

data(wrld_simpl)

wrld_simpl <- st_as_sf(wrld_simpl)
contr_c_am <- wrld_simpl %>%
  filter(SUBREGION ==13) %>%
  filter(FIPS != "MX") %>%
  select(NAME) 


# Create a raster of population (sorry for the bad example spatial distribution)
r <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90, res=0.1)
values(r) <- runif(ncell(r), 0, 100)

# keep only raster around the region of interest
r_small <- crop(r, extent(contr_c_am))
plot(r_small)
plot(st_geometry(contr_c_am), add = T)

raster_contr_c_am <- rasterize(contr_c_am, r)

raster_contr_c_am是人口网格，并且该区域的名称另存为属性。

以某种方式，我只需要过滤一个区域中的网格点，并可能使用诸如focus（）之类的函数来查找附近的总体。

focal(raster_contr_c_am, matrix(1,5,5),sum, pad = T, padValue = 0)

然后，我需要找到每个区域内哪个网格点的值最高，并保存其坐标。

我希望我的解释不要太混乱，

感谢您的帮助！

Answer 1

这里是一个示例，它在定义区域的形状上进行迭代，然后使用区域内的栅格值和focal()函数来找到最大值。

library(raster)
library(maptools)
library(sf)
library(dplyr)

set.seed(123)

data(wrld_simpl)

wrld_simpl <- st_as_sf(wrld_simpl)
contr_c_am <- wrld_simpl %>%
    filter(SUBREGION ==13) %>%
    filter(FIPS != "MX") %>%
    select(NAME) 

# Create a raster of population (sorry for the bad example spatial distribution)
r <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90, res=0.1)
values(r) <- runif(ncell(r), 0, 100)

# keep only raster around the region of interest
r_small <- crop(r, extent(contr_c_am))
raster_contr_c_am <- rasterize(contr_c_am, r_small)

# function to find the max raster value using focal
# in a region 
findMax <- function(region, raster) {
    tt <- trim((mask(raster, region))) # focus on the region
    ff <- focal(tt, w=matrix(1/25,nc=5,nr=5))
    maximumCell <- which.max(ff) # find the maximum cell id
    maximumvalue <- maxValue(ff) # find the maximum value
    maximumx <- xFromCell(ff, maximumCell) # get the coordinates
    maximumy <- yFromCell(ff, maximumCell)
    # return a data frame
    df <- data.frame(maximumx, maximumy, maximumvalue)
    df
}

numberOfShapes <- nrow(contr_c_am)
ll <- lapply(1:numberOfShapes, function(s) findMax(region = contr_c_am[s,], raster = r_small))
merged <- do.call(rbind, ll)
maxpoints <- st_as_sf(merged, coords=c('maximumx', 'maximumy'), crs=crs(contr_c_am))

library(mapview) # optional but nice visualization - select layers to see if things look right
mapview(maxpoints) + mapview(r_small) + mapview(contr_c_am)

我制作了一个sf对象，以便可以与其他空间对象一起绘制。使用mapview包，我得到了。