简介：面向气象专业学生、科研人员及Python数据分析学习者，这是一套紧扣气象统计方法课程实习的代码与数据包，帮助解决从数据读取、清洗到分析可视化的实操入门问题，尤其适合缺少编程经验的初学者快速建立完整分析流程。压缩包共12个文件，整体仅80KB，包含8个按实习顺序编号的Python脚本、2个DAT格式气象数据文件、1个CTL数据描述文件及1个TXT文本数据文件；Python脚本为主要学习对象，覆盖数据加载、Pandas与Numpy统计运算、相关与回归分析、时间序列处理及Matplotlib/Seaborn绘图等关键环节，配套数据文件可直接运行验证。目前已有1570人学习下载，资源体积小、结构清晰，适合作为课程作业参考、实习报告素材或自学气象数据统计的起步模板。脚本按实习任务递进组织，便于逐段对照学习；用户可结合实际站点数据替换练习，快速迁移到自己的分析任务中。
1. 气象统计方法实习：为什么用 Python 重做一遍传统方法
气象统计方法实习（Python版）这门课或这份材料，核心不是教你 Python 语法，而是把气候态、距平、趋势、相关回归、显著性检验这一整套经典流程，从 Fortran 和 Excel 里搬到 Python 生态里跑通。拿到这套实习材料的人，通常是手里有一份站点或格点气象数据，要在几天内交出一份能说服老师的分析报告。这篇笔记按我最常用的路径来写：先配环境、读数据，再做基础统计量，再处理相关与回归，最后单独讲最容易翻车的几个排查点。适合气象、地学相关专业的学生，也适合刚接触 netCDF 和逐日资料的科研新手。你需要的不是更多理论，而是一条少踩坑的落地路径。
2. 环境准备与数据读入：把实习数据变成干净的 Xarray 对象
2.1 Python 环境与气象生态选型：装哪些库、装多新
气象统计实习最忌惮的就是把时间浪费在环境配置上。我的建议是直接用 conda 建一个独立环境，不要动系统基础环境，避免后面装库时污染其他项目。Python 版本选 3.10 而不是最新的 3.12/3.13，原因是 xarray、netcdf4、eofs 这些库在 3.10 上轮子最全，没必要追新版本给自己找不痛快。
conda create -n met python=3.10 -y conda activate met conda install -c conda-forge xarray dask netcdf4 numpy pandas scipy matplotlib -y pip install eofs pymannkendall这套组合里，xarray 负责处理带维度坐标的格点场，dask 让超大数组可以延迟计算，netcdf4 是读取 nc 文件的后端；eofs 用来做 EOF 分解，pymannkendall 用来做 Mann-Kendall 趋势检验。如果不想装 dask，xarray 也能跑，但处理几十 GB 的再分析资料时会明显吃力。
这里有个容易踩的坑：conda 和 pip 混用容易导致依赖版本冲突。我的习惯是核心库用 conda 装，纯 Python 包用 pip 装，装完之后尽量少在两个源之间交叉升级。Windows 用户尤其建议用 conda 创建环境，因为 conda 会自动管理 PATH，省去手动配置 Python 环境变量的过程。装完记得先跑一句which python确认已经切到 met 环境，避免后面实际跑的是系统自带的老版本 Python。
2.2 读入 netCDF 实习数据：三行代码把数据读成 DataArray
实习资料最常见的格式是 netCDF，也就是 .nc 文件。xarray 读这种文件非常顺手，但前提是先把变量名、维度名和单位搞清楚。拿到一个陌生文件，第一件事不是写计算，而是打印看结构。
import xarray as xr # 读取 netCDF 数据，engine 指定后端 ds = xr.open_dataset("era5_temp_1991_2020.nc", engine="netcdf4") print(ds) # 先看变量名、维度名和全局属性 # 按变量名取数，得到一个 DataArray temp = ds["t2m"] print(temp.coords) # 检查 time / latitude / longitude 的具体名称engine="netcdf4"是告诉 xarray 用 netCDF4 库去解析文件，如果你已经通过 conda 装了 netcdf4，这个是稳定选项。如果不指定 engine，xarray 会自动猜测，但遇到扩展名不规范的文件时会报错。print(ds)这一步非常关键，它能直接告诉你维度名是latitude还是lat，经度范围是 0 到 360 还是 -180 到 180。这两个信息决定了后面所有切片和运算怎么写。
如果实习数据是 CSV 或纯文本格式，就用 pandas 读入，常见做法是pd.read_csv("station_data.csv", parse_dates=["time"])。固定宽度的 .dat 文件则用pd.read_fwf。无论是 nc 还是 csv，读完之后都要统一成 xarray 的 DataArray 格式，因为后续的 groupby、resample、corr 操作在 xarray 里最顺手。
2.3 预处理三件事：时间索引、缺失值与单位换算
数据读进来不等于可以直接用。我一般会强制自己做三个检查：时间维度对不对、缺失值有多少、单位是不是国际单位。这三件事不做，后面所有统计结果都可能错得悄无声息。
import xarray as xr temp = xr.open_dataset("era5_temp_1991_2020.nc", engine="netcdf4")["t2m"] # 1. 统一时间索引：如果数据是逐日，先重采样成逐月，减少计算量 temp_monthly = temp.resample(time="1ME").mean("time") print(temp_monthly["time"].values[:5]) # 2. 缺失值概览：打印每个格点的有效样本数 print(temp_monthly.isnull().sum("time").compute()) # 3. 单位换算：开尔文转摄氏度 print("原始单位:", temp.attrs.get("units")) temp_monthly_c = temp_monthly - 273.15resample(time="1ME")表示按月重新采样，1ME是每月最后一天作为时间标签，这是较新 pandas 版本里的写法；如果你用的是旧版本，可能要用"1M"。后面必须跟.mean("time")，否则 resample 只做了分组还没聚合。temp_monthly.isnull().sum("time")会生成一个 lat-lon 场，显示每个格点有多少个月是缺失的，如果某些格点缺失超过一半，后面做相关分析时就要特别小心。
单位换算看起来很傻，但却是血泪经验：很多再分析资料的温度变量是开尔文，实习报告里却要求画摄氏度图。如果不看attrs["units"]直接画图，图上数值全是 280 到 310，老师一眼就能看出来你连单位都没处理。这里还要注意_FillValue的问题——netCDF 文件里的缺测值通常由mask_and_scale参数自动处理，默认会转成 NaN，所以一般不需要手动处理。但如果你用的是 pandas 读文本文件，缺测值常被写成 -9999，需要手动替换成 NaN。
3. 基础统计量：气候态、距平与趋势，从公式到 DataFrame
3.1 气候态怎么算：用 groupby 替代三十年循环
气候态的定义是把某一段时间内同一月份的气温取平均，最常用的是 WMO 标准的 30 年基准期，比如 1991 到 2020。很多初学者的第一反应是写两层 for 循环，遍历年份和月份，然后一层层累加。这个写法不是不行，而是又慢又容易错，代码动不动就三四十行。xarray 里一个groupby就能解决。
# 截取 1991-2020 作为基准期 base = temp_monthly_c.sel(time=slice("1991-01-01", "2020-12-31")) # 逐月气候态：按 time.month 分组，对每组求平均 clim = base.groupby("time.month").mean("time") print(clim)groupby("time.month")会自动为每个月份生成一组，算出 12 个平均值，坐标名自动变成month。这比手写循环快得多，而且当数据量变大时，dask 还能自动并行。需要注意，如果原始数据是逐日资料，一定要先重采样成逐月再做气候态，否则一个月的 30 个值都会被当成独立样本，气候态的方差会被严重低估。
还有一类常见错误：有人用完整时间段求平均当作气候态，而不是用基准期。如果数据本身跨越 1979 到 2023 年，直接对整个时段求月平均得到的“气候态”混合了不同时期的气候背景，做距平分析时会有系统性偏差。实习报告里写的“气候态”应该明确指向你选择的基准期，通常截取 1991 到 2020，这个细节是判断你是否理解概念的关键。
3.2 距平与标准化距平：两种归一化不要混着用
距平是最常用的气象统计量，公式就是一个减法：每个时刻的值减去对应月份的气候态。标准化距平则再除以逐月标准差，得到无量纲的倍数。两者的使用场景完全不同。
# 距平：每个时刻减去对应月份气候态 # 注意月份坐标对齐 anom = temp_monthly_c - clim.sel(month=temp_monthly_c["time.month"]) # 标准化距平：距平除以逐月标准差 monthly_std = base.groupby("time.month").std("time") std_anom = anom / monthly_std.sel(month=temp_monthly_c["time.month"]) print(std_anom)clim.sel(month=temp_monthly_c["time.month"])这一步是把 12 个月的气候态值扩展到与原始时间序列等长，然后才能做减法。有人会直接写anom = temp_monthly_c - clim，这在维度不匹配时会触发广播，得到完全错误的结果。标准差的基准期选择也有讲究：我一般用base也就是 1991 到 2020 的逐月标准差，这样与气候态保持一致。如果实习报告要求做合成分析，标准化距平能消除不同地区变率差异带来的影响；如果只是看某个区域温度随时间的变化，用绝对距平更直观。两种归一化的结果不要混着用在一张图上，坐标单位都不一样。
3.3 线性趋势与 MK 检验：scipy 一行和 pymannkendall 一循环
趋势分析是实习报告里的重头戏。线性趋势用最小二乘拟合一条直线，斜率代表每月的增温速率。Mann-Kendall 检验则是非参数方法，对非正态分布和缺失值更稳健，两者配合使用是比较完整的做法。
import numpy as np from scipy import stats import pymannkendall as mk # 对某个站点的月平均温度序列做趋势分析 ts = anom.sel(latitude=30, longitude=120).values # 一维时间序列 x = np.arange(len(ts)) # 线性趋势：返回斜率、截距、r值、p值、标准误 slope, intercept, r_value, p_value, std_err = stats.linregress(x, ts) print(f"线性趋势: {slope:.4f} ℃/月, p={p_value:.4f}") # M-K 趋势检验 result = mk.original_test(ts) print(f"MK趋势: {result.trend}, p={result.p}")stats.linregress返回的五个值里，slope 的单位是“每单位 x 的变化量”。如果 x 是月份序号，那么 slope 就是每月变化多少摄氏度；如果要换算成十年趋势，需要乘 120。p 值小于 0.05 通常认为趋势显著。mk.original_test返回一个 namedtuple，trend字段给出increasing、decreasing或no trend，p是显著性水平。
对格点数据做趋势时，常见做法是把 lat 和 lon 两个维度压平，写一个双重循环对每个格点调用上面这段代码。循环前记得先用.transpose("lat", "lon", "time")把时间放到最后一维，这样内存访问更连续，速度会快很多。另一个容易忽略的问题：如果原始数据是逐月资料，M-K 检验会受季节循环干扰，结果往往是“不显著”。此时应该先用 3.2 的标准化距平去除季节循环，再跑 MK 检验，或者直接改用mk.seasonal_test。
4. 相关、回归与显著性验证：把统计关系“钉死”
4.1 空间相关场：corrwith 的两个前提条件
相关分析是气象统计实习里最常被“做出假结果”的环节。最常见的一个需求是：用一个气候指数序列（比如某个区域平均海温距平）与整个空间场的每个格点做相关，画出一张相关分布图。xarray 提供了一个非常方便的方法corr，但用之前必须满足两个前提：时间轴严格对齐，缺测值处理干净。
import xarray as xr import numpy as np # 假设 idx 是一个与时间等长的一维序列，比如区域平均海温距平 # temp_monthly_c 是 (time, lat, lon) 的月平均气温距平场 # 先统一时间：两个对象只取公共时间段 common_time = temp_monthly_c["time"].values idx_da = xr.DataArray(idx, dims="time", coords={"time": common_time}) # 方法一：xarray 自带的 corr，沿 time 维计算 corr_map = xr.corr(temp_monthly_c, idx_da, dim="time") print(corr_map.min(), corr_map.max())xr.corr沿dim="time"计算每个格点与 idx 序列的皮尔逊相关系数，结果是 (lat, lon) 场。它要求两个对象的 time 坐标完全对齐，差一个时刻都会让结果变成 NaN。如果数据是逐月但 idx 是逐年的，必须先把 idx 重复扩展到逐月，或者把温度场也降采样到逐年，保证样本一一对应。
第二个前提是缺测值。xr.corr在计算时遇到任意一个缺失值，整个格点的相关系数都会返回 NaN。如果数据里有大量格点缺测，图会一片空白。我的习惯是先对缺测格点做插值，或者只对有效样本数超过阈值（比如总时长的一半）的格点保留结果。这一步做完后，再叠加上显著性检验的结果，就离一张能放进报告里的图不远了。
4.2 一元回归与 t 检验：只给 r 值等于没做完
相关系数只能告诉你两个变量有没有线性关系，不能告诉你这个关系是否可信。气象统计实习里，只画一张相关场而不做显著性检验，是老师最容易扣分的地方。一元回归和相关系数本质上是一体两面，但回归还能给出趋势的斜率，物理意义更直接。
import numpy as np from scipy import stats # 对单个格点做一元线性回归，并计算 t 统计量 def regress_1d(y, x): slope, intercept, r_value, p_value, std_err = stats.linregress(x, y) n = len(x) # 由相关系数换算 t 统计量 t_stat = r_value * np.sqrt((n - 2) / (1 - r_value**2)) return slope, p_value, t_stat # 对空间场循环 x = np.arange(len(idx)) slopes = np.full((temp.shape[0], temp.shape[1]), np.nan) ps = np.full((temp.shape[0], temp.shape[1]), np.nan) for i in range(temp.shape[0]): for j in range(temp.shape[1]): y = temp[i, j, :] if np.all(np.isnan(y)) or np.std(y) == 0: continue slope, p, t = regress_1d(y, x) slopes[i, j] = slope ps[i, j] = p print(np.nanmin(ps), np.nanmax(ps))这段代码里有两个要点：一是if np.std(y) == 0这个判断，如果某个格点的值全相同，linregress 会直接报错或返回 NaN，需要提前跳过；二是 t 统计量可以直接由 r 值计算出来，公式是t = r * sqrt((n-2)/(1-r^2))，但这里直接用 linregress 返回的 p 值更省事。
地学里常用的显著性水平是 0.05 和 0.01。需要注意，空间场每个格点都在做一次检验，几千个格点里即使没有真实关系，也会有大约 5% 的格点 p 值小于 0.05。这就是多重检验问题。严谨的做法是在图上标注通过显著性检验的格点数量，或者用 FDR 方法校正。实习报告里如果能主动提到这一点，是非常加分的。
4.3 多元回归与共线性：什么时候放弃最小二乘
当你要同时考虑多个因子对一个气象变量的影响时，比如海温和气压场一起解释降水，就需要多元线性回归。多元回归能用 statsmodels 一行跑完，但结果很容易被“共线性”坑。
import statsmodels.api as sm # X 是 (n, m) 的预测因子矩阵，y 是 (n,) 的响应序列 # 例如 X 两列分别是海温距平和气压距平 X = sm.add_constant(X) # 加截距项，必须做 model = sm.OLS(y, X).fit() print(model.summary()) print("条件数:", model.condition_number)sm.add_constant(X)是在矩阵最前面加一列 1，对应回归方程的截距。如果不加，statsmodels 默认做无截距回归，很多新手在这里栽过跟头。model.summary()会输出一份完整报告，里面有每个系数的 p 值、R²、F 统计量。condition_number是条件数，通常超过 30 就提示存在较严重的共线性。
共线性的典型症状是：两个预测因子本身高度相关（相关系数超过 0.8），这时最小二乘的系数会互相抵消，单个系数的标准误变得很大，p 值全不显著。实习作业里遇到这种情况，最简单的处理是只选其中一个因子做一元回归，或者对两个因子做标准化后求平均合成一个指数。强行保留两个共线因子，报告里解释系数时也会自相矛盾。如果你确实需要同时考虑多个因子，可以用岭回归或主成分回归，但对实习而言，说明“为什么放弃多元回归”比硬上更得体。
5. 五个常见排查点：从报错到错结果
5.1 现象：整个相关场接近 1，空间上几乎全红
我第一次做空间相关场时，画出来的图全部是深红色，相关系数接近 1，当时还以为是找到了完美关系。后来检查才发现是两个序列的时间轴没对齐——一个是逐月数据从 1 到 12 月，另一个是逐年数据从 1991 到 2020，xarray 在广播时把两个对象按位置对齐了，完全不是同一时刻的对应关系。
原因是corr要求两个对象的时间坐标完全一致，如果坐标不一致但长度相同，xarray 不会报错，只是按位置强算。解决方法是先打印两个对象的 time 坐标，用np.intersect1d找出公共时间，再用.sel(time=common)截取。还有一个习惯值得养成：算完相关系数后，先打印corr_map.max()，如果接近 1，立刻怀疑时间对齐问题，而不是高兴。
5.2 现象：气候态场比单年场还“乱”
原本应该平滑的逐月气候态，画出来之后相邻月份差异巨大，甚至比某一年的实际值还要“野”。这个问题几乎都出在基准期的选择上：有人对完整时间段（比如 1979 到 2023）做了平均，而不是对 1991 到 2020 基准期平均。
原因是气候态的物理意义是某个基准期的平均状态。如果数据包含异常的暖年或冷年，整段平均会把它们混进去。解决方法是严格按base = temp.sel(time=slice("1991-01-01", "2020-12-31"))截取后再 groupby。另外，如果你用的是逐日资料且没有先做月平均，同一个月内 30 个日值会被当成独立样本来平均，也会导致气候态被极端日值拉偏。先 resample 到月，再做气候态，顺序不能反。
5.3 现象：MK 检验永远不显著
把一段气温序列丢给mk.original_test，p 值常年大于 0.05，哪怕线性回归已经显示明显增温趋势。这让很多人怀疑 pymannkendall 是不是装错了。其实原因是 MK 检验的原假设是序列独立且同分布，而气温序列有极强的季节性周期，1 月冷 7 月热，这种周期波动被 MK 检验当成了噪声。
解决方法是先做季节循环去除。用 3.2 节的标准化距平代替原始序列再跑 MK，或者直接用mk.seasonal_test，它能对每个季节分别检验。实习报告里如果数据是逐月，我建议优先用seasonal_test，它的输出里同样有 trend 和 p 字段，不需要额外改动。
5.4 现象：温度场数值在 250 到 310 之间，单位却标着摄氏度
拿到一份 nc 文件，打印出来的温度值在 280 左右，画出来的等值线图全是几百，明显不是摄氏度量纲，但变量的 units 属性里写的是 K。更坑的是，有些数据集在 units 属性里写的就是 Kelvin，但数据值已经减过 273.15，完全凭运气。
原因是不同来源的数据集标准不统一。解决方法是打印temp.attrs["units"]，用temp - 273.15做换算；如果单位是 K，同时数值也接近 300，那就确认是开尔文。如果单位写着 C 但数值是 300，那就是数据本身有问题，应该回溯数据源。这个检查放在预处理的第一步，后面所有统计量才不会带病工作。
5.5 现象：Windows 下读不了含中文路径的 nc 文件
在 Windows 上把实习数据放在D:\气象数据\目录下，用 xarray 读取时直接报错，提示文件不存在，但路径明明是对的。原因是 netCDF 库在 Windows 下对非 ASCII 路径的编码支持不好，中文路径会被错误解析。
解决方法是把数据文件移动到纯英文路径下，比如D:\met_data\，文件名也用英文。如果条件不允许移动，可以先用pathlib.Path转成绝对路径再传给 xarray。还有一个小技巧：用os.listdir列出目录时，如果文件名带中文，最好在读取前打印一下实际字符串，确认不是乱码。这个问题不算难，但卡住过不少人，提前把工作目录定成纯英文能省下不少时间。
6. 把实习报告做成能复现的结果：EOF、功率谱与交叉验证
基础统计量做完之后，实习报告如果想拿高分，通常需要加一两个进阶分析。我一般会推荐同学们从三个方向里选一个：EOF 分解、功率谱分析、交叉验证。它们分别回答“空间上最主要的变化型是什么”“时间上有没有显著周期”“统计关系是不是过拟合”。
EOF 用 eofs 库写起来很简洁，但有一个容易翻车的点：EOF 的空间模态符号是任意的，同一个模态乘 -1 仍是合法解。画图时必须检查时间系数和空间模态的符号匹配，否则会出现“空间模态是正的，时间系数却是负的”这种物理上矛盾的图。交叉验证的做法是把数据分成前 20 年和后 10 年，用前 20 年做回归，再用后 10 年检验预测效果。这个步骤在实习报告里写出来，老师会觉得你真的理解统计检验的意义。功率谱分析则可以用scipy.signal.welch对去趋势、去季节循环后的序列做谱估计，看是否存在显著的周期峰。
这三项工作没有想象中复杂，核心代码量都在二三十行以内，但能显著提升报告质量。我现在拿到任何一套实习数据，都会先花五分钟做一张时间序列和空间场的“扫雷图”，看看有没有明显异常值、时间轴是否连续、单位是否合理，然后再决定下一步用哪种统计方法。这个习惯帮我避开了无数次返工，希望你也能用上。
本文还有配套的精品资源，点击获取
