含NaN值的网格插值绘图出现断裂，寻求数据平滑方法

从大型数据集中生成透视表，提取y（压力）、x（纬度）、z（AOU值）构建网格后执行插值，但绘制的等值线图在纬度78-80附近出现数据断裂点，推测行方向插值未正确执行，需要有效的数据平滑方法。

原始代码 #

import numpy as np
import matplotlib.pyplot as plt

aou_df = df_1994.pivot_table(index='CTDPRS', columns = 'LATITUDE', values='AOU')
aou_df = aou_df.interpolate(method='linear', limit_area='inside', axis =0 )

##Plotting AOU 1994
y = [4.2, 4.7, 4.8, 4.9, 5.4, 9.1, 9.6, 9.7, 10.0, 10.1,
     ...
     3568.2, 3608.6, 3818.6, 3824.9, 3866.7, 3979.1, 4013.4, 4133.1,
     4159.3, 4287.3]

x = [72.13, 73.0, 73.49, 73.98, 74.5, 75.0, 75.45, 75.75, 75.94,
     76.62, 77.33, 77.78, 78.14, 78.15, 78.98, 79.98, 80.15, 80.16,
     80.33, 80.71, 81.24, 81.58, 82.47, 83.17, 84.06, 84.85, 85.89,
     87.16, 88.06, 88.79, 88.86, 88.95, 89.02, 90.0]

z = [[-12.29372749, np.nan, np.nan, ..., np.nan, np.nan, np.nan],
     [np.nan, np.nan, -43.41465869, ..., np.nan, np.nan, np.nan],
     [np.nan, -54.49999783, np.nan, ..., np.nan, np.nan, np.nan],
     ...,
     [np.nan, np.nan, np.nan, ..., np.nan, np.nan, np.nan],
     [np.nan, np.nan, np.nan, ..., np.nan, 55.87256821, np.nan],
     [np.nan, np.nan, np.nan, ..., 55.39665852, np.nan, 55.05005376]]

xi, yi = np.meshgrid(x,y,indexing='ij')

plt.figure(figsize=(25,10))
levels = np.linspace(-135,135)

plt.contourf(xi,yi,z, cmap = 'jet', levels=levels,vmin=-135, vmax=135)
plt.gca().invert_yaxis()
plt.gca().invert_xaxis()
cbar = plt.colorbar(ticks=(-135,-110,-85,-60,-35,0,35,60,85,110,135), extend= 'both')
cbar.set_label('AOU', fontsize=18)
cbar.ax.tick_params(labelsize=18)
plt.xlabel('LAT',fontsize=18)
plt.ylabel('Pressure (dbar)' ,fontsize=18)
plt.ylim(bottom = 1000)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.plot(x,range(len(x)),'gD', clip_on=False, markersize=10)

问题可视化 #

解决方案 #

1. 修正插值顺序：先纬度方向再压力方向 #

当前仅对压力轴（axis=0）插值，但纬度轴（axis=1）的稀疏缺失值是断裂核心。调整插值顺序，先填充纬度方向的缺失：

aou_df = df_1994.pivot_table(index='CTDPRS', columns='LATITUDE', values='AOU')
# 先对纬度列方向插值，限制在数据有效范围内
aou_df = aou_df.interpolate(method='linear', limit_area='inside', axis=1)
# 再对压力行方向插值
aou_df = aou_df.interpolate(method='linear', limit_area='inside', axis=0)
# 提取z值
z = aou_df.values

2. 用径向基函数（RBF）做全局插值 #

线性插值对稀疏数据效果有限，RBF适合海洋这类不规则采样数据，能更好地拟合空间分布：

from scipy.interpolate import Rbf

# 提取所有非NaN的原始采样点
mask = ~np.isnan(z)
x_points = np.tile(x, len(y))[mask.flatten()]
y_points = np.repeat(y, len(x))[mask.flatten()]
z_points = z[mask]

# 创建RBF插值模型，multiquadric适合平滑曲面
rbf = Rbf(x_points, y_points, z_points, function='multiquadric')
# 在现有网格上生成插值结果
z_interp = rbf(xi, yi)

# 替换原z值绘图即可
plt.contourf(xi, yi, z_interp, cmap='jet', levels=levels, vmin=-135, vmax=135)

3. 生成密集网格+三次插值 #

78-80纬度区间采样点稀疏，生成更密集的纬度网格后用三次插值填充：

from scipy.interpolate import griddata

# 生成密集纬度网格（72到90，步长0.1）
xi_dense = np.linspace(72, 90, 181)
yi_dense = y
xi_grid, yi_grid = np.meshgrid(xi_dense, yi_dense, indexing='ij')

# 整理原始采样点
points = []
values = []
for i in range(len(x)):
    for j in range(len(y)):
        if not np.isnan(z[j][i]):
            points.append((x[i], y[j]))
            values.append(z[j][i])

# 三次插值填充密集网格
z_dense = griddata(points, values, (xi_grid, yi_grid), method='cubic')

# 用密集网格绘图
plt.contourf(xi_grid, yi_grid, z_dense, cmap='jet', levels=levels, vmin=-135, vmax=135)

4. 后处理高斯平滑 #

如果插值后仍有局部断裂，用高斯滤波做后处理平滑：

from scipy.ndimage import gaussian_filter

# sigma控制平滑强度，1.5左右兼顾平滑和细节
z_smoothed = gaussian_filter(z_interp, sigma=1.5)

plt.contourf(xi, yi, z_smoothed, cmap='jet', levels=levels, vmin=-135, vmax=135)

内容的提问来源于stack exchange，提问作者Maria Cristina Alvarez