r - 多边形内的人口密度

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)

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

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

谢谢你的帮助！

标签： rspatialrastershapefile

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包，我明白了。