咨询：R语言中是否有适配Shapefile的等值线插值函数

你当前代码用geom_density_2d_filled生成的填充是基于点的密度分布，默认会在设定坐标范围内生成矩形填充，不会自动贴合Shapefile（毛里求斯岛屿）的边界。要实现适配Shapefile的等值线填充，需要先对离散点做空间插值生成覆盖整个Shapefile区域的栅格数据，再裁剪到Shapefile边界内，最后完成绘图。

以下是具体实现步骤和代码：

所需额外包 #

需要新增gstat（空间插值）、raster（栅格处理）两个包，先安装加载：

install.packages(c("gstat", "raster"))
library(ggplot2)
library(sf)
library(ggspatial)
library(dplyr)
library(gstat)
library(raster)

完整修改代码 #

# 读取数据
data_points <- read.csv("data_points.csv")
mauritius_shape <- st_read("mus_admbnda_adm0_2020_v2.shp")

# 重命名坐标列并过滤数据（移除重复的filter行）
data_points <- data_points %>%
  rename(lon = longitude, lat = latitude) %>%
  filter(lon >= 57.3, lon <= 57.8, lat >= -20.5, lat <= -19.95)

# 将点数据转为sf格式，确保坐标系与Shapefile一致
data_sf <- st_as_sf(data_points, coords = c("lon", "lat"), crs = st_crs(mauritius_shape))

# 创建覆盖Shapefile的栅格网格（插值目标区域）
grid <- st_make_grid(mauritius_shape, cellsize = 0.005, what = "centers") %>%
  st_sf() %>%
  st_filter(mauritius_shape)  # 仅保留Shapefile范围内的网格点

# 执行IDW空间插值（替换value为你实际要绘制的变量列名）
# 若需更精准的克里金插值，可替换为krige函数，需先拟合变异函数
interp_result <- idw(formula = value ~ 1, locations = data_sf, newdata = grid, idp = 2)

# 将插值结果转为栅格，再生成等值线多边形
r <- rasterFromXYZ(st_coordinates(interp_result) %>% cbind(interp_result$var1.pred))
contour_poly <- rasterToContour(r, nlevels = 10) %>% st_as_sf()

# 裁剪等值线多边形到Shapefile边界
contour_poly_clipped <- st_intersection(contour_poly, mauritius_shape)

# 绘图
ggplot() +
  geom_sf(data = mauritius_shape, fill = "white", color = "black") +
  geom_sf(data = contour_poly_clipped, aes(fill = level), color = NA) +
  geom_sf(data = contour_poly_clipped, color = "black", size = 0.3) +
  scale_fill_viridis_c(name = "插值值") +
  coord_sf(xlim = c(57.3, 57.8), ylim = c(-20.5, -19.95), expand = FALSE) +
  theme_minimal() +
  labs(title = "适配毛里求斯边界的等值线图", x = "经度", y = "纬度") +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering)

关键说明 #

1. 空间插值：示例用IDW（反距离加权）插值，适合离散点的空间预测；若数据存在明显空间相关性，克里金插值（krige）精度更高，需先通过variogram拟合变异函数。
2. 网格创建：st_make_grid生成覆盖Shapefile的细密网格，确保插值结果完整覆盖岛屿区域。
3. 边界裁剪：st_intersection将等值线多边形与Shapefile做交集运算，仅保留岛屿内部的填充区域，实现边界适配。
4. 变量替换：代码中的value需替换为你数据中实际要绘制等值线的变量列名（如温度、浓度等）。

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