如何为外部点匹配克里金插值生成的等时线/等值带

我在硕士论文中开展克里金空间插值工作，通过普通克里金插值绘制农业从近东向欧洲传播的等时线地图，并用ggplot完成可视化。当前需要给未参与插值的外部红色点数据集新增一列，标记每个点所属的等时线区间（例如巴利阿里群岛的红点对应(7000,7500]区间）。

尝试将ggplot中的等时线转为多边形后，用st_intersect()或st_contains()进行点匹配，但无法从ggplot地图中生成等时线多边形，尝试过ggplot_build提取图层转sf对象但未成功。以下是使用公开数据复现的完整代码：

library(gstat)
library(sf)
library(readr)
library(tidyverse)
library(rnaturalearth)

##Data
no2 <- read_csv(system.file("external/no2.csv", package = "gstat"), show_col_types = FALSE)
no2_sf <- st_as_sf(no2, crs = "OGC:CRS84", coords = c("station_longitude_deg", "station_latitude_deg")) 
no2_sf <- st_transform(no2_sf, 32632)


##Kriging
#Interpolation raster
no2_bbox <- st_bbox(no2_sf)
cell_size <- 10000
x <- seq(no2_bbox$xmin, no2_bbox$xmax, by=cell_size)
y <- seq(no2_bbox$ymin, no2_bbox$ymax, by=cell_size)
no2_grid <- expand.grid(x=x, y=y)
no2_grid$tmp <- 1
plot(no2_grid$x, no2_grid$y, pch=19, cex=0.1)
no2_grid <- st_as_sf(no2_grid, coords = c("x","y"), crs = st_crs(no2_sf))
st_crs(no2_grid) <- st_crs(no2_sf)

#Variogram
no2_sample_variogram <- gstat::variogram(NO2~1, no2_sf)
plot(no2_sample_variogram, plot.numbers = TRUE)
no2_model_variogram <- vgm(psill = 16, "Exp", range = 200000, nugget = 1)
plot(no2_sample_variogram, no2_model_variogram)
no2_fit_variogram <- fit.variogram(no2_sample_variogram, no2_model_variogram)
plot(no2_sample_variogram, no2_fit_variogram)

#Ordinary Kriging
no2_sf <- no2_sf[!duplicated(no2_sf$geometry),] #check which observation were removed
no2_kriging <- gstat::krige(NO2~1, no2_sf, no2_grid, no2_fit_variogram)

no2_kriging$x <- st_coordinates(no2_kriging)[,1]
no2_kriging$y <- st_coordinates(no2_kriging)[,2]

##Ggplot
points <- data.frame(x = runif(10, min = no2_bbox$xmin, max = no2_bbox$xmax), y = runif(10, min = no2_bbox$ymin, max = no2_bbox$ymax))
points <- st_as_sf(points, coords = c("x","y"), crs = st_crs(32632))
  
germany  <- ne_countries(scale = "medium", returnclass = "sf", country = "Germany")
germany <- st_transform(germany, 32632)

ggplot()+  geom_contour_filled(data = no2_kriging, aes(x = x, y = y, z=var1.pred))+  geom_sf(data = germany,  fill = "transparent", color = "black")+  geom_sf(data = points, size = 0.5, color = "red")

不要从ggplot图层反向提取等值带，直接基于克里金插值结果生成等值带多边形，再与外部点做空间匹配，步骤如下：

1. 安装并加载isoband包 #

isoband是geom_contour_filled背后的核心包，可以直接从网格数据生成等值带多边形，确保结果与ggplot可视化完全一致：

install.packages("isoband")
library(isoband)

2. 从克里金插值结果生成等值带sf对象 #

将插值后的网格数据转换为矩阵，再生成等值带多边形：

# 提取网格的x、y坐标和插值预测值矩阵
x_vals <- unique(no2_kriging$x)
y_vals <- unique(no2_kriging$y)
z_matrix <- matrix(no2_kriging$var1.pred, nrow = length(y_vals), ncol = length(x_vals))

# 设置等值带区间（与ggplot的binwidth保持一致，这里用自动计算，也可手动指定）
# 手动指定区间示例：breaks <- seq(min(no2_kriging$var1.pred), max(no2_kriging$var1.pred), by = 5)
breaks <- isobreaks(z_matrix, binwidth = 5) # binwidth需和geom_contour_filled的参数匹配

# 生成等值带多边形并转为sf对象
iso_polygons <- isobands(x_vals, y_vals, z_matrix, breaks, breaks[-1])
iso_sf <- iso_to_sfg(iso_polygons) %>% 
  st_sfc(crs = st_crs(no2_kriging)) %>% 
  st_sf() %>% 
  mutate(interval = names(iso_polygons)) # 添加区间标签，格式与ggplot一致

3. 空间匹配外部点与等值带区间 #

用st_join()完成点与多边形的空间连接，自动为每个点添加所属的区间：

# 空间连接，保留点的所有属性，新增interval列
points_with_interval <- st_join(points, iso_sf, join = st_intersects)

# 查看匹配结果
head(points_with_interval)

4. 验证匹配准确性（可选） #

将生成的等值带多边形和带区间的点一起绘制，确认匹配正确：

ggplot() +
  geom_sf(data = iso_sf, aes(fill = interval)) +
  geom_sf(data = germany, fill = "transparent", color = "black") +
  geom_sf(data = points_with_interval, color = "red", size = 1) +
  labs(fill = "NO2 Interval")

关键注意事项 #

• 确保breaks的binwidth与geom_contour_filled的参数完全一致，保证区间划分完全匹配
• 该方法生成的区间标签格式（如(10,15]）与ggplot完全相同，直接满足你的需求
• 避免了从ggplot图层提取数据的繁琐操作，结果更可靠

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