简介一份面向气象数据分析初学者与Python学习者的轻量实操示例包聚焦如何用Python完成气象数据的读取、清洗、分析与可视化展示覆盖从CSV站点数据导入到温度等要素变化图绘制的完整链路。包内共15个文件包括10个气象观测CSV样本、4个Python脚本及1个说明文档压缩包整体仅14KB结构小巧、便于快速上手。已有9028人学习下载适合想通过短小代码了解pandas、matplotlib、seaborn等库在气象场景中实际用法的读者。脚本按不同城市分组组织数据示例中既包含用pandas读取站点记录并绘制时间变化曲线也涉及scipy时间序列分析与周期性观察可通过对照多个城市样本理解地理差异对气象要素的影响。项目目录清晰代码与数据一一对应运行后即得可视化结果边看边改即可迁移到自己的数据集是一份入门气象数据分析与可视化流程的实用参考。1. 气象数据分析与可视化先从NetCDF到一张能用的天气图接到一份气象数据分析与可视化的活儿甲方给的是全国逐日气温NetCDF文件要求出三张图今年与去年的距平对比、夏季热浪过程的温度演变、冬季冷空气南下的空间分布。很多人拿到数据第一步就打开matplotlib画图结果是坐标轴乱掉、缺测值把色条拉爆、时间序列里混进多年数据。这篇文章讲的就是如何用Python把这套气象数据分析与可视化流程走通从读取NetCDF原始数据到清洗缺测、校准时间坐标再到计算距平和区域平均最后用Cartopy画带底图的空间场用matplotlib动画呈现时间演变让图表能放进报告或大屏。适合刚把python入门过的分析岗新人也适合气象、农业、环境行业想自建分析流程的从业者。2. 环境与数据准备装好Python库之后先读懂数据再动手做气象分析最先要过的关不是算法而是数据读取。气象行业的原始数据大多是NetCDF格式一个文件里同时存在时间、纬度、经度、气压层等多个维度和平时处理的表格结构完全不是一回事。这一章先把环境怎么装讲清楚再规范读取流程最后把数据清洗的几个边界条件交代完。2.1 依赖环境为什么是xarray加netCDF4加cartopy我常用的环境组合是Anaconda基础环境加三个核心库xarray负责多维数组和NetCDF文件的读写底层驱动是netCDF4库绘图部分用matplotlib加cartopy。之所以不推荐直接用netCDF4库读数据是因为它返回的是底层数组对象手动管理维度坐标非常麻烦xarray把这层包了一层读进来就是带维度名的DataArray后续按时间、按区域切数据都在明面上。pip install xarray netCDF4 matplotlib cartopy在Windows上建议顺手把cftime也装上后面处理带hours since 1900-01-01这种参考时间格式的坐标会用到pip install cftime装完后在python交互式环境里验证一下各库版本能否配合。实际工作中遇到过xarray升级到新版后和高版本cartopy在投影计算上的兼容性问题表现是运行不报错但出图时地图边界是空的这种黑匣子问题最耗时间。建议用一个固定版本的虚拟环境跑完整套流程别频繁升级。cftime这个库单独说明一下CF约定时间坐标是气象数据的标准写法但pandas原生不认识它需要cftime做转换层后面在避坑章节会展开讲。2.2 NetCDF读取三步维度、变量、属性逐项确认拿到任何一份气象数据文件不要急着转DataFrame也不要急着画图。先做三个确认动作看维度、看变量、看全局属性。这三个动作能省掉后面90%的排错时间尤其是当你接手的是别人传过来的数据集时。import xarray as xr ds xr.open_dataset(temp_daily_2023.nc) print(ds.dims) # 维度信息time/lat/lon分别多长 print(ds.data_vars) # 变量列表具体有哪些物理量 print(ds.attrs) # 全局属性数据来源、单位、参考时间这段代码的逻辑是先把数据集对象打开然后分别打印三部分元数据。dims告诉你每个维度的长度比如time是365、lat是181、lon是360说明这是一份全球1度分辨率的逐日数据data_vars列出变量常见的有temp、precip、pressure也可能叫t2m、tp这类短名需要查属性确认attrs里最关键的是单位说明气温是开尔文还是摄氏度降水是毫米还是千克每平方米差一个数量级图就全错了。参数说明上xr.open_dataset的第一个参数是文件路径还可以加engine参数显式指定读取引擎比如enginenetcdf4。如果文件特别大比如几十GB的再分析资料建议加上chunks参数做懒加载这一步对机器的内存压力会小很多。ds xr.open_dataset(temp_daily_2023.nc, chunks{time: 30})加了chunks之后文件不会一次性读进内存而是按time维度每30天一个块真正计算时才去读。处理ERA5这类再分析数据时这个参数几乎是必加的否则16GB内存的机器根本扛不住全球逐小时数据。2.3 数据清洗边界缺测值、时间戳和经纬度单位读取流程走完后下一步是清洗。气象数据的脏和业务数据的脏不太一样常见的有三种情况。第一种是缺测值占位缺测值通常用9999、1.0e36这类极端值填充如果没处理就参与统计均值和方差会被拉到一个离谱的数值上。第二种是时间戳格式NetCDF里的time维度经常是hours since 1900-01-01 00:00:00这种参考时间格式需要转成datetime类型才能做按月聚合。第三种是经纬度坐标单位有些数据集的lat/lon单位是弧度而不是度不转换直接切片会切出不存在的位置。# 处理缺测值占位 ds ds.where(ds[temp] ! -9999) # 统一时间坐标为datetime类型 ds xr.decode_cf(ds) # 确认经纬度范围 print(float(ds.lat.min()), float(ds.lat.max()))处理缺测值时where会把等于9999的位置置为NaN后续计算默认会跳过NaN这比手动填0要安全得多。decode_cf是把CF约定的时间编码转成可读的datetime转完后可以直接用ds[temp].sel(time2023-01-15)这种写法取数。经纬度范围打印一下如果min和max是0.5到1.5这种量级说明单位是弧度需要乘以180再除以圆周率换算回度。清洗这一步看起来琐碎但它决定了后面所有图会不会翻车。常见的情况是数据源说明里写着温度单位是K实际文件里存的是摄氏度画出来的全国7月平均气温只有5度色标全蓝错得离谱。所以每次拿到新数据集我习惯先打印一个具体点的数值和常识对一下数量级再往下走。3. 核心分析气温、降水、气压的时间序列与空间场计算环境就绪、数据干净之后才进入真正意义上的气象分析。这一章围绕三个最常见诉求展开时间序列上算月平均和距平空间场上做区域平均和季节平均场以及把极端事件量化成简单指数。每一步都给出最小可复现代码并说明参数调整的边界。3.1 时间序列用resample算月平均用groupby算气候态距平温度序列分析最基础的两个操作是降采样和距平计算。降采样把逐日数据聚合成逐月减少噪音距平表示某个时刻相对气候态的偏移是判断偏暖还是偏冷的核心指标python数据分析与可视化里这块也是最常用的。import xarray as xr ds xr.open_dataset(temp_daily_2023.nc) temp ds[temp].where(ds[temp] ! -9999) # 逐日数据降采样为逐月平均 monthly_mean temp.resample(time1ME).mean() # 计算逐月气候态假设文件内含多年数据按月份分组平均 climatology temp.groupby(time.month).mean(dimtime) # 距平 当月值 - 气候态 anomaly monthly_mean - climatology理解这段代码的关键在于resample和groupby两个操作的分工。resample是按时间频率重采样1ME表示按月取末值作为分组锚点常用写法还有1MS表示按月初。mean是聚合函数表示组内平均也可以用max、min换成月最高温、月最低温。groupby更灵活它按时间坐标的month分量分组把多年数据中所有1月聚到一起求平均得到的就是1月气候态dimtime明确表示缩减时间维度只保留month作为新的分组标签。参数上要注意的是time.month返回的是1到12的整数所以climatology的维度是month长度12。如果你手里的数据只有单年算不了气候态可以退一步用多年历史数据单独算一个climatology.nc文件再和当年数据做对齐运算。这个对齐需要保证两个文件网格分辨率完全一致否则xarray会按坐标自动插值这在气象分析里常常是隐患最好在数据准备阶段就把分辨率统一。距平的工程化价值在于它消除了季节周期。对比1月和7月的原始温度没有意义但对比它们的距平有物理含义。做可视化时距平图用的色标通常是对称的蓝红两色这种图的视觉冲击力比原始温度图强很多也是气象报告里最常出现的图表类型。3.2 空间场计算区域平均、季节平均场和多层差值当分析目标是整个华北平原今夏平均气温比常年高了多少就不能只看单点需要做空间聚合。空间聚合的常见做法是先按经纬度切片选出区域再对lat和lon两个维度求平均。# 选区域华北平原大致在35—40N, 110—120E region temp.sel(latslice(35, 40), lonslice(110, 120)) # 区域时间序列对空间维求平均 region_series region.mean(dim[lat, lon]) # 多年同一季节的空间平均场 seasonal_mean temp.sel(timetemp[time.season] JJA).mean(dimtime)sel方法用slice做切片时要求坐标轴是单调递增或递减的如果lat坐标在数据里是南纬在下、北纬在上的递增排列直接切片即可如果是几十年前数据的北纬在上排列需要先sortby处理。slice(35,40)表示取35到40度之间的所有格点注意这个区间是闭区间边界值都会被包含。mean(dim[lat,lon])的效果是把二维空间压成一个点得到一条纯粹的时间序列适合做后续的趋势回归。季节平均场这里的写法利用了season这个派生坐标xarray会自动从标准时间坐标生成季节标签DJF是冬季MAM是春季JJA是夏季SON是秋季。如果需要限定夏季三个月也可以写得更明确summer_months temp.sel(timetemp[time.month].isin([6, 7, 8]))这段代码等价于先筛选月份再求平均逻辑上更直白适合新手理解。空间场分析还有一个常用操作是沿气压层差值比如分析500hPa和850hPa的高度差来推断冷暖平流做法是用interp对level维度线性插值interp_level temp.interp(level[500, 850])interp是xarray内置的线性插值方法目标level序列放在列表里插值后对象的维度会变得和原始level维度一致。要注意的是如果原始数据里level坐标不是等间距的线性插值的误差会放大这时应该考虑用scipy的带权插值函数做更高阶处理。这一步在无资料地区的气象分析里经常单独拿出来评估。3.3 简单气候指数把极端高温/低温事件量化出来分析做完常规平均还有一个高频需求是把极端事件量化。气候指数本质上是对原始序列做阈值判断再统计频次或累计量。这里以高温日数和热浪过程为例这两个指标在夏季高温评估和农业灾害风险里非常常用。import numpy as np # 定义高温日日最高温超过35℃数据为开尔文时需换算 threshold 308.15 # 35℃ 308.15K hot_flag (temp threshold).astype(int) # 连续高温日数用累计计数识别连片时段 run_length hot_flag.groupby( (hot_flag ! hot_flag.shift(time1)).cumsum() ).cumsum()高温日数的定义极其简单布尔比较加类型转换就够。关键点是单位换算数据是开尔文35℃对应308.15K换算错了所有指数都会错得离谱。连续高温日数的计算稍微绕一点先构造一个01标志序列再用shift判断前后两天是否连续不连续则分组计数加1最后用cumsum输出连续长度。这个逻辑是气象统计里识别过程事件的通用做法同样适用于连续降水日数、连续无降水日数。实际业务中做这类时间序列分析最容易出问题的不是算法而是边界条件数据第一天的前一天是什么状态我们并不知道。所以计算run_length时第一天的连续值天然不可信使用时要考虑截掉序列首尾各一个点或用多年数据把边界补齐。这个坑在后面避坑章节里会再提一遍。4. 可视化落地静态天气图到动态气象面板分析结果是数字可视化才让它变成能对外沟通的图。这一章按三个层次来写先用cartopy画一张带底图的静态空间分布图再用matplotlib动画把时间演变串起来最后用Panel搭一个可拖拽、可切换参数的交互看板。可视化部分的坑往往不在画图本身而在投影和坐标转换上这两个问题会在第5章统一排查。4.1 用Cartopy画一个带底图的温度空间分布图单张静态图是最基础的产出一张合格的空间分布图需要包含四个要素地图底图、等值线填色、色标、投影信息。cartopy负责前三者的地理框架matplotlib负责渲染细节。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature temp_jan ds[temp].sel(time2023-01-15) fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:) im temp_jan.plot.pcolormesh( axax, transformccrs.PlateCarree(), cmapRdYlBu_r, cbar_kwargs{label: 2m temperature (K)} ) ax.set_extent([70, 140, 15, 55]) ax.set_title(2023-01-15 Surface Temperature) plt.savefig(temp_2023_jan15.png, dpi300, bbox_inchestight)这里最关键的是要理解投影坐标和数据坐标是两套系统。PlateCarree是最简单的等距圆柱投影经纬度直接映射到平面坐标对于中国区域的1度网格足够用。transformccrs.PlateCarree()告诉matplotlib数据本身是经纬度坐标画布ax的投影是另一个投影时cartopy会自动做转换。pcolormesh比contourf更适合格点数据它不做插值直接按网格画色块速度也快。cmap用RdYlBu_r这个色标的红蓝方向是暖色表示高温冷色表示低温和气象行业的读图习惯一致。set_extent控制地图显示范围这里的经纬度范围正好覆盖中国全境。savefig的dpi建议300报文章或者报告里的图这个清晰度才够。如果数据是规则1度经纬度网格pcolormesh可以直接用如果投影比较复杂比如极地立体投影数据坐标转换后会出现空白区域处理办法在第5章的避坑里会讲到。4.2 制作逐日温度演变动画matplotlib.animation的正确用法静态图只能看一个时刻逐日动画才能看到过程的演进。制作方法是用matplotlib.animation的FuncAnimation逐帧更新数据把每一帧存成PNG再合成视频。实际项目中我常用两种产出方式直接保存mp4或者生成HTML5的动画嵌入网页里后者在演示时非常方便。from matplotlib.animation import FuncAnimation temp ds[temp].where(ds[temp] ! -9999) fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.5) def update(frame): ax.clear() ax.add_feature(cfeature.COASTLINE, linewidth0.5) data temp.isel(timeframe) im data.plot.pcolormesh(axax, transformccrs.PlateCarree(), cmapRdYlBu_r, vmin250, vmax310, add_colorbarFalse) ax.set_title(fTime: {str(data.time.values)[:10]}) return im anim FuncAnimation(fig, update, frameslen(temp.time), interval200) anim.save(temp_2023_animation.gif, writerpillow, dpi120)FuncAnimation的核心是update函数这个函数在每一帧执行一次先clear掉上一次的绘图再重新绘制新一帧的数据。这里的注意点有两个。一是vmin和vmax手动固定不要让它自动缩放否则每帧的色标范围都在变动画看起来会闪烁二是add_colorbarFalse动画里每一帧都加色标会非常慢颜色图例统一放在最后单独输出或固定。interval200表示每帧间隔200毫秒即每秒5帧。保存gif用pillow writer保存mp4则用ffmpeg writer需要系统装有ffmpeg。如果你只需要在Jupyter里看效果不用保存文件update函数返回图像对象后用HTML(anim.to_html5_video())渲染即可。还有一个性能优化点值得注意——如果数据维度很大每帧都重新绘制整个网格很慢可以重用上一次的quadmesh对象只在set_array里替换数据速度能快3倍以上。4.3 用Panel搭一个可拖拽的气象数据小看板走到这一步单张图和动画都能产出了。如果分析结果要交给非技术同事看或者要放到科室的大屏上轮播就需要一个交互式的可视化小面板。这里我会用Panel框架它和Jupyter绑定得好不需要额外写前端代码就能搭出可拖拽、可切换参数的看板。import panel as pn pn.extension() ds xr.open_dataset(temp_daily_2023.nc) selector pn.widgets.Select(name观测日期, optionslist(ds.time.values)[:10]) pn.depends(selector.param.value) def plot_map(date): data ds[temp].sel(timedate) fig plt.figure(figsize(8, 5)) ax plt.axes(projectionccrs.PlateCarree()) data.plot.pcolormesh(axax, transformccrs.PlateCarree(), cmapRdYlBu_r) plt.close(fig) return fig layout pn.Column(selector, plot_map) layout.servable()用Panel搭看板的思路是用一个widget组件绑定日期pn.depends装饰器让plot_map函数在参数变化时自动重算pn.Column组合成纵向布局。servable方法可以让这个面板通过panel serve script.py直接起一个本地服务浏览器打开就能交互。做交互看板时要注意一个性能问题每次拖动日期下拉框都会重新运行plot_map如果底层数据是懒加载的大型NetCDF建议把需要展示的变量先降采样到1度或2.5度网格或者用缓存机制存住最近计算过的结果。我还习惯把色标范围固定理由和动画完全一样用户拖动时颜色跳动会让人误以为数据有问题。最后这个面板里如果要展示多个气象变量可以把Select换成多选或RadioButton然后用循环生成多个输出。项目的落地形态就是一个小型气象可视化工具它的边界在于交互逻辑不宜做得很复杂真正复杂的分析计算写在Panel外面。5. 避坑气象数据处理的5个常见翻车点这一章把前面各环节里最容易让初学者翻车的五个场景单独拎出来按现象、原因、解决的顺序梳理清楚。这些问题我基本都在实际项目中踩过有些甚至踩了两遍才总结出规律。5.1 cftime和pandas时间索引打架现象读取NetCDF后执行resample或sel(time2023-01-15)报错提示time坐标不是datetime类型或者index类型不匹配。原因CF约定的时间坐标是hours since 1900-01-01 00:00:00xarray解码后生成的是cftime.DatetimeNoLeap类型不是标准的pandas.DatetimeIndex。当数据里包含非标准日历如无闰年日历时xarray不会自动转成pandas类型直接用字符串切片就会失败。解决读取时显式使用xr.decode_cf(ds)如果日历仍是非标准类型就把time坐标手动转换后写回。转换用ds[time] ds.indexes[time].to_datetimeindex()。如果数据量太大建议在清洗阶段一次性转换并保存成一个新的nc文件后续所有分析都用这个清洗后的文件避免每个脚本都要重复处理。5.2 经纬度维度名写反地图直接变形现象画出来的图经纬度范围完全不对中国区域的地图变成横向拉伸的长条形或者色块和海岸线对不上。原因不同数据源的维度命名不一致有的是lat/lon有的是latitude/longitude还有的是y/x而且数组的行列顺序可能是经纬度反的。用plot.pcolormesh时x和y参数如果传反了图自然就歪。解决先打印ds.dims确认维度名和排列顺序然后在取数时用sel方法按名字切片而非按位置取数比如sel(latitudeslice(15,55), longitudeslice(70,140))。如果数据本身是y/x这种网格坐标先重命名维度再处理ds ds.rename({y: lat, x: lon})。养成按维度名操作的习惯后基本不会再出现维度写反的翻车。5.3 缺测值混进统计计算现象区域平均结果出现NaN或者距平图上一大片区域变成空白而不是正常的色块。原因缺测值在数据里表现为特殊的_FillValue例如9999如果不做掩膜处理NaN会通过均值传播到所有下游计算。解决在进入任何统计计算之前先统一处理缺测值。用ds[temp] ds[temp].where(ds[temp] ! -9999)把占位值替换成NaN之后所有聚合函数默认skipnaTrue自动跳过NaN。需要特别注意的地方是进行空间平均时如果某个格点在整段时间都缺测平均结果仍会是NaN这时可以用region.mean(dim[lat,lon], skipnaTrue)显式跳过。更稳妥的做法是统计每个格点的有效观测数把低于阈值的格点直接标为NaN。5.4 Cartopy投影边界撕裂现象使用非PlateCarree投影如Albers等积投影或Lambert投影绘图时中国东北或南海区域的海岸线出现横线撕裂或色块边缘参差不齐。原因数据坐标是经纬度网格投影到目标投影后原本一条经线上连续的格点被切到了画布两侧或某条边界上cartopy无法自动判断格点的连续性。解决先用data.plot.pcolormesh(axax, transformccrs.PlateCarree())明确告诉绘图函数数据原始坐标是经纬度再在ax上设置目标投影。如果撕裂依然存在把数据先插值到目标投影的规则网格上再画图。插值用scipy.interpolate.griddata目标网格用投影坐标下的np.meshgrid生成。对于中国区域我还习惯把ax的边界设置成圆角矩形既好看也能掩盖某些边界处的毛刺。5.5 文件句柄未关闭内存越用越大现象运行多个分析脚本后内存占用持续上升不退出python进程就不释放甚至出现OSError: [Errno 24] Too many open files。原因xr.open_dataset默认是懒加载文件句柄不会立即释放。在循环里反复打开文件却不关闭句柄和懒加载数据会累积在内存里。尤其在读取多个NetCDF文件合成一个时间序列时最容易触发。解决养成用上下文管理器的习惯或者显式调用close。with xr.open_dataset(temp_daily_2023.nc) as ds: result ds[temp].mean(dimtime) # 文件在此处自动关闭 # 多个文件合并时合并完成后立即释放句柄 ds_list [xr.open_dataset(f) for f in file_list] combined xr.concat(ds_list, dimtime) for ds in ds_list: ds.close()这两个例子覆盖了日常使用和批处理两种场景。前者的with块在代码退出后自动关闭文件后者的循环close确保合并操作的临时句柄被释放。如果在实际运行中内存还是控制不住就要考虑用open_dataset(..., chunks{})配合dask做分布式处理而不是盲目加内存。6. 进阶把分析流程封装成可重复执行的小工具走到这一步读、洗、算、画四条链路你都跑通了。接下来值得做的一件事是把这套流程封装成一个命令行可执行的Python脚本让换数据源、换年份、换区域时不用再改一堆参数。我一般会按函数粒度拆分load_data负责读文件和清洗cal_temp_index负责时间序列和指数计算plot_map负责画图main函数负责串起整个流程。函数之间只通过DataArray对象传参不共享全局变量。这样做的最大好处是临时接到一个把去年数据重新分析一遍的需求时只需要改main里的文件路径和一个年份参数十分钟内能出全套图。def main(data_path: str, year: int, region: tuple) - None: ds load_data(data_path) monthly cal_monthly_stats(ds, year) anomaly cal_anomaly(monthly) plot_map(anomaly, region, outputfanomaly_{year}.png)这个封装是不是过度设计要看场景。如果你只是自己写一次性的分析脚本函数拆分反而增加负担但如果要让同事复用或者要定时产出一组图这个封装就是值得的。参数校验和异常处理建议在main入口做比如文件路径不存在时给出明确提示而不是抛一个看不懂的堆栈。另外一个我常用的习惯是画完图后把计算的中间结果同步导出成CSV方便非Python用户直接拿Excel做二次分析。所谓可视化大屏最后的落地也往往是这些CSV加PNG再加少量交互组件组合出来的底层的计算逻辑并不复杂。回看这套流程最让我记忆深刻的教训是气象数据分析中80%的错误不是分析本身的问题而是数据格式、单位、缺测和投影这些周边问题。先把数据边界摸清楚分析的信心才有基础。希望这篇笔记能帮你少走一遍这些弯路。本文还有配套的精品资源点击获取
