简介PyMICAPS 是一款面向气象数据处理与可视化的开源 Python 工具源码包适合气象科研人员、预报员及 Python 开发者。它支持 GRIB、NetCDF、ASCII 等多种气象数据格式提供插值、统计、时间序列分析等预处理功能并内置等压线图、风场图、散点图等常用图表还能完成地图投影转换与动态动画展示。压缩包共 95 个文件大小约 6.74MB主要包含 py 源码如 MicapsData、Contour、Projection、Map 等模块、txt 示例数据与说明、png 效果图以及 shp/shx/dbf 地图边界数据、yml/xml 配置和 zbak 备份文件目录结构清晰便于按需查阅。目前已有 97 人学习浏览。该资源附带完整源码、示例数据、说明文档和边界地图用户可直接运行演示或基于现有模块扩展自定义功能用于理解 MICAPS 数据解析流程、学习气象可视化开发也可作为课程设计或科研项目的实用工具箱。1. 气象数据里绕不开的 MICAPS该由 Python 来处理了MICAPS 在国内气象业务里几乎是无处不在的存在模式输出、实况拼图、预报订正太多系统最终都以这一类格式落地。PyMICAPS 的思路很直接——把读 MICAPS、写 MICAPS、再把网格场画成天气图这套动作用 Python 的方式完整重做一遍让没有 Micaps 客户端环境的开发者也能在自己机器上处理这批数据。它真正解决的痛点不是“画图难”而是格式杂站点数据、格点数据、文本头、二进制体、缺测值、经纬度方向这些细节叠加起来足以让一个数据分析项目卡住一整天。适合有 Python 基础、正在做气象服务或数据平台的工程师也适合想把 MICAPS 数据接入数据分析与可视化链路的研究生和运维开发。上手门槛不高恰好也是做开源文档贡献和二次封装的好题材。2. PyMICAPS 的数据对象先分清站点和格点再谈读写2.1 第 3 类站点数据与第 4 类格点数据的解码差异MICAPS 的格式族按 diamond 类型区分。日常最常见的是第 3 类和第 4 类第 3 类是站点观测、探空和离散点数据文件里一行对应一个站的若干要素第 4 类是均匀格点场可以来自模式输出也可以是插值后的实况分析。PyMICAPS 在解码时的约定并不复杂——站点数据被组织成记录列表每条记录带站号、经纬度、高度和要素数组格点数据被组织成一个带描述信息的二维数组头部里存起始经纬度、经纬距和格点数。两类数据在后续处理里行为完全不同站点数据要先做插值才能画等值线格点数据可以直接进入 Matplotlib 或 Cartopy。这也是我拿到文件先问一句“这是 diamond 几”的原因。按格点去读站点文件得到的场必然乱序反过来按站点去读格点结果更没法用。PyMICAPS 的处理思路与此一致先识别文件头里的 diamond 类型再给对象挂上对应属性和方法调用方不用自己维护一堆格式分支。2.2 读取 diamond 4 格点场的最小示例常见的做法是先把整个文件读成字节流或文本行头部参数解析出来数据体交给 NumPy 组织成二维数组。下面是我在本地跑通的最小代码import numpy as np from pymicaps import read_grid # 常见实现里提供统一入口 # 读入一个 diamond 4 文件 grid read_grid(Z_SURF_C_BABJ_20240101120000_P_RT_000.TXT) # 数据体和网格信息分离 data grid.data # 二维 ndarrayshape 为 (ny, nx) lat0, lon0 grid.lat0, grid.lon0 # 起始纬度、经度 dlat, dlon grid.dlat, grid.dlon # 纬向、经向格距 nx, ny grid.nx, grid.ny # 经向格点数、纬向格点数 missing grid.missing_value # 缺测标记常见是 9999read_grid之后data是纯数值数组缺测点保持为文件里给的缺测值。拿到四个网格参数已经能把经纬度坐标铺出来了。这一步关键不是代码多炫而是确认坐标原点和步长是否与数据描述一致。有些文件 dlat 是负值表示纬度从北往南递减不处理的话后面绘图时纬度轴会倒挂。2.3 为什么不用 GRIB 或 NetCDF 替代处理气象数据的人常问既然有 GRIB 和 NetCDF为什么还要专门写一个 PyMICAPS因为历史积压和业务交换链路里MICAPS 仍然大量存在。GRIB 长于压缩和标准化编码NetCDF 适合多维科学计算而 MICAPS 胜在简单——一个文件就是一个场文本格式可以直接看。PyMICAPS 的价值是把这种简单格式重新接回现代 Python 生态而不是要求业务方把全链路都改成 NetCDF。格式组织结构典型用途处理建议diamond 1/2站点填图地面观测、高空风站点列表需插值diamond 3离散站点要素降水、温度实况站号关联空间插值diamond 4均匀格点场模式输出、分析场网格对象直接绘图GRIB/NetCDF多维数组科研、再分析统一转成 xarray 再对接写回只在两个场景下做一是给下游保留 MICAPS 格式的老系统投喂数据二是临时生成一个能被 Micaps 客户端打开的文件做交叉验证。写回时头部说明里的时次、时效、层次必须原样带回否则客户端打开后产品类型是错的。手写容易漏字段尽量用库或工具函数完成。2.4 头部参数解析的常见坑diamond 4 的头部信息各版本差异不小。早年很多文件是 4 行说明加数据后来部分业务文件把时次说明合并成两行。解析时不能按固定行号硬编码最好按行首关键字识别。我见过不止一次因为第 3 行是多余的备注字段导致程序把备注里的数字当成起始经度整个场偏移了几十公里。PyMICAPS 这类库一般会做兼容处理但如果你在改造旧脚本务必先打印前 5 行确认结构。3. 把 MICAPS 网格变成可用数据坐标生成、裁剪与插值3.1 经纬度网格生成与方向检查读进来的二维数组本身不携带经纬度坐标坐标是用头部参数算出来的。这里我习惯先做一个方向检查# 用起始经纬度、格距和格点数构造一维坐标 lons lon0 np.arange(nx) * dlon lats lat0 np.arange(ny) * dlat # 方向校验多数文件纬度从南到北经度从西到东 if lats[-1] lats[0]: data data[::-1, :] lats lats[::-1]这个检查很廉价但能避免后面大部分绘图问题。部分模式输出的 MICAPS 数据格点从北往南存不翻转的话画出来的等值线是倒置的叠加到地图上时偏差可以达到数百公里。除了看首尾纬度还可以拿一个已知城市的经纬度去数组里反查数值与站点实测对比。方向正确这件事越早确认越好等到出图再发现排查成本高得多。3.2 按区域裁剪与缺测处理实际项目中经常只关心某省或某流域。裁剪时不能直接对数组做切片因为切片保留的是格点序号不是经纬度范围。更稳的办法是构造经纬度掩膜再取数# 示例长江中下游区域 [28N, 34N, 108E, 123E] lat_min, lat_max 28.0, 34.0 lon_min, lon_max 108.0, 123.0 mask ( (lats[:, None] lat_min) (lats[:, None] lat_max) (lons[None, :] lon_min) (lons[None, :] lon_max) ) sub_data np.where(mask, data, np.nan)掩膜法的好处是 shape 不变后续叠加边界、画等值线时坐标系完全一致。不少人在这一步直接做数组切片一旦数据源换了分辨率或起始经纬度区域就漂移了。缺测值也要在裁剪前统一转成np.nan否则等值线会把 9999 当作真实数值参与计算出现整片异常大值区。3.3 站点数据插值与分辨率对齐当数据是 diamond 3 时画等值线必须先做空间插值。PyMICAPS 通常依赖 SciPy我常用griddata的 cubic 或 linear 方法。温度场这类平滑要素用 cubic降水用 linear 更好——降水空间变化剧烈且大量为零值cubic 容易插出负值或虚假波纹linear 至少不会发明数值。from scipy.interpolate import griddata # stations: (n, 2) 的经纬度数组; values: (n,) 的要素值 grid_lons, grid_lats np.meshgrid(lons, lats) interp_data griddata( stations, values, (grid_lons, grid_lats), methodlinear, fill_valuenp.nan )场景推荐方法说明温度、高度场cubic平滑适合等值线降水、阵风linear避免负值和虚假极值快速预览nearest计算快但锯齿明显插值完成后建议做一次范围校验结果最小值不应明显低于原始站点最小值。出现异常时要先检查站点经纬度是否把纬度和经度顺序传反。这类错误在交互式脚本里最容易藏起来因为图能画出来只是整体偏移。3.4 与 NumPy 广播相关的隐藏坑mask 构造里用了lats[:, None]和lons[None, :]这是为了让纬度做行方向扩展、经度做列方向扩展。如果漏掉维度扩展直接比较结果会变成逐元素比较或者抛出形状不匹配。另一个相关问题是网格方向不一致lat 从大到小、lon 从小到大np.meshgrid默认的indexingxy在这种组合下会产生转置的坐标矩阵。建议统一用np.meshgrid(lons, lats, indexingij)让返回的第一个数组形状与 data 的行列一一对应后续画图填色不用再担心转置问题。4. 用 Matplotlib 与 Cartopy 绘制天气图等值线、填色与风羽4.1 不依赖 Micaps 客户端的产品图生成路径PyMICAPS 的绘图能力解决的场景是业务服务器上没有图形界面没有 Micaps 客户端也要每天按时生成产品图。绘图链路通常是解码网格、叠加地图、画等值线和填色、保存 PNG。下面给出一个最小可用的骨架。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np # 假设 data 是二维格点lons/lats 是二维网格坐标已经对齐 fig plt.figure(figsize(12, 9)) ax plt.axes(projectionccrs.PlateCarree())直接在投影坐标系里画图的要点是所有经纬度数组都作为数据坐标传入Cartopy 负责把它们投影到地图上。不要手动把经纬度转成米制坐标那会引入不必要的误差也让代码更难维护。4.2 等值线与填色的参数实践# 填色层 cf ax.contourf( lons, lats, data, levelsnp.arange(0, 51, 5), # 以 0~50 为例 cmapRdYlBu_r, alpha0.85, extendboth ) # 等值线层 cs ax.contour( lons, lats, data, levelsnp.arange(0, 51, 5), colorsblack, linewidths0.6 ) ax.clabel(cs, inlineTrue, fontsize8, fmt%.0f) # 地图要素与范围 ax.coastlines(resolution50m) ax.add_feature(cfeature.BORDERS, linewidth0.4) ax.set_extent([105, 125, 25, 40], crsccrs.PlateCarree()) plt.colorbar(cf, axax, pad0.02, labelmm) plt.savefig(pymicaps_grid.png, dpi200, bbox_inchestight)等值线层级用np.arange生成确保填色和线层的 levels 完全一致否则图例与色标对不上。extendboth处理超出色标范围的数值降水产品里零值密集时尤其有用。颜色表的选择上温度常用RdYlBu_r降水常用Blues。如果你要把产品接到可视化大屏或 Web 页面上保存时加transparentTrue生成透明背景 PNG叠加到深色底图上更协调。4.3 站点数据画风羽与落区diamond 2 高空风数据画风羽用 Matplotlib 的barbsdiamond 3 站点要素画填点图时要注意经纬度顺序。barbs的输入是 u/v 分量数组和经纬度网格同形状。如果站点是散点可以先插值到网格再画也可以直接用站点经纬度画不插值的填色散点图。# 站点风场示例lon_s, lat_s 是站点经纬度u, v 是风分量 ax.barbs( lon_s, lat_s, u, v, length5, barb_incrementsdict(half2, full4, flag20) )barb_increments控制风羽的刻度比例。地面风通常 half2、full4、flag20高空风如果需要更细的分辨率可以改成 half1、full2、flag10。这个参数不声明时 Matplotlib 有默认值但默认值不一定符合气象业务习惯画出来很容易被预报员质疑。4.4 缺测区与零值区的显示控制降水产品画等值线时零值区会被无数条线叠满。常见做法是先做一次数据掩膜把小于 0.1 的值置为np.nan再传给contourf。这样零值区留白线和色标都集中在有效降水区。缺测区域的处理同理插值后的边缘区域往往是大片 nan直接画会出现色标空洞。可以先用np.ma.masked_invalid包装数据让 Matplotlib 自动跳过无效区域同时把add_feature的陆地和湖泊要素叠在数据层下面保证边界清晰。5. 从单文件脚本到批处理流水线格式封装与结果校验5.1 把解码结果封装成 NetCDF避免重复解析MICAPS 文件是文本存储解析一次的开销在单个文件上可以忽略但几百个时次累积起来就不一样了。我一般会把当天的 MICAPS 集合先解码然后用xarray统一封装import xarray as xr ds xr.Dataset( data_vars{temp: ((time, lat, lon), temp_stack)}, coords{ time: time_list, lat: lats, lon: lons, }, ) ds.to_netcdf(preprocessed_20240101.nc)temp_stack的形状是(ntime, nlat, nlon)这样后续做时间序列分析、区域平均甚至训练模型时都不需要再碰 MICAPS 原始文件。封装 NetCDF 时记得把缺测值写成encoding属性统一为np.float32可以压缩体积同时保留缺测语义。5.2 解码正确性的交叉验证方法最可靠的验证不是看单张图而是取一个已知站点的经纬度在解码后的格点场里反查最近格点的值与站点实测或官方产品对比# 以武汉站为例30.5N, 114.3E target (114.3, 30.5) jd np.argmin(abs(lons - target[0])) wd np.argmin(abs(lats - target[1])) print(nearest grid value:, data[wd, jd])这个输出如果和站点观测差在合理范围内说明坐标方向、格距和读入逻辑整体正确如果差异极大优先检查经纬度方向是否需要翻转。另一种做法是把read_grid读出的头部信息与文件前几行人工核对确认起始经纬度和格点数量没有解析错位。最后用np.nanmax和np.nanmin检查全场的数值范围模式降水场出现负的极大值往往意味着缺测没有处理干净。以上三步跑完再进入批处理循环就不会把有问题的数据一路带到产品图里。本文还有配套的精品资源点击获取
