气象统计实习Python代码全解析：从数据读取到EOF制图

简介：面向气象数据统计分析学习者的Python实习资源包，聚焦利用Python完成气象数据的读取、清洗、统计计算与可视化，适合高校气象、环境相关专业学生及初步接触气象编程的研究者。压缩包共12个文件，包括8个Python脚本（对应实习1至实习8的完整代码）、2个dat格式气象数据文件、1个ctl描述文件及1个txt数据说明，整体仅80KB，轻量而条理清晰。已有1570人学习下载。脚本覆盖Pandas数据处理、NumPy数值计算、Matplotlib/Seaborn绘图、时间序列分析等核心知识点，并涉及平均值、标准差、相关与回归等统计方法的具体实现；配套的H500.DAT、MA.DAT等样例数据与描述文件，可帮助快速复现实验结果。通过研读源码、对照运行，能够系统掌握气象统计方法从数据预处理到结果可视化的完整流程，有效节省课程实习中的调试时间，也适合作为自学气象Python编程的入门参考。

1. 先算清楚再画图：这份实习代码包到底解决什么

拿到气象统计方法实习的 Python 代码包，我建议先别急着跑，先想清楚这份资源替你解决了哪一步。气象统计实习卡人的往往不是统计公式本身，而是数据从读入到结果解释的整条链路：NetCDF 怎么读、气候态用哪一段作基准、距平正负怎么合理解释、EOF 算出来怎么判断模态真假，最后还得画出一张能放进报告里的图。这份 Python 版实习资源，正好把这条链路按可复现的脚本串起来了。适合三类人：正在写实习报告的气象与地学本科生、需要快速验证研究思路的研究生，以及想把旧统计流程替换成 Python 的业务人员。

2. 把数据和环境备齐：NetCDF 读取、站点表整理与缺测处理

实习数据最常见就两种形态：格点再分析资料的 NetCDF 文件，或者一张按站点排列的 CSV 站点表。前者空间连续、便于做 EOF 和地图展示，后者贴近观测、适合做单站时间序列和站点差值图。无论哪种，第一步都是把数据读成规整的 DataFrame 或二维数组，并且把缺测问题暴露出来，而不是等算到一半才报错。

2.1 环境最小清单：装到什么程度就够了

Python 3.8 以上版本都能跑，我一般建议直接用 3.10 或 3.11，兼容性最省心。核心库就是 xarray、pandas、numpy、scipy、matplotlib，再加一个 netCDF4 作为 xarray 读取 NetCDF 的后端，画图用 cartopy。做显著性检验时会用到 pymannkendall，这是个轻量库，专门算 M-K 趋势检验。把下面这份 requirements.txt 一次性装完就行。

# 实习环境依赖，一行装完 xarray pandas numpy scipy matplotlib netCDF4 cartopy pymannkendall

装的时候注意两点。第一，cartopy 对底层 geos 库有依赖，Windows 上偶尔装不上，这时候直接去下载对应 Python 版本的 wheel 包离线安装，比折腾 conda 省时间。第二，pymannkendall 不是必须的，如果装不上，手动实现 M-K 检验也就二十行代码，不影响实习进度。我自己的习惯是先装上面的基础七件套，跑通数据读取之后再补最后一个。

2.2 格点数据读取：用 xarray 而不是 netCDF4 裸读

很多初学者拿 netCDF4 库一组变量名一个循环地读，代码又长又容易错。xarray 的好处是变量和坐标直接映射，数据维度的信息一目了然，转 DataFrame 做缺测统计也方便。读文件第一件事不是取数，而是打印数据集结构，看清楚变量名、维度名和坐标范围。

import xarray as xr # 打开格点降水/气温文件，很多实习资料的变量名叫pre或tmp ds = xr.open_dataset("pre.nc") print(ds) # 先看维度名与坐标名，lat可能叫latitude，lon可能叫longitude lat = ds["lat"].values lon = ds["lon"].values pre = ds["pre"].isel(time=0).squeeze() # 取第一个时次并去掉多余维度 # 检查经度范围，如果是0~359.75后面画图前必须处理 print("经度范围:", lon.min(), lon.max()) # 转成DataFrame，方便统计缺测 pre_df = pre.to_dataframe(name="precip") print("缺测格点数:", pre_df["precip"].isna().sum())

这里的 isel(time=0) 是按位置索引取第一个时次，squeeze() 会把长度为一的维度自动压缩掉。to_dataframe 会把坐标变成行索引，便于直接查缺测比例。为什么要先打印数据集？因为不同来源的数据坐标名差异很大，有的叫 latitude，有的叫 lat，有的经度是 0 到 360 度，有的已经是 -180 到 180 度，不先摸清边界，后面画图必然翻车。经度 0 到 360 的问题在第 5 章避坑清单里有专门处理。

2.3 站点数据清洗：日期对齐与缺测过滤

站点数据的坑主要在日期列。很多原始 CSV 里日期是字符串，直接用 resample 会报错，必须先转成 pandas 的时间类型。站点表最常见结构是一列日期、一列站点编号、一列观测值，按日期和站点编号透视成长表，每一行是一个日期，每一列是一个站点，这样后续算气候态、距平、相关都方便。

import pandas as pd df = pd.read_csv("station.csv", parse_dates=["date"]) # 日期列直接解析成时间 df = df[["date", "id", "precip"]] # 只留业务上要用的三列 # 透视成“日期 x 站点”的宽表 df_pivot = df.pivot_table(index="date", columns="id", values="precip") # 删掉缺测比例超过5%的站点，thresh表示最少要有多少有效观测 df_clean = df_pivot.dropna(axis=1, thresh=int(len(df_pivot) * 0.95)) # 日尺度重采样到月平均，M表示月末作为时间戳 df_month = df_clean.resample("M").mean() print(df_month.head())

pivot_table 里 index 和 columns 不能写反，写反了会得到一堆 NaN。dropna 的 thresh 参数是“这一列至少保留多少个非空值”，我用整张表长度的 95% 做门槛，也就是允许站点缺测 5%，超过就删掉整列。resample("M") 里的 M 是月末标签，如果希望时间戳落在月初，可以改成 MS，两者对月均值计算没有影响，只是索引显示不同。重采样之前一定要确认日期是 DatetimeIndex，否则会直接报 “Only valid with DatetimeIndex” 之类的错误，这个坑很常见。

3. 气候态、距平与线性趋势：三个高频统计量的算与判

气候态、距平、线性趋势是气象统计分析实习里出现频率最高的三个量，也是后续合成分析、EOF、相关分析的基础。这三个量本身计算不复杂，但基准期怎么选、趋势单位怎么换算、滑动平均窗口设多大，都是实习报告评审老师喜欢追问的点。这一章把计算代码和参数取舍放在一起讲。

3.1 气候态与距平：基准期直接决定距平正负

气候态就是某一段时期逐月或逐年的平均值，最常见的是 1981 到 2010 年这个 WMO 标准气候期。距平是原始值减去气候态，正距平代表偏多/偏暖，负距平代表偏少/偏冷。基准期选得不一样，同一年的距平符号可能直接反转，所以报告里必须写明基准期。计算时先按月份分组求平均，再逐月做差。

# 取1981-2010年作为气候基准期 base = df_month.loc["1981":"2010"] clim = base.groupby(base.index.month).mean() # 得到1-12月逐月气候态 # 逐月循环计算距平，避免groupby跨年混淆 anom = pd.DataFrame(index=df_month.index, columns=df_month.columns) for m in range(1, 13): anom.loc[df_month.index.month == m] = df_month.loc[df_month.index.month == m] - clim.loc[m] print("1月气候态:", clim.loc[1]) print("距平均值(应接近0):", round(anom.mean().mean(), 6))

这里我刻意没有用 groupby.transform 一行写法，而是用循环按月份逐次做差。原因是 transform 在处理多层索引或存在缺失月份时容易搞混对齐关系，循环虽然多写两行，但每一行做了什么清清楚楚，适合实习报告这种需要反复核对结果的场景。代码最后检查距平均值，正常情况下整段序列距平的平均值应该非常接近 0，如果偏差明显，说明基准期没有正确对齐，或者是数据本身存在跳跃性不连续。

3.2 线性趋势与滑动平均：参数不当容易过度平滑

线性趋势用 scipy.stats.linregress 最省事，它一次性返回斜率、截距、相关系数、p 值和标准误。实习报告里最容易漏掉的是斜率单位：数据是逐月序列，斜率就是“每个月的增量”，要换算成“每十年的增量”需要乘 120，而不是乘 10。

import numpy as np from scipy.stats import linregress # 以序列长度作为自变量 x = np.arange(len(df_month)) # 以某站点或区域平均序列作为因变量 series = df_month["站A"].dropna() x_valid = x[series.index] # 剔除缺测后的时间索引 slope, intercept, r_value, p_value, std_err = linregress(x_valid, series) print("趋势(每月):", slope) print("趋势(每十年):", slope * 120) print("p值:", p_value) p = p_value < 0.05 # False不显著

linregress 返回的 p 值默认是双尾检验，p 小于 0.05 通常认为趋势显著。业务报告习惯同时给“倾向率”和“显著性水平”，只说趋势大小不给 p 值等于没做完。另一个坑是缺测：linregress 遇到 NaN 会直接返回 NaN，所以先 dropna 再回归，同时把 x 也筛掉对应位置，否则 x 和 y 长度对不上。

滑动平均的经验法则是：逐月数据用 5 或 7 个月窗口，逐日数据用 9 到 31 天窗口，窗口越大曲线越光滑但真实信号衰减越厉害。实现上关键是 center=True，这样平滑后的曲线不会向右偏移。

# 5个月滑动平均，center=True保持曲线与原始序列对齐 smooth = df_month["站A"].rolling(window=5, center=True, min_periods=3).mean() # 验证：滑动平均序列与原始序列的相关系数不应低于0.9 print("相关系数:", df_month["站A"].corr(smooth))

window 是窗口长度，center=True 表示窗口中心对齐当前位置，min_periods=3 表示窗口内至少有 3 个有效值才计算，这样首尾不会出现一大片 NaN。有同学图省事不写 min_periods，结果序列前后各缺两个月，画图时首尾断线，报告观感很差。窗口选太大曲线倒是平滑了，但极端事件的信号也被抹掉，这一点第 5 章避坑里会再强调。

4. 显著性检验与EOF分解：实习报告里最硬的骨头

合成分析、相关分析、EOF 分解，这三块基本决定了实习报告的上限。很多同学公式背得熟，一到代码就发怵，其实代码反而是最好复现的部分。真正需要花心思的是两个问题：样本独立性怎么处理，以及 EOF 模态怎么解释才不牵强。

4.1 t 检验与合成分析：检验之前先想清楚自由度

合成分析是气象里最常用的显著性检验场景：选出一组厄尔尼诺年、一组正常年，比较两组降水距平差异是否显著。实现用 scipy.stats.ttest_ind 即可。要注意的坑不是 Welch 还是 Student，而是样本独立性，详见第 5 章。

from scipy.stats import ttest_ind # 强厄尔尼诺年与正常年的降水距平，列名是年份 el_nino = anom.loc[["1997", "2015", "1982"], "站A"] normal = anom.loc[["2000", "2001", "2005"], "站A"] t_stat, p_value = ttest_ind(el_nino, normal, equal_var=False, nan_policy="omit") print("t统计量:", t_stat, "p值:", p_value) # 结果为显著的标志 sig_flag = p_value < 0.05

equal_var=False 对应 Welch 检验，两组方差不等时更稳健。nan_policy="omit" 会在计算时自动剔除 NaN，省去手动清理的麻烦。如果做的是逐格点 t 检验，建议把 t 统计量也保存下来，报告里可以附上 t 值分布图，比只贴一张 p 值图更有说服力。还需要留意的是，年份列表必须和 anom 的索引类型一致，如果索引是 Timestamp，这里要改成 anom.loc[["1997-01-01","2015-01-01"]]，否则直接报 KeyError。

M-K 趋势检验在实习里也常用，它不要求序列服从正态分布，比线性回归检验更稳健。pymannkendall 库的接口很友好：

from pymannkendall import original_test result = original_test(series.dropna()) print("趋势方向:", result.trend) # increasing / decreasing print("p值:", result.p) print("Sen斜率:", result.slope)

result.trend 直接给出 “increasing” 或 “decreasing” 字符串，Sen 斜率是趋势的稳健估计，适合报告里直接引用。需要强调的是，M-K 检验的结论要写“通过 M-K 检验，序列存在显著上升趋势”，而不是“序列存在线性上升趋势”，两种检验的假设不同。

4.2 用 SVD 手写 EOF 分解：比 sklearn 更透明

EOF 分解在数学上就是矩阵的奇异值分解：把距平场 X 拆成时间系数 PC 和空间型 EOF。很多教程推荐用 sklearn 的 PCA，但 PCA 默认做了中心化和归一化，返回的对象还要再拿属性，不如直接调 numpy.linalg.svd 来得透明。这里给出一个实习可用的完整实现。

import numpy as np # anom 是 (时间, 格点) 的距平场，行=时间，列=空间点 X = anom.values.astype(np.float64) # 实习阶段对缺测做简单补零；严肃研究建议用插值或只选完整格点 X = np.nan_to_num(X, nan=0.0) # 中心化：每列减去该格点时间均值，EOF要求输入是距平或标准化场 X = X - X.mean(axis=0) # 做奇异值分解：X = U @ diag(S) @ VT U, S, VT = np.linalg.svd(X, full_matrices=False) # 前3个模态的空间型，VT的行向量即EOF模态 eof = VT[:3, :] # 时间系数：U的前3列乘以对应奇异值 pc = U[:, :3] * S[:3] # 方差贡献率：S的平方占比 variance_frac = (S**2) / np.sum(S**2) print("前3个模态方差贡献:", variance_frac[:3])

X 的形状必须确认是 (时间, 格点)，如果反了算出来的 EOF 就成了时间序列，画图时所有格点全是零，这是最高频的操作失误。np.linalg.svd 的 full_matrices=False 表示做经济型分解，只返回有用部分，实习数据规模下计算毫无压力。EOF 空间型是无量纲的，不要直接写“EOF1 的单位是 mm”，正确说法是“EOF1 的空间分布型”。想要看有物理单位的空间型，把该模态的 PC 与原始距平场做重建，也就是 np.outer(pc[:, 0], eof[0, :])，重建结果与原始距平同单位，这个技巧在 6.2 小节会用到。

模态数怎么取？实习报告一般取前 3 个就够。判断模态是否“可分”用 North 准则：相邻模态的方差贡献误差区间发生交叠时，这两个模态在统计上不可区分。误差近似为 λ * sqrt(2/N)，N 是时间样本数。如果 EOF1 贡献 22%，N=60，误差约为 22% * sqrt(2/60) ≈ 4%，而 EOF2 恰好是 18%，两者区间交叠，那说明 EOF1 和 EOF2 可能是一对耦合模态，不能硬拆。这段逻辑写进报告是明显的加分项。

4.3 从 EOF 结果到物理解释：不要硬套遥相关

EOF1 在降水场上通常表现为全区一致型，对应大尺度背景场异常；EOF2 常见南北反相型，对应副热带高压的位置变化。解释时尊重符号约定：空间型数值乘时间系数，PC 为正时，空间型正值区表示偏多/偏暖，负值区表示偏少/偏冷。经常有人把空间型和 PC 单独拿出来各讲各的，空间型画了一堆，PC 序列也画了一堆，就是不说明“哪个阶段正异常最明显”，报告读起来像两张不相关的图。

这里给一个既严谨又不费力的解释套路：先讲模态方差贡献，再讲空间型正负大值中心的位置，再挑 PC 序列里绝对值最大的 2 到 3 个年份做典型个例，最后用一句话概括。例如“EOF1 解释方差 34%，空间型呈全区一致分布，PC 在 2015 年达到极大正位相，对应当年降水整体偏多”。这样一套逻辑下来，评审老师想追问的点基本都被你提前堵住了。

5. 避坑清单：五个气象统计实习的典型翻车现场

这一章列五个我在实际处理气象数据时反复踩过的坑，全部来自真实数据，每一条都按“现象、原因、解决”来讲。建议把这份清单打印出来贴显示器旁边，写报告前逐条核对一遍。

5.1 读取与数据格式的坑

坑一：从 NetCDF 读出的经纬度全是 0 到 359.75，画图时亚洲出现在地图最右侧。

现象是画出来的等值线图整体错位，中国区域被劈成两半。原因是很多格点数据集的经度用 0 到 360 表示，而 Cartopy 默认用 -180 到 180。解决方法是先把经度坐标重排：把经度大于 180 的格点减去 360，然后重新按新经度排序。

lon = np.where(lon > 180, lon - 360, lon) # 按新的经度顺序重排数组 order = np.argsort(lon) lon = lon[order] field = field[:, order]

注意执行顺序：先换经度，再排序，最后才能用 np.meshgrid 去画图。很多同学先画图再想办法调坐标，就晚了。如果数据本身是全球格点，重排后还要确认起始经度是否正好落在 -180 上，相差一点点不影响实习报告，但不影响不代表可以忽略，画出来边界整齐是第一步。

坑二：CSV 里的日期列是字符串，resample 直接报 TypeError。

现象是 resample 或 groupby 时报 “Only valid with DatetimeIndex”。原因是 pandas 不知道字符串是日期。解决就是读入时显式指定 parse_dates：pd.read_csv("data.csv", parse_dates=["date"])。如果读进来已经是 DataFrame，也可以用df["date"] = pd.to_datetime(df["date"])。这个坑本身不值钱，但特别消耗时间，尤其是日期格式写成 “1997/01/01” 或者 “1997-1-1” 时，pd.to_datetime 能自动解析，但 parse_dates 偶尔会漏，建议两种方式都写一句。

5.2 统计计算里的坑

坑三：EOF 分解前没处理缺测，SVD 结果全是 NaN。

现象是打印 PC 和 EOF 矩阵，里面全是 NaN，没有任何有用数字。原因是 SVD 不允许输入矩阵含缺失值，NaN 进入线性代数运算后会直接传染整个矩阵。解决方法是进入 SVD 前先检查：np.isnan(X).sum()，如果缺测格点不多，用该格点时间均值填充；如果缺测比例较高，只选择完整格点参与分解。实习阶段用 nan_to_num 补零可以接受，但报告里要写明补零处理，否则审阅者会质疑结论。我曾经为省事直接补零，结果 EOF1 的方差贡献被明显压低，因为零方差被拉进了第一主轴。

坑四：t 检验结果一片显著，但业务上根本没法解释。

现象是合成分析里 p 值小于 0.05 的格点密密麻麻，但降水异常的空间分布东一块西一块，没有物理意义。原因是气象时间序列有很强的自相关，样本量 n 并不等于独立样本数。有效样本数大约为 n_eff ≈ n * (1 - r) / (1 + r)，r 是时间序列的一阶自相关系数。n=30、r=0.6 时，有效样本数只有 7.5，自由度被严重高估，t 检验自然容易“假显著”。解决方法是先用自相关函数估算 n_eff，再用 n_eff 结果做检验；或者改用蒙特卡洛置换检验，把年份标签随机打乱 1000 次，看观测统计量落在零分布的什么位置。实习报告里能主动提到有效样本数，已经是超出平均水平的严谨程度了。

5.3 画图与出图的坑

坑五：Cartopy 画海岸线时长时间卡住或直接报错。

现象是 add_feature 或 coastlines() 执行后界面一直空白，最后报一个下载或文件缓存错误。原因是 Cartopy 首次使用会联网获取自然地球数据，网络环境不好时就会卡住。解决方法是提前把缓存数据准备好，或者改用离线 shapefile。只在实习报告里画简单地图的话，用ax.coastlines(linewidth=0.5)这一个命令就够了，它依赖的数据量小，出图快，比加载全套 LAND 和 OCEAN 要素稳定得多。需要填陆地底色时再考虑 add_feature，但要做好首次运行较慢的心理准备。

6. 用地图把结果讲清楚：Cartopy 填图、打点与批量出图

统计结果最终要靠地图呈现，实习报告里最常见的组合是填色图加显著性打点。填色图的配色选择有讲究：距平场用红蓝双色（RdBu_r），0 在中间，正值红色、负值蓝色；原始量场用连续渐变暖色更合适。这一章给出一套直接可用的模板。

6.1 一张可复用的填色图模板

下面代码里，field 是待绘制的二维场，lon 和 lat 是一维坐标，LON、LAT 是网格化后的二维坐标。显著性打点用网格化数组做布尔掩码，坐标对应关系不会错位。

import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 网格化经度纬度，shape与field一致 LON, LAT = np.meshgrid(lon, lat) fig, ax = plt.subplots(figsize=(8, 5), subplot_kw={"projection": ccrs.PlateCarree()}) ax.coastlines(linewidth=0.5) ax.add_feature(cfeature.LAND, color="#eeeeee", zorder=0) # levels=21让0值位于色条中间，RdBu_r适合距平场 im = ax.contourf(LON, LAT, field, levels=21, cmap="RdBu_r", transform=ccrs.PlateCarree()) plt.colorbar(im, ax=ax, shrink=0.8, label="precip anomaly (mm)") # 显著性打点，只画p<0.05的格点 sig_mask = p_value < 0.05 ax.scatter(LON[sig_mask], LAT[sig_mask], s=1, c="black", marker=".", transform=ccrs.PlateCarree()) # 横纵坐标刻度按20度间隔设置，避免全挤在一起 ax.set_xticks(np.arange(60, 141, 20)) ax.set_yticks(np.arange(0, 61, 10)) ax.set_title("EOF1 Precip Anomaly") plt.savefig("eof1.png", dpi=300, bbox_inches="tight")

坐标刻度密集的问题是新手画地图最常遇到的：经纬度从 0 到 359，matplotlib 默认把每个整数值都放上去，整条横轴全是数字。用 set_xticks 手动指定间隔，一次解决。打点用 scatter 而不用 plot，是因为 scatter 对散点数据更高效，s=1 控制点的大小，marker="." 避免点太大糊成一片。savefig 务必加 dpi=300 和 bbox_inches="tight"，否则插入 Word 时清晰度不够，或者图周边留白过多。

6.2 结果自查与批量出图：算完先验证再上报告

EOF 分解完之后一定要做一步自查：用前 K 个模态重建原场，算重建场与原场的空间相关系数。这个数字能直接告诉你模态截断是否合理。

# 用EOF1重建距平场并与原始距平场比较 recon = np.outer(pc[:, 0], eof[0, :]) corr = np.corrcoef(recon.flatten(), anom.values.flatten())[0, 1] print("EOF1重建-原场空间相关:", corr)

如果 EOF1 的解释方差只有 15%，那么这个相关系数通常在 0.3 到 0.5 之间，说明一个模态不足以代表整体场。报告里可以补充累加前 3 个模态后的重建，相关系数一般能到 0.7 以上。这一步的价值在于，它把“EOF 模态方差贡献”和“重建场对原始场的还原能力”两个概念拧在了一起，写出来的结论更站得住脚。

批量出图是实习报告收尾时最节省时间的一招。用循环把前 3 个模态一次性画进一张多子图里，文件名按序号输出，后面写报告时想引用哪张就引用哪张。

fig, axes = plt.subplots(1, 3, figsize=(13, 4), subplot_kw={"projection": ccrs.PlateCarree()}) for i, ax in enumerate(axes): ax.contourf(LON, LAT, eof[i].reshape(LON.shape), levels=15, cmap="RdBu_r", transform=ccrs.PlateCarree()) ax.coastlines(linewidth=0.5) ax.set_title(f"EOF{i+1}") plt.savefig("eof_modes_123.png", dpi=150)

eof 数组是一维展开的格点序列，必须 reshape 成 LON.shape 才能填图，这个细节漏掉会直接报维度不匹配错误。从那以后我每拿到一套气象实习数据，都会先花十分钟把经度范围、缺测比例、时间跨度摸一遍，再开始写计算代码，这个习惯帮我省掉无数次返工。希望这份笔记能帮你少走同样的弯路，祝实习顺利。

本文还有配套的精品资源，点击获取