在 R 中给定边界内使用 sf 的点的平滑密度图

Question

我正在尝试为 R 中的多个点创建平滑地图，但我在这里没有找到完美的解决方案。

library(mapchina)
library(sf)
library(dplyr)
library(ggplot2)

# Create some sample data
sf_beijing = china %>% 
  filter(Code_Province == '11') %>% 
  st_transform(4326)

sf_points = data.frame(
  lat = c(39.523, 39.623, 40.032, 40.002, 39.933, 39.943, 40.126, 40.548),
  lon = c(116.322, 116, 116.422, 116.402, 116.412, 116.408, 116.592, 116.565)
) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Plot the boundary for Beijing and the points
ggplot() +
  geom_sf(data = sf_beijing, fill = NA) + 
  geom_sf(data = sf_points, color = 'red') + 
  theme_test()

此外，我发现这个解决方案可以为点创建平滑的地图sf。该解决方案的问题是平滑后的地图没有完全填满北京边界，并且一些平滑部分超出了边界。

ggplot() +
  stat_density_2d(data = sf_points, 
                  mapping = aes(x = purrr::map_dbl(geometry, ~.[1]),
                                y = purrr::map_dbl(geometry, ~.[2]),
                                fill = stat(density)),
                  geom = 'tile',
                  contour = FALSE,
                  alpha = 0.8) +
  geom_sf(data = sf_beijing, fill = NA) + 
  geom_sf(data = sf_points, color = 'red') + 
  scale_fill_viridis_c(option = 'magma', direction = -1) +
  theme_test()
ggsave('p1.png', width = 7, height = 8)

我的问题是：有没有办法为这些点创建平滑地图，并且平滑地图完美填充外部边界（没有空白区域，也没有“侵入”）？

Answer 1

我想提出以下方法。它相当复杂，可能有更有效的解决方案，但我认为它有效。

加载包

library(mapchina)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 2.2-2
#> Loading required package: spatstat.core
#> Loading required package: nlme
#> Loading required package: rpart
#> spatstat.core 2.3-0
#> Loading required package: spatstat.linnet
#> spatstat.linnet 2.3-0
#> 
#> spatstat 2.2-0       (nickname: 'That's not important right now') 
#> For an introduction to spatstat, type 'beginner'
library(ggplot2)

创建一个多边形和一些示例数据。请注意，我设置了预计的 CRS，因为包需要它spatstat（见下文）。

sf_beijing = china %>%
  dplyr::filter(Code_Province == '11') %>%
  st_transform(32650)

sf_points = data.frame(
  lat = c(39.523, 39.623, 40.032, 40.002, 39.933, 39.943, 40.126, 40.548),
  lon = c(116.322, 116, 116.422, 116.402, 116.412, 116.408, 116.592, 116.565)
) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(32650)

将点转换为ppp对象。检查?ppp并参考其中的更多详细信息。

ppp_points <- as.ppp(sf_points)

转换sf_beijing为owin+ 添加窗口到ppp_points. 查看?Window更多详细信息。

Window(ppp_points) <- as.owin(sf_beijing)

阴谋

par(mar = rep(0, 4))
plot(ppp_points, main = "")

平滑点

density_spatstat <- density(ppp_points, dimyx = 256)

转换density_spatstat成星星对象。检查https://r-spatial.github.io/stars/index.html了解更多详细信息。

density_stars <- stars::st_as_stars(density_spatstat)
#> Registered S3 methods overwritten by 'stars':
#>   method            from
#>   st_crs.SpatRaster sf  
#>   st_crs.SpatVector sf

转换density_stars成sf对象

density_sf <- st_as_sf(density_stars) %>%
  st_set_crs(32650)

阴谋

ggplot() +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  scale_fill_viridis_c() +
  geom_sf(data = st_boundary(sf_beijing)) +
  geom_sf(data = sf_points, size = 2, col = "black")

^{由reprex 包(v2.0.0)创建于 2021-08-04}

平滑值是使用该spatstat包估计的，并且它们非常适合原始边界。dimyx如果您确实需要填补微小的空白，请增加 的值。检查?density.ppp其中的参考文献以获取更多详细信息。