R语言akima包interp()插值结果查询指定(x,y)对应z值的方法

你用akima::interp()生成的是规则格网插值结果，核心结构为：长度为nx的升序经度格网向量、长度为ny的升序纬度格网向量、存储每个格网点高程的nx*ny矩阵。要实现z_i = f(x_i,y_i)的任意点查询，不要用单维度匹配最近值的写法——那种写法完全忽略纬度方向的位置差异，结果误差极大。

最稳妥的实现是基于规则格网做双线性插值，不需要额外依赖包，直接写通用查询函数即可：

# 定义任意点高程查询函数
get_z <- function(target_x, target_y, interp_result) {
  # 提取格网信息
  x_grid <- interp_result$x
  y_grid <- interp_result$y
  z_matrix <- interp_result$z
  
  # 边界判断，超出插值范围直接返回NA
  if (target_x < min(x_grid) | target_x > max(x_grid) | 
      target_y < min(y_grid) | target_y > max(y_grid)) {
    return(NA_real_)
  }
  
  # 定位目标点在x、y格网中相邻的两个格网点位置
  x_idx <- findInterval(target_x, x_grid)
  y_idx <- findInterval(target_y, y_grid)
  
  # 处理刚好落在格网最右/最上边界的情况
  if (x_idx == length(x_grid)) x_idx <- x_idx - 1
  if (y_idx == length(y_grid)) y_idx <- y_idx - 1
  
  # 提取相邻四个格网点的坐标与高程
  x1 <- x_grid[x_idx]
  x2 <- x_grid[x_idx + 1]
  y1 <- y_grid[y_idx]
  y2 <- y_grid[y_idx + 1]
  z11 <- z_matrix[x_idx, y_idx]
  z12 <- z_matrix[x_idx, y_idx + 1]
  z21 <- z_matrix[x_idx + 1, y_idx]
  z22 <- z_matrix[x_idx + 1, y_idx + 1]
  
  # 双线性插值计算目标点z值
  x_weight <- (target_x - x1) / (x2 - x1)
  y_weight <- (target_y - y1) / (y2 - y1)
  z_interp <- (1 - x_weight) * (1 - y_weight) * z11 +
    x_weight * (1 - y_weight) * z21 +
    (1 - x_weight) * y_weight * z12 +
    x_weight * y_weight * z22
  
  return(z_interp)
}

使用示例 #

如果你按自己的习惯把插值结果存在s.smooth对象中，调用时把传入的interp_result参数替换为s.smooth即可。

library(akima)
# 你原来的插值代码
s <- interp(x, y, z, nx=100, ny=100)

单点查询 #

直接传入目标经纬度和插值结果对象即可：

# 查询(116.39, 39.90)位置的高程
z_result <- get_z(target_x = 116.39, target_y = 39.90, interp_result = s)

批量查询 #

如果你有一个存储了多组待查询坐标的数据框，列名分别为x_i、y_i，可以直接批量生成对应高程列：

# 待查询坐标数据框示例
points <- data.frame(
  x_i = c(116.39, 116.41, 116.42),
  y_i = c(39.90, 39.92, 39.89)
)

# 批量计算插值高程
points$z_i <- mapply(
  FUN = get_z,
  target_x = points$x_i,
  target_y = points$y_i,
  MoreArgs = list(interp_result = s)
)

补充说明 #

• 函数自带边界判断，如果待查询坐标超出原始插值的经纬度范围，会返回NA，不建议随意做范围外的高程外推，结果无实际参考价值
• 如果需要更高的查询精度，可以适当调大interp()里的nx、ny参数，格网密度越高，插值查询的结果越平滑，对应计算耗时也会小幅上升
• 双线性插值的结果是连续平滑的，和interp()生成的曲面完全适配，比最近邻匹配的结果精度高一个量级

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