简介这份资源是面向气象、大气科学及相关专业学生的Python气象统计实习配套材料对应「气象统计实习一」的完整实践内容适合正在学习气象数据处理与可视化、需要动手完成课程作业或实验报告的学习者。压缩包共74个文件以72张png结果图、1个py脚本和1个dat数据文件为主整体约4.38MB其中脚本负责数据读取与统计计算dat文件提供原始气象场数据png则对应气候场、均方差、距平场等分析结果便于对照代码与输出。目前已有1298人学习下载说明该实习案例在同类课程中具有一定参考价值。通过这份材料读者可以了解气象数据从加载、清洗到统计分析与图形呈现的完整流程借助现成脚本和结果图快速核对计算思路理解气候场、距平场、均方差等统计量的含义与展示方式为后续气象预报、气候研究等方向的数据分析实践提供可复用的参考模板。1. 气象统计实习一.rar一份能直接跑通的气象Python实操包如果你手头正好有一份叫「气象统计实习一.rar」的压缩包却不确定里面是什么、能不能直接拿来练手那这篇笔记就是写给你的。它本质上是一套面向气象数据处理入门到进阶的实习材料核心围绕气象Python展开覆盖了从原始观测数据读取、质量控制、统计量计算到可视化出图的一条龙流程。适合气象、大气科学、环境科学方向的学生也适合刚转行做气象数据分析、需要快速补齐实操能力的开发者。我见过太多人拿到这类资源后卡在环境配置和编码问题上最后把压缩包扔进硬盘吃灰——其实只要把几个关键环节打通这份材料能省掉你至少两周的摸索时间。2. 环境准备与数据读取把 .rar 里的东西跑起来2.1 为什么气象Python首选这套组合气象数据有个特点格式杂、维度多、缺测值满天飞。常见的站点观测是文本表格再分析资料是 NetCDF 或 GRIB卫星产品又可能是 HDF。纯靠 pandas 读表格能做一部分但一旦涉及多维数组和投影信息就力不从心。所以气象Python的事实标准组合是 xarray netCDF4 numpy pandas matplotlib再配上 cartopy 做地图投影。这套组合的好处是xarray 把带标签的多维数组封装得跟 pandas 一样顺手netCDF4 负责底层解码cartopy 解决经纬度网格画图时的投影变形问题。「气象统计实习一.rar」里的脚本大概率就是基于这套栈写的。我一般会先建一个干净的虚拟环境避免和系统里已有的包打架。conda 和 venv 都行气象圈用 conda 更多因为 cartopy 的依赖GEOS、PROJ用 conda 装最省心。# 创建独立环境Python 版本建议 3.9 到 3.11太新反而有些气象包没轮子 conda create -n meteo_stat python3.10 -y conda activate meteo_stat # 核心包一次装齐cartopy 走 conda-forge 渠道 conda install -c conda-forge xarray netcdf4 numpy pandas matplotlib cartopy scipy -y # 如果要用 GRIB 格式再加 cfgrib 和 eccodes conda install -c conda-forge cfgrib eccodes -y装完之后别急着跑脚本先验证一下 cartopy 能不能正常画图因为它的投影库最容易出问题。跑一句python -c import cartopy.crs as ccrs; print(ccrs.PlateCarree())不报错就说明底层依赖没问题。2.2 读取站点观测与格点数据的两种写法实习材料里通常有两类数据一类是站点观测的 CSV 或 Excel另一类是格点化的 NetCDF。读取方式完全不同但目标一致——转成带时间维和空间维的数组。站点数据用 pandas 读重点是把日期列解析成 datetime 索引把缺测值常见的是 9999、32766 或空字符串替换成 NaN。这一步不做后面算平均值全是坑。import pandas as pd import numpy as np # 读取站点观测假设文件是逗号分隔含日期、站号、经纬度、气温 df pd.read_csv(station_obs.csv, encodingutf-8) # 日期列转 datetime并设为索引 df[date] pd.to_datetime(df[date], format%Y-%m-%d) df df.set_index(date) # 气象数据里 9999 常代表缺测统一替换成 NaN df[temp] df[temp].replace([9999, 32766, -9999], np.nan) # 按站号分组算月平均resample 的 M 表示按月重采样 monthly df.groupby(station_id)[temp].resample(M).mean() print(monthly.head())格点数据用 xarray 读一行就能打开而且自带坐标信息。关键是理解它的维度顺序通常是 time、level、lat、lon但不同来源可能不一样读进来第一件事是print(ds)看清楚。import xarray as xr # 打开 NetCDFdecode_times 自动解析时间坐标 ds xr.open_dataset(reanalysis_temp.nc, decode_timesTrue) # 先看结构确认维度名和变量名 print(ds) # 假设变量叫 t2m取 2020 年夏季做区域平均 summer ds[t2m].sel(timeslice(2020-06-01, 2020-08-31)) region_mean summer.sel(latslice(20, 40), lonslice(100, 120)).mean(dim[lat, lon]) print(region_mean.values)提示xarray 的sel用标签索引isel用位置索引。经纬度如果是递减的从北到南slice(20, 40)会返回空得先sortby(lat)或者反过来写slice(40, 20)。这个坑我踩过不止一次。3. 统计量计算与质量控制别让缺测值毁掉结果3.1 气象常用统计量的正确算法实习材料里一定会涉及均值、方差、相关系数、趋势系数这些基础统计量。但气象数据有两个特殊性一是缺测值不能简单当 0 处理二是时间序列往往有季节循环算趋势前得先去掉季节项否则趋势会被季节波动淹没。以线性趋势为例很多人直接对原始序列做最小二乘拟合结果斜率被季节变化带偏。正确做法是先算距平减去气候态或者用去季节后的序列。气候态一般用多年同期的平均值比如 30 年气候态。# 假设 ds 里有一个时间序列变量 ts维度是 time # 先算气候态按“年积日”分组求平均 clim ts.groupby(time.dayofyear).mean(dimtime) # 算距平原始值减去对应日期的气候态 anom ts.groupby(time.dayofyear) - clim # 对距平序列做线性趋势用 scipy 的最小二乘 from scipy import stats import numpy as np # 把时间转成数值年为单位方便回归 years anom[time].dt.year (anom[time].dt.dayofyear - 1) / 365.25 slope, intercept, r_value, p_value, std_err stats.linregress(years, anom.values) print(f趋势系数: {slope:.4f} 单位/年, p值: {p_value:.4f})这里groupby(time.dayofyear)是关键它把每年同一天的数据归到一组算出来的气候态才有物理意义。如果数据跨了闰年2 月 29 日会多出来xarray 默认会处理但最好检查一下有没有异常。3.2 质量控制三个必须做的检查气象观测数据没有经过质量控制就直接统计结论基本不可信。实习材料里如果有质控环节通常包括界限值检查、内部一致性检查和空间一致性检查。我一般至少做三件事第一界限值检查。气温超出 -60 到 60 摄氏度、风速为负、相对湿度超过 100% 的一律标记为可疑。第二重复值检查。同一站同一时刻出现多条记录可能是传输重复要去重。第三时间连续性检查。相邻两个时次差值过大比如气温一小时跳变 15 度以上要么是锋面过境要么是仪器故障需要结合天气形势判断。# 界限值检查把超限的值设为 NaN df.loc[(df[temp] -60) | (df[temp] 60), temp] np.nan df.loc[df[rh] 100, rh] np.nan df.loc[df[wind_speed] 0, wind_speed] np.nan # 重复值检查按站号和日期去重保留第一条 df df.reset_index().drop_duplicates(subset[station_id, date], keepfirst).set_index(date) # 时间连续性检查算相邻时次差值超过阈值的标记 df[temp_diff] df.groupby(station_id)[temp].diff().abs() df.loc[df[temp_diff] 15, temp] np.nan注意质控阈值不是固定的不同气候区差别很大。沙漠地区日温差能到 30 度用 15 度阈值会误杀大量正常数据。实际项目里我会按站点历史分位数动态定阈值实习阶段用固定值先跑通流程即可。3.3 缺测值处理插补还是剔除缺测值怎么处理取决于后续分析目的。如果只是算个区域平均少量缺测直接跳过skipnaTrue影响不大。但如果要做时间序列的谱分析或滞后相关缺测会打断序列必须插补。常见插补方法有三种线性插值适合缺测段很短连续缺 1 到 3 个时次的情况气候态插补适合缺测段较长但不需要精确值的情况空间插补用邻近站点的值来补适合站点密集的区域。我一般优先用线性插值超过 6 小时才考虑气候态。# 线性插值limit3 表示最多连续插 3 个缺测 df[temp] df.groupby(station_id)[temp].apply( lambda x: x.interpolate(methodlinear, limit3) ).reset_index(level0, dropTrue) # 剩余的长缺测段用该站该月的气候态填充 monthly_clim df.groupby([station_id, df.index.month])[temp].transform(mean) df[temp] df[temp].fillna(monthly_clim)这段代码里groupby加apply的组合是处理站点数据的常用套路保证插值在每个站点内部独立进行不会跨站串数据。transform(mean)返回的是与原索引对齐的序列直接 fillna 不会错位。4. 可视化与出图让统计结果能直接放进报告4.1 时间序列图与趋势线气象统计结果最终要落到图上。时间序列图是最基础的但有几个细节决定图能不能用横轴日期格式要可读缺测段要断开而不是连成直线趋势线要标注显著性。import matplotlib.pyplot as plt import matplotlib.dates as mdates fig, ax plt.subplots(figsize(10, 4)) # 画原始序列缺测处自动断开 ax.plot(anom[time], anom.values, colorsteelblue, linewidth0.8, label距平) # 画趋势线 trend_line intercept slope * years ax.plot(anom[time], trend_line, colorcrimson, linewidth1.5, labelf趋势 {slope:.3f}/年) # 日期格式主刻度按年次刻度按季度 ax.xaxis.set_major_locator(mdates.YearLocator(2)) ax.xaxis.set_major_formatter(mdates.DateFormatter(%Y)) ax.xaxis.set_minor_locator(mdates.MonthLocator(bymonth[1,4,7,10])) ax.axhline(0, colorgray, linewidth0.5, linestyle--) ax.set_xlabel(年份) ax.set_ylabel(气温距平 (℃)) ax.legend() plt.tight_layout() plt.savefig(trend.png, dpi300)mdates.YearLocator(2)表示每两年一个主刻度数据跨度短的话改成 1。dpi300是投稿或报告的基本要求屏幕预览用 100 就行。4.2 空间分布图cartopy 投影与色标格点数据的空间分布图必须用 cartopy否则经纬度直接当直角坐标画高纬度地区会严重变形。核心是三步建投影、加地图特征、画填色图。import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(8, 6)) # 用 PlateCarree 投影适合中低纬度区域 ax plt.axes(projectionccrs.PlateCarree()) # 画填色图transform 指定数据的原始投影 cf ax.contourf(ds[lon], ds[lat], region_mean.values, levels20, cmapRdBu_r, transformccrs.PlateCarree()) plt.colorbar(cf, axax, orientationvertical, shrink0.8, label气温 (℃)) # 加海岸线和国界 ax.add_feature(cfeature.COASTLINE, linewidth0.5) ax.add_feature(cfeature.BORDERS, linewidth0.3, linestyle:) # 设置范围 ax.set_extent([100, 120, 20, 40], crsccrs.PlateCarree()) plt.savefig(spatial.png, dpi300)levels20控制等值线密度太多会糊成一片太少则看不出梯度。cmapRdBu_r是气温距平图的经典配色红色暖、蓝色冷反向加_r。如果画的是降水改用BrBG或YlGnBu更合适。提示cartopy 第一次画图会下载海岸线数据如果网络不通会卡住。可以提前手动下载 Natural Earth 数据放到缓存目录或者用cfeature.NaturalEarthFeature指定本地路径。5. 避坑与排查那些让脚本跑不起来的常见问题5.1 编码报错中文路径和 CSV 编码现象pd.read_csv报UnicodeDecodeError或者路径里有中文时xr.open_dataset直接找不到文件。原因Windows 下 CSV 常见 GBK 编码pandas 默认 UTF-8中文路径在某些库的底层 C 接口里处理有问题。解决读 CSV 时显式指定encodinggbk或encodingutf-8-sig。路径尽量用英文或者用pathlib.Path包一层再传进去。如果压缩包解压后文件名是乱码用 7-Zip 重新解压并指定编码。5.2 经纬度顺序反了图上下颠倒现象画出来的空间图南北颠倒或者sel选区域返回空数组。原因NetCDF 里纬度可能是从 90 到 -90 递减存储而 cartopy 和大部分绘图库默认纬度递增。解决读进来先检查ds[lat].values[0]和[-1]如果第一个大于最后一个执行ds ds.sortby(lat)。lon 同理但 lon 还涉及 0 到 360 和 -180 到 180 的转换用ds[lon] (ds[lon] 180) % 360 - 180再排序。5.3 时间坐标解析失败decode_times 报错现象xr.open_dataset时decode_timesTrue抛出ValueError提示无法解析时间单位。原因有些 NetCDF 的时间单位写的是hours since 1900-01-01 00:00:00但实际存储的是浮点数且缺少calendar属性xarray 无法自动推断。解决先decode_timesFalse打开手动查看ds[time].attrs然后自己构造时间索引。或者用xr.decode_cf(ds, decode_timesTrue)配合use_cftimeTrue参数让 xarray 用 cftime 库处理非标准日历。5.4 内存溢出大文件一次性读入现象打开一个几 GB 的再分析文件ds[t2m].values直接把内存吃满进程被系统杀掉。原因.values会把整个数组加载到内存如果文件里有几十年的逐日数据轻松超过 10 GB。解决用 xarray 的惰性加载和分块。打开时加chunks{time: 100}然后按时间片循环处理每次只取一小段。或者用ds[t2m].sel(time2020).values先筛选再取值减少实际加载量。# 分块读取每次处理 100 个时次 ds xr.open_dataset(large_file.nc, chunks{time: 100}) result [] for i in range(0, len(ds[time]), 100): chunk ds[t2m].isel(timeslice(i, i100)).mean(dim[lat, lon]) result.append(chunk.values)5.5 绘图中文乱码方框代替汉字现象图上的中文标签全变成方框。原因matplotlib 默认字体不含中文字形。解决在脚本开头指定中文字体Windows 用SimHeiLinux 用WenQuanYi Micro HeimacOS 用Arial Unicode MS。同时把axes.unicode_minus设为 False否则负号也会显示异常。plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False6. 进阶技巧把实习脚本改造成可复用的分析流程实习材料里的脚本通常是「一次性」的改个数据路径就得动代码。真正干活时我会把它重构成配置驱动的小流程把路径、变量名、时间范围、统计方法全抽到一个 YAML 或 JSON 配置里主脚本只负责读配置、调度函数。这样换一套数据只需要改配置文件不用碰代码。import yaml # config.yaml 内容示例 # data_path: data/reanalysis.nc # variable: t2m # time_range: [2020-01-01, 2020-12-31] # stat_method: mean with open(config.yaml, r, encodingutf-8) as f: cfg yaml.safe_load(f) ds xr.open_dataset(cfg[data_path]) var ds[cfg[variable]].sel(timeslice(*cfg[time_range])) if cfg[stat_method] mean: result var.mean(dim[lat, lon]) elif cfg[stat_method] trend: # 调用前面写的趋势计算函数 result compute_trend(var)另一个实用技巧是加日志。实习脚本出错了只能靠 print 猜改成 logging 模块后每个步骤的耗时、数据形状、缺测比例都能记录下来排查问题时不用重新跑一遍。import logging logging.basicConfig(levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s) logging.info(f数据形状: {var.shape}, 缺测比例: {np.isnan(var.values).mean():.2%})还有一个验证习惯每次算完统计量随手用另一个独立方法交叉验证。比如用 xarray 的polyfit算趋势再用 scipy 的linregress算一遍两个结果对不上就说明中间某步出了问题。这个习惯帮我逮住过好几次维度对齐错误——xarray 广播时如果坐标没对齐它会自动做笛卡尔积结果形状悄悄变大不检查根本发现不了。从那以后我每次拿到新的气象数据都强制走一遍「打印结构 → 检查坐标顺序 → 验证缺测编码 → 小样本试算」的流程确认无误再跑全量。希望帮到你。本文还有配套的精品资源点击获取