简介2020年启动的NASA IMPACTS大西洋沿岸暴风雪微观物理学与降水研究项目利用P-3B和ER-2飞机、微物理探测器、地面雷达及GPM、GOES-16等多颗卫星协同观测是东海岸暴风雪三十年来首次综合研究。这里提供的是该项目航测数据与Python可视化程序包面向气象数据爱好者、Python数据分析学习者及科研入门者。资源共28个文件以21个CSV航迹气象导航数据为核心辅以2个Python脚本、3张PNG航迹图、说明文档与许可文件整体压缩包约21.47MB。各CSV文件按航班日期组织便于按任务检索flight5.py和impacts.py可读取CSV自动绘制P-3B、ER-2单机航迹及双机联合对比图便于快速复原2020年2月5日中西部浅雪带前锋带调查等典型飞行场景文件目录含datasets和img等分组结构清晰。目前已有191人学习适合希望结合真实科研项目理解气象观测数据、数据清洗与可视化流程的读者。 把NASA IMPACTS暴风雪研究程序的原始数据变成一套完整可用的可视化图像这件事听起来硬核但在我看来它恰恰是数据分析里最有成就感的一类工作数据源足够专业图里每一条色带、每一个回波结构都能直接回答一个气象学问题而不是为了好看而好看。NASA IMPACTSInvestigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms是一个专门针对美国东海岸暴风雪云系微观物理与降水过程的研究计划飞行平台会带着机载雷达和云微物理探头直接穿进冬季风暴拿到的是非常规整的科学观测数据。这篇内容就以“nasa_impacts”这个Python可视化程序为例完整拆解如何把这些NetCDF格式的遥感数据变成可复现、可扩展、可分享的三维可视化成果。1. NASA IMPACTS到底在飞什么数据1.1 这架“会飞的实验室”采集了什么IMPACTS计划的主力观测平台以高空飞机和低空飞机配合为主机载设备基本可以分成两类。一类是主动遥感设备比如气象雷达和云雷达它们从飞机上向下或向前发射电磁波通过接收回波来反演云层内部的反射率因子、多普勒速度、谱宽等参数。另一类是原位测量设备包括各种云粒子成像仪、液态水含量探头它们负责记录飞机实际穿过的那个位置到底有什么粒子、粒子长什么样、浓度是多少。这两类数据合在一起就构成了暴风雪研究中非常关键的“同一时刻、不同尺度”的观测视角雷达看到的是云系宏观结构粒子探头看到的是微观相态细节。两者都需要一帧一帧地组织起来再沿着飞行轨迹串成时间序列。所以在设计nasa_impacts这个程序时我一开始就明确了目标不能只画一张图而是要能按飞行架次、按时间窗、按垂直高度和水平区域批量产出图形否则几十趟航班的数据靠手工在绘图软件里点来点去效率完全跟不上。1.2 为什么可视化这一步比想象中更重要很多人可能会觉得气象研究不就是为了出论文里的那张示意吗其实不是。暴风雪系统是强非均匀、快速演变的天气过程飞机在几个小时内穿过的路径相当于在风暴的不同部位“打了几针采样”。没有高效可视化你很难知道哪一段数据对应的是对流单体前缘哪一段是层状云降雪区哪一段又是飞机颠簸导致的仪器噪声。可视化在这里承担的是“数据质量检查、物理过程理解、成果沟通表达”三重职能。我见过不少刚拿到数据的人直接跑统计模型去求相关系数结果因为数据里混入了雷达旁瓣污染得出了完全离谱的结论。正确做法是先画一批反射率剖面图和水平扫描图用眼睛确认回波形态是否连续、粒子相态分布是否符合常识再进入定量分析。所以这篇文章不只是在讲一个画图程序而是在讲一个气象数据工作流里绕不开的基础设施。2. 环境与数据读取用Python把NetCDF变成可操作的数组2.1 工具链选型和理由nasa_impacts在工具选型上非常常规但组合起来相当顺手。核心依赖是这几个xarray带标签的多维数组库读取NetCDF时能把坐标、变量名、单位全部保留这是处理复杂气象数据时最有价值的一点。netCDF4底层NetCDF读写接口xarray会调用它作为后端。matplotlib出图主力配合cartopy处理地图投影和海岸线。numpy/scipy做极坐标转笛卡尔坐标、插值、滤波的基础。cartopy地图投影库画飞行轨迹和扫描图时离不开。可选装pyart专门面向气象雷达的Python库做衰减订正、坐标变换时能省很多事。环境建议直接用conda创建避免源码编译Cartopy时踩坑。实测下来这套组合在Python 3.11下兼容性不错安装命令如下conda create -n nasa_impacts python3.11 conda activate nasa_impacts conda install -c conda-forge numpy xarray netcdf4 matplotlib cartopy scipy pyart2.2 读取NetCDF时的第一个坑先打印再动手拿到一个IMPACTS数据文件第一反应不要急着画图先用xarray打开并打印全局信息和变量结构。这一步能解决后面80%的困惑。import xarray as xr ds xr.open_dataset(IMPACTS_radar_20220213.nc) print(ds) # 常见变量 # reflectivity 反射率因子单位通常是 dBZ # doppler_velocity 多普勒径向速度 # spectrum_width 谱宽 # 常见坐标 # time, range, azimuth, elevation, latitude, longitude从实际经验看NASA这类公开档案里的文件变量命名还算规范但也存在不同架次之间维度顺序不一致的情况。比如有的文件把time放在第一维有的把range放在第一维如果直接按固定下标取值换一个架次就很容易报错或画出镜像错位图。用xarray的好处在于可以靠变量名取数比如ds[reflectivity].isel(time0)代码里不依赖硬编码下标程序的健壮性会好很多。模块内部我会建议做一个“字段标准化”函数把不同文件里可能的别名统一映射成标准字段名例如将REF、refl、reflectivity统一成reflectivity再把坐标统一排序为(time, range, azimuth)。这样后续绘图函数只需要关心标准字段名不用在每个函数里再判断一次变量存在性。3. 绘制第一张反射率图看懂极坐标扫描的可视化逻辑3.1 机载雷达扫描和地面雷达的差异气象雷达有两种基本扫描方式地面业务雷达常做平面位置扫描PPI也就是固定仰角、转一圈方位角得到一个圆锥面的回波。机载雷达则灵活得多除了PPI还可以做垂直方向的位置扫描RHI沿某个方位角上下扫得到垂直剖面。在设计可视化程序时我专门为这两种模式做了分离PPI模式画水平地图RHI模式画高度剖面。更复杂的情况是飞机本身在运动所以每一帧扫描的雷达中心位置都在变化画图时必须把经纬度坐标算对而不是假设雷达位于某个固定站点。3.2 从NetCDF到第一张水平扫描图下面这段代码实现的是读取某一时刻的反射率数据并转换到地图坐标系下进行绘制import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np def plot_ppi(ds, time_idx0, varreflectivity): refl ds[var].isel(timetime_idx) lat ds[latitude].isel(timetime_idx) lon ds[longitude].isel(timetime_idx) fig plt.figure(figsize(12, 8)) ax fig.add_subplot(1, 1, 1, projectionccrs.PlateCarree()) ax.set_extent([lon.min() - 0.5, lon.max() 0.5, lat.min() - 0.5, lat.max() 0.5], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth1.0) ax.add_feature(cfeature.STATES, linewidth0.6) im ax.pcolormesh(lon, lat, refl, cmapturbo, vmin-30, vmax40, transformccrs.PlateCarree()) plt.colorbar(im, axax, labelReflectivity (dBZ), fraction0.025, pad0.04) ax.set_title(fIMPACTS Reflectivity PPI - Time {time_idx}, fontweightbold) return fig, ax这里有个关键点pcolormesh需要二维坐标数组通常NetCDF里的经纬度在机载雷达场景下就是二维的直接用即可。如果遇到一维坐标要先做np.meshgrid处理。颜色映射我用的是turbo视觉效果比jet更平滑比viridis更贴近业务雷达图习惯比較适合在公众号、PPT和论文初稿里快速展示。但如果你要把图投到某个要求严格的期刊还是建议查一下期刊对色标的具体要求有的期刊明确要求用感知均匀的色标。3.3 极坐标转笛卡尔RHI剖面图的实现垂直剖面图能直观展示雪带从低层到高层的演变。机载雷达的RHI扫描通常以飞机为中心记录不同距离和仰角上的回波。绘制时先把极坐标(距离, 仰角)转换成笛卡尔坐标(水平距离, 高度)再用pcolormesh画出。def plot_rhi(ds, az_idx0, varreflectivity, max_range_km60): rng ds[range] / 1000.0 # 转为公里 elev np.deg2rad(ds[elevation].isel(azimuthaz_idx)) az np.deg2rad(ds[azimuth].isel(azimuthaz_idx)) R, E np.meshgrid(rng, elev) x R * np.cos(E) * np.cos(az) z R * np.sin(E) ref ds[var].isel(azimuthaz_idx) fig, ax plt.subplots(figsize(14, 6)) ax.pcolormesh(x, z, ref, cmapviridis, vmin-30, vmax40) ax.invert_yaxis() ax.set_xlabel(Horizontal Distance (km)) ax.set_ylabel(Altitude (km)) ax.set_title(IMPACTS RHI vertical cross-section) return fig, ax有一点必须提醒机载雷达的range和地面雷达一样是从雷达天线出发的斜距不是水平距离。如果你的数据显示了地形回波或非常明显的底部杂波先别急着怀疑仪器先检查一下是否忘记做高度订正。飞机本身就有高度因此以海平面或飞机所在高度作为参考零点的差异会直接影响剖面图的解释。建议在程序里加一个参数reference_altitude_km默认取飞机GPS高度实际画图时再根据需求平移。4. 程序化出图飞行轨迹、剖面和系列帧的高效组织4.1 自动化输出目录结构真正实用的可视化程序不是只有一个画图函数而是能把一批文件、一批时间点自动处理成整齐的图集。我建议项目里保留如下目录结构nasa_impacts/ ├── config.yaml ├── nasa_impacts/ │ ├── load.py │ ├── plot_ppi.py │ ├── plot_rhi.py │ ├── flight_track.py │ └── make_quicklook.py ├── data/ │ └── 20220213/ ├── outputs/ │ ├── ppi/ │ ├── rhi/ │ ├── track/ │ └── overview/这样做的目的很直接数据源不乱动输出按图和类型分离等哪天需要做对比图或者做成Web服务时直接读对应目录就行不用再满硬盘翻文件。在config.yaml里可以指定数据文件列表、要绘制的变量名、色标范围、输出格式、DPI等参数。整个程序只需要一个入口命令python -m nasa_impacts.make_quicklook --config config.yaml --flight 202202134.2 飞行轨迹叠加时间着色除了雷达扫描图飞行轨迹本身也是重要可视化对象。把飞机轨迹按时间顺序连成线并按时间或高度着色能快速看出飞机是在风暴哪个高度层穿插。这里用LineCollection会比较方便from matplotlib.collections import LineCollection def plot_track(ds_track, t, lat, lon, alt, cmap_nameplasma): points np.array([lon, lat]).T.reshape(-1, 1, 2) segments np.concatenate([points[:-1], points[1:]], axis1) lc LineCollection(segments, cmapcmap_name, linewidth2) lc.set_array(alt[:-1]) fig, ax plt.subplots(figsize(10, 10), subplot_kw{projection: ccrs.PlateCarree()}) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.STATES) ax.add_collection(lc) ax.autoscale() plt.colorbar(lc, axax, labelAltitude (km)) return fig, ax飞行轨迹图是检查数据质量问题的一把好手。如果轨迹出现不连续的跳变通常说明GPS定位信号有丢失或者数据插值出了问题这种异常在统计指标里很难发现但画在图上非常扎眼。4.3 循环生成“帧序列”为了观察风暴演变我会按固定时间间隔生成多帧PPI或RHI图再组合成动画或GIF。程序里写一个循环import imageio.v2 as imageio from pathlib import Path writer imageio.get_writer(impacts_loop.gif, fps3) for t in range(0, len(ds[time]), 5): fig, ax plot_ppi(ds, time_idxt) fname foutputs/ppi/frame_{t:04d}.png fig.savefig(fname, dpi120, bbox_inchestight) plt.close(fig) writer.append_data(imageio.imread(fname)) writer.close()这里建议保存PNG帧之后再用imageio合成GIF不要边画边往内存里塞大数组否则几十分钟的飞行数据足够把16GB内存吃满。主力分析时还是保存单帧PNG更实用GIF只用于汇报和快速浏览。5. 实际工程中的坑给nasa_impacts项目的七条避坑记录这几条都是我在实际做可视化程序时踩过的有些是气象领域特有的有些是任何NetCDF处理项目都会遇到。写在这儿希望你拿到新数据时不会被同样的石头绊住。5.1 维度顺序不一致不同架次数据里变量的维度顺序不完全一致有的把time放最前面有的把range放中间。肉眼看不出来但脚本里的固定索引会直接取错。用xarray的isel和sel可以避免但前提是先打印ds确认维度名。更保险的做法是写一个normalize_dims(ds)函数把标准的坐标名重新排序后返回保证生产线稳定。5.2 单位不是看上去的单位反射率数据的单位通常是dBZ但有些文件里存的是线性反射率Z数值范围大不相同。如果直接用默认色标轻则图像饱和度过高重则完全看不出回波结构。读取数据后务必检查units属性如果值是mm^6/m^3这类线性单位先做10 * np.log10(z)转成dBZ再绘图。5.3 缺测值没有掩膜NetCDF文件里缺测点通常用NaN或一个大数表示比如-9999。pcolormesh遇到NaN会透明不画但遇到-9999会画出极值颜色。处理方式是在绘制前统一处理refl refl.where(refl -100, np.nan)与此同时也要检查经纬度数组里是否混入缺测值那会导致投影点跑到非洲去。5.4 色标选择比想象中重要我早期喜欢用jet后来发现等降水量线、论文审稿人、视觉正常的人都不太买账。现在默认用两种大范围回波用turbo学术剖面图用viridis或pyart内置的Reflectivity colormap。如果目标读者是气象业务人员NWS Reflectivity色标可能更符合他们的看图习惯这个可以把matplotlib.colors.ListedColormap配一组颜色分界点来实现。5.5 Cartopy与现有坐标不匹配如果直接在ax.pcolormesh里传入从NetCDF读出的经纬度且没有指定transformccrs.PlateCarree()图上所有数据点会落在投影坐标系的原点附近画出来就是一团乱麻。记住一条原则数据坐标用transform参数声明地图框用projection参数声明两个坐标系不同时Cartopy会自动转换但这层转换必须在绘图调用里显式写出来。5.6 大数据文件内存占用失控一整架次飞行雷达原始数据可能包含上千个时间步每个时间步又有多条扫描线。可以读入内存但不要一下把整个数据集的所有变量都载入。我通常用xr.open_dataset(..., chunks{time: 20})启用Dask分块绘图前只取所需变量和时间切片。处理完一批后立即del fig, ax和plt.close(fig)否则Matplotlib在循环里积累的Figure对象会越占越多。5.7 文件按架次分组时的时间标准化IMPACTS数据文件的命名通常带日期和架次号但不同年份的数据归档方式可能不同。项目里建议在配置文件里维护一个字典把架次名映射到文件路径再用pandas.to_datetime解析时间坐标。这样后面做时间轴统一、跨架次拼接时不会因为UTC偏移和本地时区的问题而对不上号。6. 从命令行脚本到可视化工具包的演进建议一旦把单张图画通下一步就是让整个流程更像一个“程序”而不是一坨在Jupyter Notebook里反复粘贴的代码。我在项目里做过的三个改进非常值得推荐。第一是加QuickLook模式。所谓QuickLook就是每处理一个数据文件时自动生成一张概览图左半边是飞行轨迹右半边是不同高度层的回波拼图底部再附一行时间标签和文件名。这种图不是给最终报告用的而是给数据检查用的好处是几乎不用看日志就能判断哪些架次的数据质量有问题。第二是增加配置文件。所有色标范围、变量名称、飞行日期、区域范围都放在config.yaml里避免每次改场景都需要改动代码。比如某个数据文件把经纬度变量的名字从latitude改成了lat只需要修改配置文件里的lat_var: lat不需要动任何绘图函数。第三是预留交互式前端接口。现在气象数据可视化经常被问到“能不能把数据发布到大屏上”其实完全可以先把PNG图集成到Grafana或自建Web页面里再挂上一张按时间排列的缩略图墙。实现方式也不复杂让程序每次跑完输出一个index.html把生成的图按文件列表写进一个简单网格页面浏览器打开就能看到整个架次的图集。以我个人的经验做这类NASA科研数据的可视化程序最难的不是写代码而是建立对数据的“手感”。你需要在一次次出图、一遍遍回头调整坐标轴和色标的过程中把每条回波带、每个数据缺口都变成自己能识别的信号。等你能在0.5倍速看飞行动画时直接说出“这架次北段在17公里高度有一个干层”就说明这个可视化程序已经真正长在你脑子里了。到那个阶段画图已经不是任务而是你理解风暴的一种方式了。本文还有配套的精品资源点击获取 NASA IMPACTS暴风雪数据可视化：用Python把NetCDF雷达数据变成三维气象图
简介2020年启动的NASA IMPACTS大西洋沿岸暴风雪微观物理学与降水研究项目利用P-3B和ER-2飞机、微物理探测器、地面雷达及GPM、GOES-16等多颗卫星协同观测是东海岸暴风雪三十年来首次综合研究。这里提供的是该项目航测数据与Python可视化程序包面向气象数据爱好者、Python数据分析学习者及科研入门者。资源共28个文件以21个CSV航迹气象导航数据为核心辅以2个Python脚本、3张PNG航迹图、说明文档与许可文件整体压缩包约21.47MB。各CSV文件按航班日期组织便于按任务检索flight5.py和impacts.py可读取CSV自动绘制P-3B、ER-2单机航迹及双机联合对比图便于快速复原2020年2月5日中西部浅雪带前锋带调查等典型飞行场景文件目录含datasets和img等分组结构清晰。目前已有191人学习适合希望结合真实科研项目理解气象观测数据、数据清洗与可视化流程的读者。 把NASA IMPACTS暴风雪研究程序的原始数据变成一套完整可用的可视化图像这件事听起来硬核但在我看来它恰恰是数据分析里最有成就感的一类工作数据源足够专业图里每一条色带、每一个回波结构都能直接回答一个气象学问题而不是为了好看而好看。NASA IMPACTSInvestigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms是一个专门针对美国东海岸暴风雪云系微观物理与降水过程的研究计划飞行平台会带着机载雷达和云微物理探头直接穿进冬季风暴拿到的是非常规整的科学观测数据。这篇内容就以“nasa_impacts”这个Python可视化程序为例完整拆解如何把这些NetCDF格式的遥感数据变成可复现、可扩展、可分享的三维可视化成果。1. NASA IMPACTS到底在飞什么数据1.1 这架“会飞的实验室”采集了什么IMPACTS计划的主力观测平台以高空飞机和低空飞机配合为主机载设备基本可以分成两类。一类是主动遥感设备比如气象雷达和云雷达它们从飞机上向下或向前发射电磁波通过接收回波来反演云层内部的反射率因子、多普勒速度、谱宽等参数。另一类是原位测量设备包括各种云粒子成像仪、液态水含量探头它们负责记录飞机实际穿过的那个位置到底有什么粒子、粒子长什么样、浓度是多少。这两类数据合在一起就构成了暴风雪研究中非常关键的“同一时刻、不同尺度”的观测视角雷达看到的是云系宏观结构粒子探头看到的是微观相态细节。两者都需要一帧一帧地组织起来再沿着飞行轨迹串成时间序列。所以在设计nasa_impacts这个程序时我一开始就明确了目标不能只画一张图而是要能按飞行架次、按时间窗、按垂直高度和水平区域批量产出图形否则几十趟航班的数据靠手工在绘图软件里点来点去效率完全跟不上。1.2 为什么可视化这一步比想象中更重要很多人可能会觉得气象研究不就是为了出论文里的那张示意吗其实不是。暴风雪系统是强非均匀、快速演变的天气过程飞机在几个小时内穿过的路径相当于在风暴的不同部位“打了几针采样”。没有高效可视化你很难知道哪一段数据对应的是对流单体前缘哪一段是层状云降雪区哪一段又是飞机颠簸导致的仪器噪声。可视化在这里承担的是“数据质量检查、物理过程理解、成果沟通表达”三重职能。我见过不少刚拿到数据的人直接跑统计模型去求相关系数结果因为数据里混入了雷达旁瓣污染得出了完全离谱的结论。正确做法是先画一批反射率剖面图和水平扫描图用眼睛确认回波形态是否连续、粒子相态分布是否符合常识再进入定量分析。所以这篇文章不只是在讲一个画图程序而是在讲一个气象数据工作流里绕不开的基础设施。2. 环境与数据读取用Python把NetCDF变成可操作的数组2.1 工具链选型和理由nasa_impacts在工具选型上非常常规但组合起来相当顺手。核心依赖是这几个xarray带标签的多维数组库读取NetCDF时能把坐标、变量名、单位全部保留这是处理复杂气象数据时最有价值的一点。netCDF4底层NetCDF读写接口xarray会调用它作为后端。matplotlib出图主力配合cartopy处理地图投影和海岸线。numpy/scipy做极坐标转笛卡尔坐标、插值、滤波的基础。cartopy地图投影库画飞行轨迹和扫描图时离不开。可选装pyart专门面向气象雷达的Python库做衰减订正、坐标变换时能省很多事。环境建议直接用conda创建避免源码编译Cartopy时踩坑。实测下来这套组合在Python 3.11下兼容性不错安装命令如下conda create -n nasa_impacts python3.11 conda activate nasa_impacts conda install -c conda-forge numpy xarray netcdf4 matplotlib cartopy scipy pyart2.2 读取NetCDF时的第一个坑先打印再动手拿到一个IMPACTS数据文件第一反应不要急着画图先用xarray打开并打印全局信息和变量结构。这一步能解决后面80%的困惑。import xarray as xr ds xr.open_dataset(IMPACTS_radar_20220213.nc) print(ds) # 常见变量 # reflectivity 反射率因子单位通常是 dBZ # doppler_velocity 多普勒径向速度 # spectrum_width 谱宽 # 常见坐标 # time, range, azimuth, elevation, latitude, longitude从实际经验看NASA这类公开档案里的文件变量命名还算规范但也存在不同架次之间维度顺序不一致的情况。比如有的文件把time放在第一维有的把range放在第一维如果直接按固定下标取值换一个架次就很容易报错或画出镜像错位图。用xarray的好处在于可以靠变量名取数比如ds[reflectivity].isel(time0)代码里不依赖硬编码下标程序的健壮性会好很多。模块内部我会建议做一个“字段标准化”函数把不同文件里可能的别名统一映射成标准字段名例如将REF、refl、reflectivity统一成reflectivity再把坐标统一排序为(time, range, azimuth)。这样后续绘图函数只需要关心标准字段名不用在每个函数里再判断一次变量存在性。3. 绘制第一张反射率图看懂极坐标扫描的可视化逻辑3.1 机载雷达扫描和地面雷达的差异气象雷达有两种基本扫描方式地面业务雷达常做平面位置扫描PPI也就是固定仰角、转一圈方位角得到一个圆锥面的回波。机载雷达则灵活得多除了PPI还可以做垂直方向的位置扫描RHI沿某个方位角上下扫得到垂直剖面。在设计可视化程序时我专门为这两种模式做了分离PPI模式画水平地图RHI模式画高度剖面。更复杂的情况是飞机本身在运动所以每一帧扫描的雷达中心位置都在变化画图时必须把经纬度坐标算对而不是假设雷达位于某个固定站点。3.2 从NetCDF到第一张水平扫描图下面这段代码实现的是读取某一时刻的反射率数据并转换到地图坐标系下进行绘制import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np def plot_ppi(ds, time_idx0, varreflectivity): refl ds[var].isel(timetime_idx) lat ds[latitude].isel(timetime_idx) lon ds[longitude].isel(timetime_idx) fig plt.figure(figsize(12, 8)) ax fig.add_subplot(1, 1, 1, projectionccrs.PlateCarree()) ax.set_extent([lon.min() - 0.5, lon.max() 0.5, lat.min() - 0.5, lat.max() 0.5], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth1.0) ax.add_feature(cfeature.STATES, linewidth0.6) im ax.pcolormesh(lon, lat, refl, cmapturbo, vmin-30, vmax40, transformccrs.PlateCarree()) plt.colorbar(im, axax, labelReflectivity (dBZ), fraction0.025, pad0.04) ax.set_title(fIMPACTS Reflectivity PPI - Time {time_idx}, fontweightbold) return fig, ax这里有个关键点pcolormesh需要二维坐标数组通常NetCDF里的经纬度在机载雷达场景下就是二维的直接用即可。如果遇到一维坐标要先做np.meshgrid处理。颜色映射我用的是turbo视觉效果比jet更平滑比viridis更贴近业务雷达图习惯比較适合在公众号、PPT和论文初稿里快速展示。但如果你要把图投到某个要求严格的期刊还是建议查一下期刊对色标的具体要求有的期刊明确要求用感知均匀的色标。3.3 极坐标转笛卡尔RHI剖面图的实现垂直剖面图能直观展示雪带从低层到高层的演变。机载雷达的RHI扫描通常以飞机为中心记录不同距离和仰角上的回波。绘制时先把极坐标(距离, 仰角)转换成笛卡尔坐标(水平距离, 高度)再用pcolormesh画出。def plot_rhi(ds, az_idx0, varreflectivity, max_range_km60): rng ds[range] / 1000.0 # 转为公里 elev np.deg2rad(ds[elevation].isel(azimuthaz_idx)) az np.deg2rad(ds[azimuth].isel(azimuthaz_idx)) R, E np.meshgrid(rng, elev) x R * np.cos(E) * np.cos(az) z R * np.sin(E) ref ds[var].isel(azimuthaz_idx) fig, ax plt.subplots(figsize(14, 6)) ax.pcolormesh(x, z, ref, cmapviridis, vmin-30, vmax40) ax.invert_yaxis() ax.set_xlabel(Horizontal Distance (km)) ax.set_ylabel(Altitude (km)) ax.set_title(IMPACTS RHI vertical cross-section) return fig, ax有一点必须提醒机载雷达的range和地面雷达一样是从雷达天线出发的斜距不是水平距离。如果你的数据显示了地形回波或非常明显的底部杂波先别急着怀疑仪器先检查一下是否忘记做高度订正。飞机本身就有高度因此以海平面或飞机所在高度作为参考零点的差异会直接影响剖面图的解释。建议在程序里加一个参数reference_altitude_km默认取飞机GPS高度实际画图时再根据需求平移。4. 程序化出图飞行轨迹、剖面和系列帧的高效组织4.1 自动化输出目录结构真正实用的可视化程序不是只有一个画图函数而是能把一批文件、一批时间点自动处理成整齐的图集。我建议项目里保留如下目录结构nasa_impacts/ ├── config.yaml ├── nasa_impacts/ │ ├── load.py │ ├── plot_ppi.py │ ├── plot_rhi.py │ ├── flight_track.py │ └── make_quicklook.py ├── data/ │ └── 20220213/ ├── outputs/ │ ├── ppi/ │ ├── rhi/ │ ├── track/ │ └── overview/这样做的目的很直接数据源不乱动输出按图和类型分离等哪天需要做对比图或者做成Web服务时直接读对应目录就行不用再满硬盘翻文件。在config.yaml里可以指定数据文件列表、要绘制的变量名、色标范围、输出格式、DPI等参数。整个程序只需要一个入口命令python -m nasa_impacts.make_quicklook --config config.yaml --flight 202202134.2 飞行轨迹叠加时间着色除了雷达扫描图飞行轨迹本身也是重要可视化对象。把飞机轨迹按时间顺序连成线并按时间或高度着色能快速看出飞机是在风暴哪个高度层穿插。这里用LineCollection会比较方便from matplotlib.collections import LineCollection def plot_track(ds_track, t, lat, lon, alt, cmap_nameplasma): points np.array([lon, lat]).T.reshape(-1, 1, 2) segments np.concatenate([points[:-1], points[1:]], axis1) lc LineCollection(segments, cmapcmap_name, linewidth2) lc.set_array(alt[:-1]) fig, ax plt.subplots(figsize(10, 10), subplot_kw{projection: ccrs.PlateCarree()}) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.STATES) ax.add_collection(lc) ax.autoscale() plt.colorbar(lc, axax, labelAltitude (km)) return fig, ax飞行轨迹图是检查数据质量问题的一把好手。如果轨迹出现不连续的跳变通常说明GPS定位信号有丢失或者数据插值出了问题这种异常在统计指标里很难发现但画在图上非常扎眼。4.3 循环生成“帧序列”为了观察风暴演变我会按固定时间间隔生成多帧PPI或RHI图再组合成动画或GIF。程序里写一个循环import imageio.v2 as imageio from pathlib import Path writer imageio.get_writer(impacts_loop.gif, fps3) for t in range(0, len(ds[time]), 5): fig, ax plot_ppi(ds, time_idxt) fname foutputs/ppi/frame_{t:04d}.png fig.savefig(fname, dpi120, bbox_inchestight) plt.close(fig) writer.append_data(imageio.imread(fname)) writer.close()这里建议保存PNG帧之后再用imageio合成GIF不要边画边往内存里塞大数组否则几十分钟的飞行数据足够把16GB内存吃满。主力分析时还是保存单帧PNG更实用GIF只用于汇报和快速浏览。5. 实际工程中的坑给nasa_impacts项目的七条避坑记录这几条都是我在实际做可视化程序时踩过的有些是气象领域特有的有些是任何NetCDF处理项目都会遇到。写在这儿希望你拿到新数据时不会被同样的石头绊住。5.1 维度顺序不一致不同架次数据里变量的维度顺序不完全一致有的把time放最前面有的把range放中间。肉眼看不出来但脚本里的固定索引会直接取错。用xarray的isel和sel可以避免但前提是先打印ds确认维度名。更保险的做法是写一个normalize_dims(ds)函数把标准的坐标名重新排序后返回保证生产线稳定。5.2 单位不是看上去的单位反射率数据的单位通常是dBZ但有些文件里存的是线性反射率Z数值范围大不相同。如果直接用默认色标轻则图像饱和度过高重则完全看不出回波结构。读取数据后务必检查units属性如果值是mm^6/m^3这类线性单位先做10 * np.log10(z)转成dBZ再绘图。5.3 缺测值没有掩膜NetCDF文件里缺测点通常用NaN或一个大数表示比如-9999。pcolormesh遇到NaN会透明不画但遇到-9999会画出极值颜色。处理方式是在绘制前统一处理refl refl.where(refl -100, np.nan)与此同时也要检查经纬度数组里是否混入缺测值那会导致投影点跑到非洲去。5.4 色标选择比想象中重要我早期喜欢用jet后来发现等降水量线、论文审稿人、视觉正常的人都不太买账。现在默认用两种大范围回波用turbo学术剖面图用viridis或pyart内置的Reflectivity colormap。如果目标读者是气象业务人员NWS Reflectivity色标可能更符合他们的看图习惯这个可以把matplotlib.colors.ListedColormap配一组颜色分界点来实现。5.5 Cartopy与现有坐标不匹配如果直接在ax.pcolormesh里传入从NetCDF读出的经纬度且没有指定transformccrs.PlateCarree()图上所有数据点会落在投影坐标系的原点附近画出来就是一团乱麻。记住一条原则数据坐标用transform参数声明地图框用projection参数声明两个坐标系不同时Cartopy会自动转换但这层转换必须在绘图调用里显式写出来。5.6 大数据文件内存占用失控一整架次飞行雷达原始数据可能包含上千个时间步每个时间步又有多条扫描线。可以读入内存但不要一下把整个数据集的所有变量都载入。我通常用xr.open_dataset(..., chunks{time: 20})启用Dask分块绘图前只取所需变量和时间切片。处理完一批后立即del fig, ax和plt.close(fig)否则Matplotlib在循环里积累的Figure对象会越占越多。5.7 文件按架次分组时的时间标准化IMPACTS数据文件的命名通常带日期和架次号但不同年份的数据归档方式可能不同。项目里建议在配置文件里维护一个字典把架次名映射到文件路径再用pandas.to_datetime解析时间坐标。这样后面做时间轴统一、跨架次拼接时不会因为UTC偏移和本地时区的问题而对不上号。6. 从命令行脚本到可视化工具包的演进建议一旦把单张图画通下一步就是让整个流程更像一个“程序”而不是一坨在Jupyter Notebook里反复粘贴的代码。我在项目里做过的三个改进非常值得推荐。第一是加QuickLook模式。所谓QuickLook就是每处理一个数据文件时自动生成一张概览图左半边是飞行轨迹右半边是不同高度层的回波拼图底部再附一行时间标签和文件名。这种图不是给最终报告用的而是给数据检查用的好处是几乎不用看日志就能判断哪些架次的数据质量有问题。第二是增加配置文件。所有色标范围、变量名称、飞行日期、区域范围都放在config.yaml里避免每次改场景都需要改动代码。比如某个数据文件把经纬度变量的名字从latitude改成了lat只需要修改配置文件里的lat_var: lat不需要动任何绘图函数。第三是预留交互式前端接口。现在气象数据可视化经常被问到“能不能把数据发布到大屏上”其实完全可以先把PNG图集成到Grafana或自建Web页面里再挂上一张按时间排列的缩略图墙。实现方式也不复杂让程序每次跑完输出一个index.html把生成的图按文件列表写进一个简单网格页面浏览器打开就能看到整个架次的图集。以我个人的经验做这类NASA科研数据的可视化程序最难的不是写代码而是建立对数据的“手感”。你需要在一次次出图、一遍遍回头调整坐标轴和色标的过程中把每条回波带、每个数据缺口都变成自己能识别的信号。等你能在0.5倍速看飞行动画时直接说出“这架次北段在17公里高度有一个干层”就说明这个可视化程序已经真正长在你脑子里了。到那个阶段画图已经不是任务而是你理解风暴的一种方式了。本文还有配套的精品资源点击获取