nc数据补值
我的NC数据的右上角缺失值，准备使用IDW插值进行补充。上图为补值前，下图为补值后：
整体逻辑为：
代码：
import numpy as np
from netCDF4 import Dataset
from scipy.spatial import cKDTree
import os
def fill_missing_idw(variable, lat, lon, power=2, max_neighbors=8):
"""
使用IDW方法填补二维数组中的缺失值。
参数:
- variable: 2D numpy数组，可能包含缺失值 (masked array)
- lat, lon: 2D numpy数组，与variable形状相同的纬度和经度
- power: IDW的幂参数
- max_neighbors: 使用的最大邻居数量
返回:
- filled_variable: 填补后的二维numpy数组
"""
filled_variable = variable.copy()
ny, nx = variable.shape
if isinstance(variable, np.ma.MaskedArray):
known_mask = ~variable.mask
missing_mask = variable.mask
known_values = variable.data[known_mask]
else:
known_mask = ~np.isnan(variable)
missing_mask = np.isnan(variable)
known_values = variable[known_mask]
known_points = np.column_stack((lon[known_mask], lat[known_mask]))
missing_points = np.column_stack((lon[missing_mask], lat[missing_mask]))
if len(known_points) == 0:
print("没有可用的已知点进行插值。")
return filled_variable
tree = cKDTree(known_points)
distances, indexes = tree.query(missing_points, k=max_neighbors, n_jobs=-1)
if max_neighbors == 1:
distances = distances[:, np.newaxis]
indexes = indexes[:, np.newaxis]
with np.errstate(divide='ignore'):
weights = 1 / distances**power
weights[~np.isfinite(weights)] = 0
interpolated_values = np.sum(weights * known_values[indexes], axis=1) / np.sum(weights, axis=1)
exact_matches = distances[:,0] == 0
interpolated_values[exact_matches] = known_values[indexes[exact_matches, 0]]
filled_variable[missing_mask] = interpolated_values
# 调试信息
n_filled = np.sum(~filled_variable[missing_mask].mask) if isinstance(filled_variable, np.ma.MaskedArray) else np.sum(~np.isnan(filled_variable[missing_mask]))
print(f"尝试填补了 {len(missing_points)} 个缺失点，成功填补了 {n_filled} 个。")
return filled_variable
def fill_missing_idw(variable, lat, lon, power=2, max_neighbors=8):
"""
使用IDW方法填补二维数组中的缺失值。
参数:
- variable: 2D numpy数组，可能包含缺失值 (masked array)
- lat, lon: 2D numpy数组，与variable形状相同的纬度和经度
- power: IDW的幂参数
- max_neighbors: 使用的最大邻居数量
返回:
- filled_variable: 填补后的二维numpy数组
"""
filled_variable = variable.copy()
ny, nx = variable.shape
if isinstance(variable, np.ma.MaskedArray):
known_mask = ~variable.mask
missing_mask = variable.mask
known_values = variable.data[known_mask]
else:
known_mask = ~np.isnan(variable)
missing_mask = np.isnan(variable)
known_values = variable[known_mask]
known_points = np.column_stack((lon[known_mask], lat[known_mask]))
missing_points = np.column_stack((lon[missing_mask], lat[missing_mask]))
if len(known_points) == 0:
print("没有可用的已知点进行插值。")
return filled_variable
tree = cKDTree(known_points)
distances, indexes = tree.query(missing_points, k=max_neighbors)
if max_neighbors == 1:
distances = distances[:, np.newaxis]
indexes = indexes[:, np.newaxis]
with np.errstate(divide='ignore'):
weights = 1 / distances**power
weights[~np.isfinite(weights)] = 0 # 处理距离为0的情况
# 计算权重和
