如何删除netCDF文件中的NA值以匹配海岸港口与最近有效波高网格点

我正在处理气象相关的波高数据，现有存储浪高信息的netCDF文件来自ERA5数据存储库，我的目标是将沿海港口与数据中最近的网格点匹配，后续在R中开展分析。申请数据时我已经指定了研究区域的经纬度边界，目前我使用ncdf4包的相关函数开展分析，所用代码如下：

当前运行代码 #

# 打开netCDF文件连接
nc <- nc_open("swell_bilmap.nc")

# 提取经纬度
lat <- ncvar_get(nc,'lat')
lon <- ncvar_get(nc,'lon')
dim(lat);dim(lon)

# 提取时间
t <- ncvar_get(nc, "time")

# 查看时间单位
ncatt_get(nc,'time')

# 将小时转换为日期时间格式，lubridate包的as_datetime()函数需要秒级输入
timestamp <- as_datetime(c(t*60*60),origin="1900-01-01")

# 提取浪高变量shts
data_shts <- ncvar_get(nc,"shts")

# 关闭netCDF连接
nc_close(nc)


# 生成所有经纬度网格组合
lonlat <- expand.grid(lon=lon,lat=lat)

# 将网格坐标转为空间对象
coord_brookings <- st_as_sf(lonlat,coords=c("lon","lat"))%>%
  st_set_crs(4326)

# 生成港口点空间对象，PORTS为存储港口经纬度的数据集
psj_brookings <- st_point(c(PORTS[1,5],PORTS[1,4]))%>%
  st_sfc()%>%
  st_set_crs(4326)


# 计算每个网格到港口的距离
coord_brookings <- mutate(coord_brookings,dist=st_distance(coord_brookings,psj_brookings))

# 生成和数据维度匹配的距离矩阵
dist_mat_brookings_shts <- matrix(coord_brookings$dist,dim(data_u)[-3])

# 提取最近网格的行列索引
mat_index_brookings_shts <- as.vector(arrayInd(which.min(dist_mat_brookings_shts), dim(dist_mat_brookings_shts)))

# 提取最近网格的时间序列
df_brookings_shts <- data.frame(shts=data_shts[mat_index_brookings_shts[1],mat_index_brookings_shts[2],],time=timestamp)

遇到的问题 #

上述代码可以输出包含最近网格点各时间步对应变量值的data frame。浪高数据中未与海洋主体重合的网格对应值为NA，当尝试匹配港口点到最近网格时，返回的结果为NA而非有效值，因为要匹配的港口点位于陆面，可视化时无颜色的区域并非没有数据，只是值为NA所以无展示，运行代码后得到的结果全是NA，仅需要保留非NA值的有效网格用于匹配。
刚接触netCDF文件，优先希望使用R实现删除netCDF中NA值的需求，也可以接受CDO、ArcGIS Pro的实现方案。

R方案为优先推荐实现方式

R实现方案 #

不需要修改原始netCDF文件，只需要在匹配逻辑中过滤掉全为NA的陆面网格即可，同时还修复了原代码中data_u的笔误问题：

# 原读取文件、经纬度、时间、shts变量的逻辑保持不变
nc <- nc_open("swell_bilmap.nc")
lat <- ncvar_get(nc,'lat')
lon <- ncvar_get(nc,'lon')
t <- ncvar_get(nc, "time")
timestamp <- as_datetime(c(t*60*60),origin="1900-01-01")
data_shts <- ncvar_get(nc,"shts")
nc_close(nc)

# 新增：筛选有效海洋网格，判断每个网格的shts值是否不全为NA
valid_mask <- apply(data_shts, c(1,2), function(x) !all(is.na(x)))
lonlat <- expand.grid(lon=lon, lat=lat)
lonlat$valid <- as.vector(valid_mask)
# 仅保留有效网格
lonlat_valid <- lonlat[lonlat$valid, ]

# 基于有效网格生成空间对象
coord_brookings <- st_as_sf(lonlat_valid, coords=c("lon","lat")) %>%
  st_set_crs(4326)

# 港口点生成逻辑保持不变
psj_brookings <- st_point(c(PORTS[1,5],PORTS[1,4])) %>%
  st_sfc() %>%
  st_set_crs(4326)

# 计算距离后取最近的有效网格
coord_brookings$dist <- st_distance(coord_brookings, psj_brookings)
nearest_grid <- coord_brookings[which.min(coord_brookings$dist), ]

# 提取最近网格对应的时间序列
nearest_lon <- st_coordinates(nearest_grid)[1]
nearest_lat <- st_coordinates(nearest_grid)[2]
lon_idx <- which(lon == nearest_lon)
lat_idx <- which(lat == nearest_lat)
df_brookings_shts <- data.frame(
  shts = data_shts[lon_idx, lat_idx, ],
  time = timestamp
)

CDO实现方案 #

如果需要提前生成剔除陆面网格的netCDF文件，可执行以下命令：

# 生成有效海洋网格掩码，shts不全为NA的网格标记为1，其余为0
cdo -expr,'mask = !all(is.na(shts))' swell_bilmap.nc mask.nc
# 基于掩码过滤原始文件，输出仅含有效海洋网格的新文件
cdo ifthen mask.nc swell_bilmap.nc swell_bilmap_valid.nc

生成的swell_bilmap_valid.nc可直接用于后续匹配分析。

ArcGIS Pro实现方案 #

1. 用多维工具集的创建NetCDF要素图层工具导入原始netCDF文件，变量选择shts，坐标对应lon、lat生成要素图层
2. 用按属性选择工具，筛选条件设置为shts IS NOT NULL，选中所有有效海洋网格
3. 右键图层选择「数据>导出要素」保存为新的要素类，后续直接用该要素类和港口点做近邻分析即可获取最近有效网格。

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