如何基于IDW插值绘制秘鲁利马地区水质pH值热力图？

所需R包 #

先加载必备的空间数据处理与绘图包：

library(sf)
library(gstat)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

1. 数据准备 #

将你的经纬度-pH数据转换为sf空间对象（替换示例中的模拟数据为你的真实数据）：

# 模拟数据（替换为你的实际数据）
set.seed(123)
ph_data <- data.frame(
  lon = runif(50, -77.2, -76.8),
  lat = runif(50, -12.2, -11.8),
  ph = rnorm(50, 7.5, 0.3)
)
# 转换为WGS84坐标系（EPSG:4326）的sf对象
ph_sf <- st_as_sf(ph_data, coords = c("lon", "lat"), crs = 4326)

2. 获取利马地图边界 #

提取秘鲁利马地区的行政边界，确保与数据坐标系一致：

# 获取秘鲁省级行政边界
peru <- ne_states(country = "peru", returnclass = "sf")
# 筛选利马区域
lima <- peru[peru$name == "Lima", ]
# 统一为WGS84坐标系
lima <- st_transform(lima, crs = 4326)

3. 生成利马区域内的插值网格 #

不再生成方形网格，而是基于利马边界生成贴合区域的网格：

# 在利马边界内生成规则网格（cellsize控制精度，单位为度）
grid <- st_make_grid(lima, cellsize = 0.01, what = "centers")
# 筛选出网格中位于利马边界内的点
grid_sf <- st_sf(grid) %>% 
  st_intersection(lima) %>% 
  rename(geometry = grid)

4. 执行IDW插值 #

使用gstat包的idw函数完成插值，调整参数适配你的数据：

# IDW插值：idp为距离权重幂次（默认2，值越大近点权重占比越高）
idw_model <- idw(formula = ph ~ 1, locations = ph_sf, newdata = grid_sf, idp = 2)
# 重命名插值结果列方便后续绘图
idw_result <- idw_model %>% 
  rename(predict_ph = var1.pred)

5. 绘制叠加地图的热力图 #

用ggplot2将插值结果与利马地图边界叠加：

ggplot() +
  # 绘制IDW插值的热力填充
  geom_sf(data = idw_result, aes(fill = predict_ph, color = predict_ph), size = 0.1) +
  # 叠加利马地图边界
  geom_sf(data = lima, fill = NA, color = "black", size = 1) +
  # 设置pH值的颜色渐变
  scale_fill_viridis_c(option = "plasma", name = "pH值") +
  scale_color_viridis_c(option = "plasma", name = "pH值") +
  # 图表标题与主题
  labs(title = "利马区域水质pH值IDW插值热力图", x = "经度", y = "纬度") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

可选优化：获取更精细的利马边界 #

如果需要更详细的城市边界，可使用osmdata包从OpenStreetMap获取：

library(osmdata)
# 基于现有利马边界生成查询范围
lima_bbox <- st_bbox(lima)
# 查询利马城市边界
lima_osm <- opq(lima_bbox) %>%
  add_osm_feature(key = "admin_level", value = "6") %>%
  osmdata_sf()
# 提取利马边界对象
lima_boundary <- lima_osm$osm_multipolygons %>%
  filter(name == "Lima")
# 后续网格生成与插值可替换为lima_boundary

内容的提问来源于stack exchange，提问作者Mirko Caja Ventura