我最早开始用Python处理气象海洋数据的时候和多数人一样先扑向数据可视化——但很快就被打蒙了。contourf出图倒是快可海岸线全歪了纬度方向也反了画风场又找不到合适的投影方式。折腾几天后我才意识到气象海洋数据与普通表格数据完全不是一回事它的可视化背后绑定着一整套数据结构和投影逻辑而数据分析方法也有自己的一套“约定俗成”。这篇内容就围绕python气象海洋这两条主线展开从数据读取预处理到常用可视化类型再到趋势分析、EOF、相关分析、谱分析等方法最后把那些不踩一遍很难记住的坑集中整理出来。适合刚入坑的学生、转行做模式后处理的科研人员以及想在气象海洋方向把Python用得更顺手的朋友。1. 气象海洋数据的“脾气”四维网格、NetCDF与生态选型1.1 为什么气象海洋数据不能用普通表格的思路处理日常数据分析中我们面对的往往是二维表格行是样本列是特征。但气象海洋领域的原始数据通常是四维的——时间time、垂直层level或depth、纬度latitude、经度longitude。比如一个再分析资料的逐日海表温度场它是“时间×纬度×经度”的三维数组如果再叠加多个深度层就变成四维。用pandas处理这种数据不是不行但会把好好的网格结构拍扁丢失坐标信息。比如你知道某一列是温度但不知道它对应哪个经纬度你做区域平均时还得自己维护一大串索引。更麻烦的是海洋气象领域广泛使用的NetCDFNetwork Common Data Form和GRIB格式内部除了数值外还带大量元数据——单位、缺测值、坐标属性、数据来源说明等。这些信息如果被拍平成表格基本就丢了。所以气象海洋Python分析的第一个关键认知就是不要用二维思维硬啃多维数据。正确做法是用xarray这一层“数据容器”把维度、坐标、变量组织起来再基于它做切片、聚合和可视化。1.2 NetCDF/GRIB格式的基础概念NetCDF是气象海洋领域最通用的自描述格式之一本质是一个多维数组容器包含几个核心概念dimensions维度比如lat、lon、time、level用大小描述。variables变量比如温度temperature、风场U分量u、海表温度sst每个变量都有它的维度。attributes属性用来描述变量信息比如单位units、缺测值_FillValue、长名称long_name等。GRIB格式在数值预报产品中非常常见它把数据按“二进制段”存储每一段都带有描述要素、层次、时效的编码。在Python里处理GRIB我一般用cfgrib配合xarray或者直接用pygrib读出来之后转换成numpy数组再做后续操作。如果你打开一个典型的SST NetCDF文件用xarray.open_dataset检查时结构大概是这样的import xarray as xr ds xr.open_dataset(sst_monthly.nc) print(ds)输出里会清晰显示各个变量的维度、坐标、属性和缺测标记比直接看二进制文件“可读性”高出一个数量级。1.3 构建一个气象海洋分析的基础Python环境我自己的常用环境组合是python3.10左右即可xarraynumpyscipymatplotlibcartopynetCDF4或h5netcdfpandas。如果涉及模式后处理还会加cfgrib、windspharm做EOF时可以用eofs或xeofs但多数情况下自己用numpy.linalg.svd写也就够了。安装的时候强烈建议直接用conda从conda-forge通道装特别是cartopy这种带底层地理数据的库用pip装有时候会缺GEOS和PROJ运行库导致运行时反复报错。一条命令到位conda create -n atm python3.10 -c conda-forge conda activate atm conda install -c conda-forge xarray numpy scipy matplotlib cartopy netCDF4 pandas提示conda install -c conda-forge cartopy会自动解决GEOS/PROJ依赖比手动编译省心太多。这是我在环境配置上最大的一个省时经验。依赖装好之后import cartopy.crs as ccrs能顺利通过就说明地图投影组件已经能用了。2. 用xarray把多维数据“拎起来”从读取到预处理的典型操作2.1 DataArray与Dataset按名字而不是按位置取数据xarray里有两种核心对象Dataset和DataArray。Dataset类似于多个变量打包在一起的容器而DataArray是单个变量的多维数组。日常操作中我会先用Dataset打开整个文件然后通过变量名取出DataArray来画图或计算sst ds[sst] # 按坐标选区域和时段而不是手工数索引 sst_region sst.sel(latslice(-10, 10), lonslice(110, 160), timeslice(2020-01-01, 2020-12-31))这里sel就是按坐标值“标签化”索引。再配合.isel按位置索引比如取第0层u_500 ds[u].isel(level0)如果你更习惯数值下标也可以直接用sst_point sst.sel(lat20, lon130, methodnearest)methodnearest会自动寻找最近的格点省去自己写距离计算的功夫。这在提取站点附近格点值的时候尤其好用。2.2 缺失值、_FillValue与数据清洗气象海洋数据里出现缺测太常见了。NetCDF文件一般用_FillValue指定一个特殊数值来表示缺测读入xarray后这个值通常会被自动转成NaN。但也有例外尤其是地形遮盖的区域、模式陆地点、某些同化产品中的缺测层。如果你发现画出的图上有奇怪的“值域”先跑一下import numpy as np print(np.nanmin(sst.values), np.nanmax(sst.values))如果出现了类似-32767这种极端值就说明_FillValue仍然以真实数值形式存在需要手动处理sst_clean sst.where(sst -100)或者显式设置sst_clean xr.where(sst -32767, np.nan, sst)注意在计算气候态、距平之前一定要先统一缺测处理方式。否则后面mean()、std()算出来的全是“污染”对象排查起来还极其隐蔽。另外xarray默认在计算时会跳过NaN吗不会全跳过比如ds.mean(dimtime)会把NaN继续传播除非你写skipnaTrue。好在xarray对大多数reduce操作默认skipnaTrue但如果自行调numpy的np.mean(data)就会被NaN坑一把。建议在做任何统计前先np.isfinite()筛选或利用skipna特性。2.3 经纬度顺序、坐标名称与标准化同一个变量在不同数据产品里坐标名字和方向是不一样的。有的用latitude/longitude有的用lat/lon甚至有的用y/x纬度可能是从小到大排也可能是从大到小排。面对这种问题我一般做两步统一坐标名称ds ds.rename({latitude: lat, longitude: lon})让纬度严格递减配合绘图习惯ds ds.sortby(lat)经度方面有些产品经度范围是0到360有些是-180到180。如果遇到不一致用这个方式统一ds ds.assign_coords(lon(((ds.lon 180) % 360) - 180)).sortby(lon)这一步的目的不是“洁癖”而是避免后续做区域切片或插值时因为边界问题把太平洋切成两半。3. 可视化第一步Cartopy投影、填色图与风场图的“标准动作”3.1 地图投影先有坐标系再谈画图气象海洋可视化绕不开地图。直接拿matplotlib的contourf画经纬度网格出来的图是完全失真的——高纬度地区极度拉宽海岸线也没有。所以需要cartopy来定义“数据坐标”与“显示投影”。一个标准的全球海表温度填色图框架如下import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.coastlines(linewidth0.5) ax.add_feature(cfeature.LAND, colorlightgray) im ax.contourf(sst.lon, sst.lat, sst.mean(dimtime), levels24, cmapRdYlBu_r, transformccrs.PlateCarree()) plt.colorbar(im, axax, shrink0.7, labelSea Surface Temperature (°C)) ax.set_global() plt.title(Global SST Climatology) plt.show()这里最关键的是transformccrs.PlateCarree()。它的意思是数据本身是经纬度坐标等距圆柱投影而你显示的projection可以是任意投影。如果你不写transformcartopy会默认数据坐标也等于投影坐标在非圆柱投影下数据会错位得离谱。我在实际工作中的习惯是先用PlateCarree快速检查数据形态再用Orthographic或Robinson出汇报图。3.2 配色、层级和色标让图“可信”而不是“好看”在气象海洋圈子里matplotlib默认的viridis适合作业但专业配色我会更倾向于温度场RdYlBu_r红蓝翻转暖色代表高值冷色代表低值降水场BrBG或者cmocean里的turbid、rain风场用cmocean的balance或者直接matplotlib的PuOr距平场统一红蓝色系居中关键点是把vmin和vmax设置成对称的vmax 2.0 clevs np.linspace(-vmax, vmax, 21) im ax.contourf(lon, lat, anomaly, levelsclevs, cmapRdBu_r, transformccrs.PlateCarree())这是很实用的建议画距平场时颜色bar中心必须是白色0值处否则别人一眼觉得你故意夸大正异常。等值线叠加时用contour加clabelcs ax.contour(lon, lat, field, levelsclevs[::2], colorsk, linewidths0.6, transformccrs.PlateCarree()) ax.clabel(cs, cs.levels, inlineTrue, fontsize8, fmt%.1f)这样既有填色的大局观又有等值线的数值精度适合在论文级图中使用。3.3 风场与流线矢量图的画法细节风场u和v分量的可视化常见的有三种方式箭头图quiver、风羽图barb、流线图streamplot。quiver的空间间隔是重点。直接画全网格图面会密密麻麻挤成一团。通常我会在sel出目标层次后每隔3~5个格点画一个箭头import numpy as np # 假设u10m, v10m为10m风场DataArray lon u10m.lon.values lat u10m.lat.values step 4 ax.quiver(lon[::step], lat[::step], u10m.values[::step, ::step], v10m.values[::step, ::step], transformccrs.PlateCarree(), scale10, width0.002)风羽图更适合给“站址/格点风”做专业呈现ax.barbs(lon[::step], lat[::step], u10m.values[::step, ::step], v10m.values[::step, ::step], transformccrs.PlateCarree(), length5, sizesdict(emptybarb0.2))流线图则更适合展示气流“轨迹感”ax.streamplot(lon, lat, u10m, v10m, transformccrs.PlateCarree(), density1, colork, linewidth0.8)但流线图有两个坑一是遇到奇异点风速为0时会画得杂乱二是区域较大时计算量猛增。所以流线图我一般只在中小区域用比如台风路径或副高外围。另外画台风或涡旋时叠加contourf风速大小sqrt(u^2v^2)配上quiver方向比单纯画矢量更直观。4. 进阶图谱剖面图与Hovmöller图的绘制思路拓展4.1 垂直剖面图把三维结构“切”开看气象海洋研究经常要回答“某一要素随高度或深度怎么变化”。比如看海洋温度沿经向的垂直剖面或者大气位势高度沿纬向的垂直剖面。方法很直白——把目标方向留作横轴垂直层作为纵轴# 选某一条经线上的温度剖面 temp_cross ds[temperature].sel(lat20, lonslice(110, 160), methodnearest) # 经度-深度二维场 fig, ax plt.subplots(figsize(10, 5)) cf ax.contourf(temp_cross.lon, temp_cross.depth, temp_cross.mean(dimtime), levels20, cmapRdYlBu_r) ax.set_ylim(1000, 0) # 深度轴转成向下 plt.colorbar(cf, labelTemperature (°C))这里最容易被忽略的就是set_ylim反转深度轴——因为海洋/大气的垂直坐标习惯上是“上小下大”或者“等压面随高度递减”你不翻转坐标轴读者会看得非常别扭。如果使用的是等压面数据比如1000hPa到100hPa同理可用ax.invert_yaxis()。4.2 Hovmöller图时间-经度纬度传播图的画法Hovmöller图在气象海洋里用来追踪信号传播比如赤道开尔文波、季节内振荡MJO、Rossby波列。它把一维空间固定纬度通常是赤道附近或5°N-5°S平均作为横轴时间作为纵轴填色展示要素演变# 先做区域平均 eq_sst sst.sel(latslice(-5, 5)).mean(dimlat) # 再做季节循环剔除后面会展开 eq_sst_anom eq_sst.groupby(time.month) - eq_sst.groupby(time.month).mean(time) fig, ax plt.subplots(figsize(9, 8)) im ax.contourf(eq_sst_anom.lon, eq_sst_anom.time, eq_sst_anom, levelsnp.linspace(-1, 1, 21), cmapRdBu_r, extendboth) ax.set_xlabel(Longitude) ax.set_ylabel(Time) plt.title(Time-Longitude Section of SST Anomaly (5°S–5°N)) plt.colorbar(im, labelSST anomaly (°C))画这种图我一般还会叠加一条“东传/西传”的参考速度线比如斜线代表信号的相速度方便读者判断传播方向。实现上就是ax.plot(lon, time_ref, k--)这类不复杂但很出效果。4.3 让多子图布局和公共色标更省事当需要比较多个时段、多个层次或多个区域时子图间保持“同一个色标范围”是硬要求。我用过最顺的方案是使用mpl_toolkits.axes_grid1或者GridSpec加cbar_aximport matplotlib.gridspec as gridspec fig plt.figure(figsize(12, 5)) gs gridspec.GridSpec(1, 2, width_ratios[1, 1], wspace0.1) ax1 fig.add_subplot(gs[0], projectionccrs.PlateCarree()) ax2 fig.add_subplot(gs[1], projectionccrs.PlateCarree()) vmin, vmax, cmap -2, 2, RdBu_r cf1 ax1.contourf(lon, lat, data1, levelsnp.linspace(vmin, vmax, 21), cmapcmap, transformccrs.PlateCarree()) cf2 ax2.contourf(lon, lat, data2, levelsnp.linspace(vmin, vmax, 21), cmapcmap, transformccrs.PlateCarree()) cbar_ax fig.add_axes([0.92, 0.15, 0.02, 0.7]) fig.colorbar(cf2, caxcbar_ax, labelAnomaly)上面这步比自动缩放色标的关键点在于它让两个子图共享levels和vmin/vmax不会出现“同一种颜色在两个图里代表不同值”的误导。5. 数据分析方法趋势、EOF、相关与谱分析中的常用套路5.1 趋势分析与显著性检验气象海洋数据分析最基础的需求就是看“某个量随时间有没有显著变化”。比如全球平均海表温度是否在升高答案是算线性趋势。具体做法是先把原始三维/四维数据在某个区域做平均得到一维时间序列然后用scipy.stats.linregress计算斜率和p值from scipy import stats series sst_region.mean(dim[lat, lon]) years series.time.dt.year.values.astype(float) slope, intercept, r_value, p_value, std_err stats.linregress(years, series.values) # 每十年变化 slope * 10 trend_per_decade slope * 10如果要在二维空间场上逐个格点计算趋势可以用xarray的polyfit或直接循环格点。但循环几千个格点会很慢。我建议用xr.polyfit(dimtime, deg1)trend ds[sst].polyfit(dimtime, deg1) sst_trend trend.polyfit_coefficients.isel(degree0) * 10注意自由度问题如果数据是月平均序列相邻月份存在自相关有效样本量低于实际月数。严谨一点的做法是用“有效样本量”修正t检验否则p值会偏乐观。简单处理可以先对月均序列做12个月滑动平均降噪再重新采样成年平均减小自相关。5.2 EOF分析手写一套更踏实经验正交函数分解EOF也叫PCA在气象海洋里太常用了——用来提取空间主导模态及其时间演变。比如北大西洋涛动NAO、El Niño的空间模态基本就是EOF第一模态这种角色。虽然eofs和xeofs库能直接算但我建议自己写一次核心过程这样你对“特征向量、时间系数、方差贡献率”的理解会扎实很多import numpy as np import xarray as xr # data形状: (time, lat, lon)先加权并去除缺测 def eof_2d(data): # 1. 数据矩阵化 (time, space) nt, ny, nx data.shape X data.reshape(nt, -1) # 2. 去掉空间点全为缺测的列 valid np.isfinite(X).all(axis0) X_valid X[:, valid] # 3. 距平 X_anom X_valid - X_valid.mean(axis0, keepdimsTrue) # 4. 对空间协方差阵做SVD分解 # X U * S * VtVt的行就是空间模态 U, S, Vt np.linalg.svd(X_anom, full_matricesFalse) # 时间系数 U * S (或直接投影) pc U * S # 空间模态EOF eof Vt # 方差贡献 variance S**2 / (X_anom.shape[0] - 1) explained variance / variance.sum() * 100 return pc, eof, explained拿到eof后可以把它重新reshape成(neof, ny, nx)的空间模态同时把time坐标绑定在pc上。在实际使用中我通常还会给数据做“纬度加权”——因为等经纬度网格在高纬的格点面积偏小直接用原始网格做EOF会过度放大高纬度信息。常用做法是乘以sqrt(cos(lat))weights np.sqrt(np.cos(np.deg2rad(lat))) data_weighted data * weights # 画图还原时除以权重这一步在分析全球尺度的海温、海平面气压场时尤其重要无数新手在这个细节上栽过跟头。5.3 相关分析与合成分析两种常用“找联系”的手段相关分析用来回答“两个要素之间是否存在线性关联”。最常见的场景是做一个指数比如Nino3.4指数与全球降水场的相关分布找到遥相关型。逐格点的Pearson相关可以这样算corr xr.corr(index_series, precip, dimtime)但这里有个问题随便拿一个指数和全球所有格点做相关样本空间巨大显著性检验必须做。最简单的方式是转成t统计量t corr * np.sqrt((n - 2) / (1 - corr**2)) p 2 * (1 - stats.t.cdf(np.abs(t), dfn-2))然后把p 0.05的区域打点标注出来。合成分析则是根据某一个离散事件比如台风年、El Niño年、热浪年来分组做平均看两个状态下的差异是否显著。做法是# 按指数年份分组 warm_years years_index[nino34 0.5] cold_years years_index[nino34 -0.5] # 取对应年份的要素场 warm_composite sst_annual.sel(timewarm_years).mean(dimtime) cold_composite sst_annual.sel(timecold_years).mean(dimtime) diff warm_composite - cold_composite之后可以用t检验对差异场做显著性检验再把显著区域用stippling打点标出。5.4 功率谱与滤波从时间序列里找周期气象海洋里有很多周期现象季节内30-60天、年循环、ENSO2-7年、准两年振荡QBO等。用功率谱估计可以发现序列中的显著周期。scipy.signal里提供了welch和periodogram。对于逐日/逐月数据我常用Welch方法因为它对噪声更稳健from scipy import signal f, Pxx signal.welch(series, fs1.0, nperseglen(series)//4) # 转换成周期 period 1 / f画图时横轴用周期单位月或者天纵轴用功率谱密度再叠加红噪声背景和显著性水平线。对于季节内振荡这类信号滤波是常用前置手段。例如做一个30-60天带通滤波from scipy.signal import butter, filtfilt def bandpass_filter(data, dt, low, high, order5): nyq 0.5 / dt lowb low / nyq highb high / nyq b, a butter(order, [lowb, highb], btypeband) return filtfilt(b, a, data, axis-1) # dt 1若数据为逐日 filtered bandpass_filter(series.values, dt1, low1/60, high1/30)这类滤波器在分析MJO、季节内振荡时几乎是标配。加上功率谱就能很清晰地判断某个频段上是否存在显著信号。6. 踩坑笔记气象海洋Python高频翻车点与排查思路6.1 纬度从北往南与从南往北图为什么“倒”了不少全球数据尤其海洋模式输出的纬度是从-90到90也就是从南极到北极排的。如果你不做sortby(lat)直接contourf图面会上下颠倒——赤道在下面南极在上面看图的人一头雾水。我的排查思路很简单先打印ds.lat.values[:5]和ds.lat.values[-5:]看首尾值。如果首尾不是北高南低就执行ds ds.sortby(lat)。另一个容易踩的是pcolormesh和contourf对二维纬度数组比如曲线网格的处理方式不同。如果遇到非矩形网格如MPAS、WRF的扭曲网格还直接用contourf(lon, lat, data)会失败或变形这时候要么插值到规则网格要么用pcolormesh加shadingauto。6.2 投影与数据坐标不一致最隐蔽的地图错位这是一个让我调了一整晚的坑用ccrs.NorthPolarStereo()投影画北极海冰忘记加transformccrs.PlateCarree()出来的海冰范围偏移到了大西洋中部。记住一个原则只要你的数据是经纬度网格transform必须是ccrs.PlateCarree()而projection可以是任何投影。两者负责的事情不一样前者告诉cartopy“数据从哪里来”后者告诉它“显示到什么地方”。6.3 时间维解码cftime与日期索引问题NetCDF里的时间一般是以“天”或“小时”为单位的数值配合units和calendar属性存储。xarray默认会尝试解码成datetime64。但如果遇到360_day日历或noleap日历常见于气候模式输出xarray就会使用cftime类型。此时sel(time2020-01-01)可能直接报错。解决方案是明确指定解码方式ds xr.open_dataset(model_output.nc, decode_timesTrue, use_cftimeTrue)或者干脆不解码时间自己手动处理ds xr.open_dataset(model_output.nc, decode_timesFalse) # 读取time变量做转换如果你在用cftime对象做季节循环、年际合成时先转成数值年份time.dt.year不一定有效可以手动year np.array([t.year for t in ds.time.values])6.4 内存不足与性能问题分块、降采样与延迟计算再分析资料动辄几个GB。如果一次性open_dataset后再做复杂运算内存很容易爆掉。我的处理思路有两个一是用open_dataset(..., chunks{time: 10})结合dask让计算延迟到需要时才执行并在多核上并行ds xr.open_dataset(large_sst.nc, chunks{time: 10})二是先做降采样或者只载入需要的区域sst_region ds[sst].sel(latslice(-30, 30), lonslice(90, 180)).load()只在最后出图/算结果需要全部数据时才调.load()否则一直保持延迟计算。如果循环处理多个文件另一种提速思路是提前把多文件合并成一个数据集ds xr.open_mfdataset(data_*.nc, combineby_coords)这时候chunks也会自动分配效率和手写循环完全不同。6.5colorbar范围被自动缩放对比图“失真”问题你在单张图上可能注意不到但一旦做前后对比、多情景对比contourf默认会按每个子图自己的最大最小值分配颜色层级导致“同样蓝色在一张图代表0.5°C在另一张图代表2°C”。这是气象海洋论文审稿人最常指出的一点。解决办法就是显式传入vmin、vmax和levels让所有子图共用一个色标范围。前面已经写过代码这里再强调一下在多子图设置中宁可手动指定范围也不要依赖自动scale。写在最后的一点个人建议如果让我给刚接触这个方向的人一个建议我的想法是与其零散地学Python语法不如从复现一张经典图开始。找一份ERSST或ERA5数据先画全球海表温度气候态再做区域平均、季节循环、距平场接着叠加EOF、趋势和显著性检验。这个流程走下来xarray的常用操作、cartopy投影设置、统计检验的基本套路基本都能覆盖到。遇到报错时最常见的不是语法问题而是“坐标系没对上”“时间维没解码”“数据里有掩膜值”这类数据层面的坑。所以排查顺序也建议从数据本身开始先print(ds)看维度与坐标再画一维时间序列确认量纲量级最后才做2D填色。这样能避免大量无效调试。数据可视化只是手段数据分析方法才是核心支撑。气象海洋数据不是普通表格它有时间、空间、垂直结构的完整性只有把数据容器、投影体系、统计工具打通才能把一个研究问题从原始数据推进到物理机制。希望这些实操经验能帮你在处理气象海洋数据时少走弯路。