R语言合并蒙特卡洛模拟点地图及计算KDE栅格均值

嘿，我注意到你已经用purrr::map_dfr搞定了50次蒙特卡洛模拟点的合并，但在处理KDE栅格均值的时候遇到了麻烦。先帮你把点合并的部分再确认下，然后重点解决栅格的问题！

一、模拟点合并的优化（已解决，补充细节） #

你提到用purrr::map_dfr合并点，其实可以给每个模拟结果加个iteration标识列，这样后续可视化或者分析时还能区分不同批次的模拟点，代码更清晰：

# 合并所有模拟点，并添加迭代编号列
all_points <- purrr::map_dfr(monte, ~., .id = "iteration")

得到的all_points就是包含所有50次模拟点的sf对象，直接用ggplot2或者基础plot()就能画出合并后的地图了。

二、KDE栅格均值计算的核心解决方案 #

首先得明确：计算多个栅格的均值，前提是所有栅格的范围、分辨率、CRS必须完全一致。下面分步骤给你拆解：

1. 先确保栅格参数一致 #

假设你已经生成了每个模拟对应的KDE栅格，存放在列表kde_raster_list里（比如用map(monte, ~st_kde(.x))生成）。先检查所有栅格的参数是否统一：

# 提取所有栅格的关键属性
raster_props <- purrr::map_df(kde_raster_list, function(r) {
  tibble(
    crs = st_crs(r)$wkt,
    xmin = st_bbox(r)[["xmin"]],
    xmax = st_bbox(r)[["xmax"]],
    ymin = st_bbox(r)[["ymin"]],
    ymax = st_bbox(r)[["ymax"]],
    res_x = st_resolution(r)[1],
    res_y = st_resolution(r)[2]
  )
})

# 检查是否所有栅格属性都和第一个一致
all(raster_props == raster_props[1, ])

如果返回TRUE，直接进入下一步；如果是FALSE，需要用st_warp()把所有栅格对齐到第一个栅格的参数：

# 对齐所有栅格到第一个栅格的范围、分辨率和CRS
kde_raster_aligned <- purrr::map(kde_raster_list, ~st_warp(.x, kde_raster_list[[1]]))

2. 两种计算均值的方法 #

方法一：基于sf数据框分组计算（易理解）

把每个栅格转换为带KDE值的sf数据框，然后按栅格单元（geometry）分组求均值：

# 将所有对齐后的栅格转为长格式数据框，保留迭代标识
kde_df <- purrr::map_dfr(kde_raster_aligned, function(r) {
  st_as_sf(r, as_points = FALSE) %>%  # 把栅格转为sf多边形（每个单元是一个多边形）
    mutate(kde_value = as.vector(r)) %>%  # 提取每个单元的KDE值
    select(geometry, kde_value)
}, .id = "iteration")

# 按栅格单元分组，计算均值
kde_mean_df <- kde_df %>%
  group_by(geometry) %>%
  summarise(mean_kde = mean(kde_value, na.rm = TRUE)) %>%
  ungroup()

# 可选：转回栅格格式（如果需要栅格输出）
kde_mean_raster <- st_rasterize(kde_mean_df, template = kde_raster_list[[1]])

方法二：用terra包高效计算（推荐）

如果你的栅格是terra::SpatRaster格式（现在sf生态更推荐配合terra处理栅格），可以直接用集合操作计算均值，效率更高：

library(terra)

# 将栅格列转为SpatRaster集合
kde_terra_collection <- sprc(kde_raster_list)

# 一键计算所有栅格的均值
kde_mean_raster <- mean(kde_terra_collection, na.rm = TRUE)

# 可选：转回sf格式的栅格
kde_mean_sf <- st_as_sf(kde_mean_raster)

三、完整流程示例（包含KDE生成） #

如果还没生成KDE栅格，这里给你补全整个流程的代码，确保能跑通：

library(sf)
library(tidyverse)
library(terra)
library(spatstat)  # 用于KDE计算，也可以用sf::st_kde

# 你的原始数据和模拟代码
mydf <- structure(list(x = c(555624, 572481, 703318, 700818, 571713, 559113, 731606, 604972), y = c(218959, 184051, 233180, 233603, 182307, 153136, 279015, 200216), acc.=c(120.35451, 128.74603, 74.63894, 73.96213, 93.63799, 257.49205, 53.59192, 190.78791)), .Names = c("x","y","acc."), class = "data.frame", row.names = c(NA, -8L))
mydf = st_as_sf(mydf, coords=c("x","y"), crs=21781)
n = 50 

move_point <- function(point, maxdistance){
  n_points <- length(point)
  angle_deg <- runif(n_points,1,360)
  distance <- rnorm(n_points,mean = 0,sd = maxdistance/2)
  angle_rad <- (angle_deg * pi) / (180)
  point_old <- st_coordinates(point)
  xoffset <- cos(angle_rad) * distance
  yoffset <- sin(angle_rad) * distance
  point_new <- point_old + matrix(c(xoffset,yoffset), ncol = 2)
  point_new <- point_new %>% split(1:nrow(.)) %>% map(~st_point(.x))%>% st_sfc()
  point_new
}

monte <- purrr::map(1:n, function(x){
  at <-move_point(mydf$geometry, mydf$acc.)
  st_as_sf(at, crs=21781)
})

# 1. 合并所有模拟点
all_points <- purrr::map_dfr(monte, ~., .id = "iteration")

# 2. 生成统一参数的KDE栅格列
# 先创建覆盖所有点的统一窗口（确保所有KDE栅格范围一致）
window <- st_bbox(all_points) %>%
  st_as_sfc() %>%
  as("owin")

# 生成每个模拟的KDE栅格
kde_raster_list <- purrr::map(monte, function(sf_points) {
  # 转换为spatstat的ppp对象
  ppp_points <- as(sf_points, "ppp")
  # 计算KDE（用Diggle带宽，也可以自定义sigma）
  kde <- density(ppp_points, sigma = "bw.diggle", window = window)
  # 转换为terra的SpatRaster并设置CRS
  raster(kde) %>%
    set.crs(st_crs(sf_points)$wkt)
})

# 3. 计算KDE栅格均值
kde_terra_collection <- sprc(kde_raster_list)
kde_mean_raster <- mean(kde_terra_collection, na.rm = TRUE)

# 4. 可视化：合并点 + 均值KDE
ggplot() +
  geom_sf(data = st_as_sf(kde_mean_raster), aes(fill = mean), alpha = 0.5) +
  geom_sf(data = all_points, size = 0.5, alpha = 0.3) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(title = "50次蒙特卡洛模拟点合并与KDE均值", fill = "均值KDE值")

内容的提问来源于stack exchange，提问作者uelf1