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

我的问题是：有没有办法为这些点创建一个平滑的地图，并且平滑的地图完全填充在外部边界内（没有白色 space 也没有“侵入”）？

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 + 添加 window 到 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 package (v2.0.0)}

于 2021-08-04 创建

平滑值是使用 spatstat 包估算的，它们非常适合原始边界。如果您确实需要填充微小的间隙，请增加 dimyx 的值。查看 ?density.ppp 及其中的参考资料以获取更多详细信息。

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

Smoothed density maps for points in using sf within a given boundary in R

r

smoothing

ggplot2

kernel-density

sf