R语言空间建模：如何在地图上展示时间序列预测的预测区间

前置准备 #

提前整理好你的时间序列预测结果数据框，至少包含点位唯一标识、经度、纬度、预测值四个核心字段，按需安装依赖包：

install.packages(c("sf", "ncdf4", "ggplot2", "ggspatial", "dplyr"))

场景1：对接Shapefile空间底图 #

• 读取Shapefile底图

  library(sf)
  library(dplyr)
  # 读取本地Shapefile，quiet=T屏蔽冗余日志输出
  shp_base <- st_read("你的底图存储路径/底图文件名.shp", quiet = T)
  # 查看底图坐标系，后续预测点要和该坐标系保持一致
  print(st_crs(shp_base))
  

• 将预测结果转换为空间点对象

  # coords参数传入你的经纬度字段名，crs和底图保持一致
  pred_point <- st_as_sf(你的预测结果数据框, coords = c("经度字段名", "纬度字段名"), crs = st_crs(shp_base))
  

• （可选）空间关联底图属性
  如果需要将点位匹配到对应的行政单元/网格单元，可执行空间连接：

  pred_point_join <- st_join(pred_point, shp_base, join = st_within)
  

场景2：对接netCDF空间数据 #

• 读取netCDF文件并提取空间信息

  library(ncdf4)
  # 打开nc文件连接
  nc_conn <- nc_open("你的nc文件存储路径/文件名.nc")
  # 提取nc文件内的坐标、属性字段，字段名需和你的nc文件内置字段匹配
  nc_lon <- ncvar_get(nc_conn, "lon")
  nc_lat <- ncvar_get(nc_conn, "lat")
  nc_attr <- ncvar_get(nc_conn, "你需要提取的属性字段名")
  # 关闭连接
  nc_close(nc_conn)
  

• 转换为空间网格对象并匹配预测点

  # 将nc的三维数组转换为数据框
  nc_df <- expand.grid(lon = nc_lon, lat = nc_lat) %>%
    mutate(attr_val = as.vector(nc_attr))
  # 转为空间对象，坐标系默认WGS84可按需调整
  nc_sf <- st_as_sf(nc_df, coords = c("lon", "lat"), crs = 4326)
  # 近邻匹配每个预测点对应的nc网格属性
  pred_point_join <- st_join(pred_point, nc_sf, join = st_nearest_feature)
  

地图可视化实现 #

基于ggplot2的sf扩展即可快速出图，示例代码如下：

library(ggplot2)
library(ggspatial)

ggplot() +
  # 绘制空间底图，填充、描边颜色可自定义
  geom_sf(data = shp_base, fill = "#f8f9fa", color = "#adb5bd", linewidth = 0.3) +
  # 绘制预测点位，颜色映射预测值
  geom_sf(data = pred_point, aes(color = 预测值字段名), size = 2, alpha = 0.8) +
  # 自定义配色方案
  scale_color_viridis_c(option = "plasma") +
  # 添加指南针、比例尺
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering, height = unit(1, "cm"), width = unit(1, "cm")) +
  annotation_scale(location = "br", unit_category = "metric") +
  # 标签与主题调整
  labs(title = "时间序列预测结果空间分布", color = "预测值") +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.title = element_blank(), panel.grid = element_blank())

常见注意事项 #

• 所有空间数据必须保持坐标系（CRS）一致，若存在偏移可使用st_transform(目标CRS代码)转换坐标系，常用WGS84坐标系的EPSG代码为4326
• 若你有多个时间步的预测结果，可添加facet_wrap(~时间步字段名)批量生成多时间步的分布面板
• 若点位数量较多，可替换geom_sf为stat_density2d_filled绘制预测值热力图，展示空间聚集特征

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