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

Question

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

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()

另外，我发现this solution为 sf 点创建平滑 map 。该方案的问题是平滑后的 map 没有完全填满北京边界，并且有些平滑部分超出了边界。

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)

我的问题是:有没有办法为这些点创建平滑 map ，并且平滑 map 完美填充外部边界(没有空白，也没有“侵入”)？

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)

将密度_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

将密度星转换为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 package 于 2021 年 8 月 4 日创建(v2.0.0)}

平滑值是使用 spatstat 包估计的，并且它们非常适合原始边界。如果您确实需要填补微小的间隙，请增加 dimyx 的值。检查?密度.ppp 和其中的引用资料以获取更多详细信息。