★ FEATURED ARTICLE 简介面向气象数据统计分析学习者的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可能叫latitudelon可能叫longitude lat ds[lat].values lon ds[lon].values pre ds[pre].isel(time0).squeeze() # 取第一个时次并去掉多余维度 # 检查经度范围如果是0~359.75后面画图前必须处理 print(经度范围:, lon.min(), lon.max()) # 转成DataFrame方便统计缺测 pre_df pre.to_dataframe(nameprecip) print(缺测格点数:, pre_df[precip].isna().sum())这里的 isel(time0) 是按位置索引取第一个时次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(indexdate, columnsid, valuesprecip) # 删掉缺测比例超过5%的站点thresh表示最少要有多少有效观测 df_clean df_pivot.dropna(axis1, threshint(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(indexdf_month.index, columnsdf_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 天窗口窗口越大曲线越光滑但真实信号衰减越厉害。实现上关键是 centerTrue这样平滑后的曲线不会向右偏移。# 5个月滑动平均centerTrue保持曲线与原始序列对齐 smooth df_month[站A].rolling(window5, centerTrue, min_periods3).mean() # 验证滑动平均序列与原始序列的相关系数不应低于0.9 print(相关系数:, df_month[站A].corr(smooth))window 是窗口长度centerTrue 表示窗口中心对齐当前位置min_periods3 表示窗口内至少有 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_varFalse, nan_policyomit) print(t统计量:, t_stat, p值:, p_value) # 结果为显著的标志 sig_flag p_value 0.05equal_varFalse 对应 Welch 检验两组方差不等时更稳健。nan_policyomit 会在计算时自动剔除 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, nan0.0) # 中心化每列减去该格点时间均值EOF要求输入是距平或标准化场 X X - X.mean(axis0) # 做奇异值分解X U diag(S) VT U, S, VT np.linalg.svd(X, full_matricesFalse) # 前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_matricesFalse 表示做经济型分解只返回有用部分实习数据规模下计算毫无压力。EOF 空间型是无量纲的不要直接写“EOF1 的单位是 mm”正确说法是“EOF1 的空间分布型”。想要看有物理单位的空间型把该模态的 PC 与原始距平场做重建也就是 np.outer(pc[:, 0], eof[0, :])重建结果与原始距平同单位这个技巧在 6.2 小节会用到。模态数怎么取实习报告一般取前 3 个就够。判断模态是否“可分”用 North 准则相邻模态的方差贡献误差区间发生交叠时这两个模态在统计上不可区分。误差近似为 λ * sqrt(2/N)N 是时间样本数。如果 EOF1 贡献 22%N60误差约为 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_datespd.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 是时间序列的一阶自相关系数。n30、r0.6 时有效样本数只有 7.5自由度被严重高估t 检验自然容易“假显著”。解决方法是先用自相关函数估算 n_eff再用 n_eff 结果做检验或者改用蒙特卡洛置换检验把年份标签随机打乱 1000 次看观测统计量落在零分布的什么位置。实习报告里能主动提到有效样本数已经是超出平均水平的严谨程度了。5.3 画图与出图的坑坑五Cartopy 画海岸线时长时间卡住或直接报错。现象是 add_feature 或 coastlines() 执行后界面一直空白最后报一个下载或文件缓存错误。原因是 Cartopy 首次使用会联网获取自然地球数据网络环境不好时就会卡住。解决方法是提前把缓存数据准备好或者改用离线 shapefile。只在实习报告里画简单地图的话用ax.coastlines(linewidth0.5)这一个命令就够了它依赖的数据量小出图快比加载全套 LAND 和 OCEAN 要素稳定得多。需要填陆地底色时再考虑 add_feature但要做好首次运行较慢的心理准备。6. 用地图把结果讲清楚Cartopy 填图、打点与批量出图统计结果最终要靠地图呈现实习报告里最常见的组合是填色图加显著性打点。填色图的配色选择有讲究距平场用红蓝双色RdBu_r0 在中间正值红色、负值蓝色原始量场用连续渐变暖色更合适。这一章给出一套直接可用的模板。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(linewidth0.5) ax.add_feature(cfeature.LAND, color#eeeeee, zorder0) # levels21让0值位于色条中间RdBu_r适合距平场 im ax.contourf(LON, LAT, field, levels21, cmapRdBu_r, transformccrs.PlateCarree()) plt.colorbar(im, axax, shrink0.8, labelprecip anomaly (mm)) # 显著性打点只画p0.05的格点 sig_mask p_value 0.05 ax.scatter(LON[sig_mask], LAT[sig_mask], s1, cblack, marker., transformccrs.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, dpi300, bbox_inchestight)坐标刻度密集的问题是新手画地图最常遇到的经纬度从 0 到 359matplotlib 默认把每个整数值都放上去整条横轴全是数字。用 set_xticks 手动指定间隔一次解决。打点用 scatter 而不用 plot是因为 scatter 对散点数据更高效s1 控制点的大小marker. 避免点太大糊成一片。savefig 务必加 dpi300 和 bbox_inchestight否则插入 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), levels15, cmapRdBu_r, transformccrs.PlateCarree()) ax.coastlines(linewidth0.5) ax.set_title(fEOF{i1}) plt.savefig(eof_modes_123.png, dpi150)eof 数组是一维展开的格点序列必须 reshape 成 LON.shape 才能填图这个细节漏掉会直接报维度不匹配错误。从那以后我每拿到一套气象实习数据都会先花十分钟把经度范围、缺测比例、时间跨度摸一遍再开始写计算代码这个习惯帮我省掉无数次返工。希望这份笔记能帮你少走同样的弯路祝实习顺利。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
