如何修改R代码从ERA5 Land NetCDF中提取最近非缺失值？

方案一：寻找最近的非缺失网格点（低算力优先） #

原代码直接取目标经纬度的最近网格点，但遇到NA时无法返回有效数据。我们可以先筛选出所有有有效值的网格点，再计算目标点到这些点的空间距离，最终选取距离最近的点提取数据。

具体步骤：

• 读取NetCDF文件中的经度、纬度及目标变量的完整空间场
• 筛选出变量值非NA的网格点及其坐标
• 计算目标点到每个有效网格点的空间距离
• 找到距离最小的网格点，提取对应数据

修改后的R代码：

library(tidyverse)
library(ncdf4)
library(geosphere) # 用于计算高精度球面距离，若用平面距离可省略

# 打开NetCDF文件
nc <- nc_open("./Weather/2019_07_18.nc")

# 目标经纬度
target_lat <- 14.67543
target_lon <- -17.4484

# 读取完整的经纬度和变量场（以t2m为例，先取第一个时间步）
lons <- nc$dim$longitude$vals
lats <- nc$dim$latitude$vals
t2m_field <- ncvar_get(nc, varid = "t2m", start = c(1,1,1), count = c(-1,-1,1))

# 将数据整理为数据框，保留非NA的点
valid_points <- expand.grid(lon = lons, lat = lats) %>%
  mutate(t2m = as.vector(t2m_field)) %>%
  filter(!is.na(t2m))

# 计算目标点到每个有效点的球面距离（单位：米）
valid_points <- valid_points %>%
  mutate(distance = distHaversine(c(target_lon, target_lat), cbind(lon, lat)))

# 找到距离最近的点
closest_point <- valid_points %>% slice_min(distance, n = 1)

# 提取该点的全时间序列t2m数据
t2m_result <- ncvar_get(nc, varid = "t2m",
                        start = c(which(lons == closest_point$lon),
                                  which(lats == closest_point$lat),
                                  1),
                        count = c(1,1,-1))

# 关闭NetCDF文件
nc_close(nc)

# 查看结果
t2m_result

说明：

• 若不需要高精度球面距离，可替换distHaversine为平面距离计算：sqrt((lon - target_lon)^2 + (lat - target_lat)^2)，省去geosphere包的依赖
• 批量处理多个目标点时，可将上述逻辑封装为函数，循环或向量化处理

方案二：空间插值填充缺失值 #

如果需要更平滑的结果，可对缺失区域进行空间插值。以下是反距离加权插值的示例：

library(tidyverse)
library(ncdf4)
library(gstat)
library(sp)

nc <- nc_open("./Weather/2019_07_18.nc")

target_lat <- 14.67543
target_lon <- -17.4484

lons <- nc$dim$longitude$vals
lats <- nc$dim$latitude$vals
t2m_field <- ncvar_get(nc, varid = "t2m", start = c(1,1,1), count = c(-1,-1,1))

# 整理为空间数据格式
valid_points <- expand.grid(lon = lons, lat = lats) %>%
  mutate(t2m = as.vector(t2m_field)) %>%
  filter(!is.na(t2m)) %>%
  coordinates(~lon+lat)

# 定义反距离加权插值模型，nmax控制参与插值的邻近点数量
idw_model <- gstat(formula = t2m~1, data = valid_points, nmax = 10)

# 对目标点进行插值
target_point <- data.frame(lon = target_lon, lat = target_lat) %>% coordinates(~lon+lat)
t2m_interpolated <- predict(idw_model, newdata = target_point)$var1.pred

nc_close(nc)

# 查看插值结果
t2m_interpolated

说明：

• 插值方法需要加载额外空间分析包，算力消耗高于方案一，适合需要连续场或批量填充缺失值的场景
• 可根据需求调整插值参数（如nmax）优化结果

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