我正在尝试使用R将我的shapefile从LHS“转换”为RHS，如下所示 .

enter image description here

原始数据（如LHS所示）取自2016年加拿大人口普查 . 它包括加拿大不列颠哥伦比亚省萨里市的所有传播区域（最小的基地之一） . 形状文件可以下载here .

将其转换为RHS的基本原理是因为典型的区域插值函数，例如 st_interpolate_aw 来自 sf 假设值在单位面积上是常数 . 我希望为不同的地理单位插入人口普查数据，而我用来做这个的包使用 st_interpolate_aw . 结果，有人建议我计算shapefile的"population ecumene" . 发给我RHS shapefile的人在GIS中做了他的分析，他描述了他的步骤如下：

  • 使用内核插补器估算shapefile上的人口密度度量

  • 提取具有合理密度水平的区域（例如50 ppl / km ^ 2）

  • 将这些提取的区域转换为整数掩码

  • 将蒙版转换为shapefile，并将其与原始shapefile“掩盖”（交叉） .

虽然这些步骤看起来非常简单，但在R中实现它是非常具有挑战性的 . 我尝试了各种空间插值/平滑技术，包括简单的平滑，idw，kriging，没有运气复制或实现结果 . 我最接近的是来自 tmaptools 库的 smooth_map 函数，如下所示 .

library(rgdal)
library(GISTools)
library(tmaptools)

## Read-in file, add `Pop_Den` column for population density
Surrey_census16_geom.sp <- readOGR("./Surrey_DA_geometry_16.shp"),
                        layer = "Surrey_DA_geometry_16")
Surrey_census16_geom.sp$Pop_Den <- Surrey_census16_geom.sp$Population / Surrey_census16_geom.sp$Shape.Area
## smooth_map doesn't seem to play well with `longlat` projections
Surrey_census16_proj.sp <- set_projection(Surrey_census16_geom.sp, projection = "eck4")
x <- smooth_map(Surrey_census16_proj.sp, "Pop_Den", unit.size = 100, smooth.raster =  T, extracting.method="grid", to.Raster = T)
xr <- x$raster

但是如果我尝试使用像 xr>50 或 focal(xr >= 50, w, sum, na.rm = T, pad = T ) 这样的代码从光栅中进行"extract"（其中w被简单地定义为3 * 3矩阵1，因为我符合该标准，而不是单个单元格 .

无论如何，我真的很感激

1）......任何人都可以指出我应该做什么来通过人口密度进行核密度估计并根据阈值提取网格/栅格单元（即修复上面的代码）

2）......任何人都可以推荐一种替代方法来实现基于人口密度“缩小”密度多边形区域的目标，以实现像RHS这样的结果

（注意：RHS图是通过比较每个多边形的总人口数与阈值来创建的 . 这里我希望按人口密度进行比较，因此实际结果可能不同）

r geospatial spatial raster