解决NDVI线性回归等值线图海洋填充问题（Xarray/Numpy/Matplotlib）

云原生阿岚

2026-7-5

解决NDVI线性回归等值线图海洋区域错误填充问题

核心问题分析

你遇到的问题本质是线性回归计算时丢失了原始海洋区域的NaN掩码：CDO掩膜后的海洋区域本应为全NaN，但np.polyfit处理含NaN的数组时，若直接对整组数据计算（未过滤NaN），会导致无值区域被错误填充为0或其他数值。关键解决方案是保留原始掩码，在回归计算前后严格维护无值区的NaN状态。

具体解决步骤&代码示例

1. 读取数据并提取原始海洋掩码

先读取CDO处理后的NetCDF数据，提取所有时间步均为NaN的区域（即海洋掩膜区），后续用于恢复无值状态：

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# 读取掩膜后的NDVI数据
ds = xr.open_dataset('your_ndvi_data.nc')
ndvi_data = ds['ndvi'].values  # 假设变量名为ndvi，维度为(时间, 纬度, 经度)
time_axis = np.arange(len(ds.time))  # 生成时间序列索引（如年数）
lat = ds.lat.values
lon = ds.lon.values

# 提取原始海洋掩码：所有时间步均为NaN的网格点
ocean_mask = np.isnan(ndvi_data).all(axis=0)

2. 矢量化计算线性回归斜率（保留NaN）

避免逐网格循环的低效，用矢量化方法对每个网格点的有效时间序列做拟合，直接跳过NaN：

# 将三维数据展平为(时间, 空间点)的二维数组
ndvi_flat = ndvi_data.reshape(ndvi_data.shape[0], -1)

# 定义单个空间点的斜率拟合函数
def calc_trend(ts):
    valid_idx = ~np.isnan(ts)
    # 至少需要2个有效点才能拟合，否则返回NaN
    if valid_idx.sum() < 2:
        return np.nan
    return np.polyfit(time_axis[valid_idx], ts[valid_idx], 1)[0]

# 对所有空间点批量计算斜率
trend_flat = np.apply_along_axis(calc_trend, 0, ndvi_flat)
# 恢复为(纬度, 经度)的二维数组
ndvi_trend = trend_flat.reshape(ndvi_data.shape[1:])

# 强制恢复海洋区域的NaN状态
ndvi_trend[ocean_mask] = np.nan

3. 绘制等值线图（自动忽略NaN区域）

用contourf绘图时，NaN区域会自动留白，无需额外处理：

plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines(linewidth=0.8)

# 绘制等值线图，extend参数控制颜色条的极值显示
contour = ax.contourf(lon, lat, ndvi_trend, 
                      cmap='RdYlGn', levels=15, 
                      extend='both', transform=ccrs.PlateCarree())

# 添加颜色条与坐标轴格式化
plt.colorbar(contour, orientation='horizontal', pad=0.05, 
             label='NDVI Linear Trend (per year)')
ax.set_xticks(np.arange(-180, 181, 60), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(ccrs.mpl.ticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(ccrs.mpl.ticker.LatitudeFormatter())

plt.title('NDVI Annual Linear Trend')
plt.show()

关键注意事项

不要用np.nan_to_num()将NaN转为0：这会直接破坏原始掩膜信息，是导致海洋着色的根本原因。
优先用原始数据掩码而非maskoceans：CDO已完成针对性掩膜，用原始全NaN区域的掩码比基于地形数据的maskoceans更精准，避免额外的误掩/漏掩。
拟合时必须检查有效点数量：若某个网格点有效数据不足2个，直接返回NaN，避免无效拟合结果。

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