简介:这份资源围绕NASA IMPACTS暴风雪研究计划展开,面向气象学、遥感与数据科学方向的学习者和研究者,帮助其直观理解美国东海岸冬季风暴的微物理探测与多平台观测数据。IMPACTS是东海岸暴风雪三十年来首次综合研究,借助P-3B与ER-2飞机、地面雷达及GPM、GOES-16、JPSS等多颗卫星协同采样,资源以Python脚本对相关航次与气象导航数据做可视化呈现。压缩包共28个文件,约21.47MB,包含21个csv观测数据、3张png航迹与对比图、2个py可视化脚本,以及license和README说明文档,结构清晰,便于按航次与平台检索。已有191人学习下载。读者可据此复现P-3B与ER-2联合航班的探测轨迹,理解前锋带浅雪带等天气过程的观测组织方式,并借助脚本快速上手NASA公开气象数据的读取、绘图与分析流程。
1. 从一份气象数据包说起:nasa_impacts 到底能拿来干什么
如果你手头正躺着一批 NASA IMPACTS 暴风雪研究计划的观测数据,却不知道怎么把它们从原始文件变成能看的图、能跑的统计,那这个nasa_impacts可视化项目就是冲着你来的。IMPACTS 全称 Investigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms,是美国航天机构牵头的一个冬季风暴野外观测计划,用 P-3 飞机搭载雷达、微波辐射计、云粒子探头等设备,在暴风雪里穿云采样,攒下了大量一手数据。这些数据对做降水微物理、雷达反演、数值模式验证的人来说是实打实的富矿,但原始格式分散、变量命名不统一、时间轴对不齐,直接上手很容易卡在第一步。nasa_impacts这个项目做的事,就是把这批数据用 Python 串起来,做成一套可复现的可视化流程,让你不用从零写解析脚本,就能把 P-3 飞行轨迹、雷达剖面、降水粒子分布画出来。它适合大气科学方向的研究生、做冬季降水分析的从业者,以及想拿真实观测数据练手的数据科学人。下面我按自己拆包复现的顺序,把这份资源从结构到跑通再到排错讲一遍。
2. 拆开数据包:目录结构、依赖与 P-3 数据组织逻辑
2.1 先看清仓库里有什么,别急着 pip install
拿到一个科研数据可视化项目,我习惯先不跑代码,而是把目录树和依赖文件翻一遍。nasa_impacts这类项目的典型结构是:一个data/或datasets/目录放样本数据或下载脚本,一个notebooks/放 Jupyter 分析流程,一个src/或包目录放可复用的解析和绘图函数,根目录下是requirements.txt、environment.yml和README。先看README里有没有写数据获取方式,因为 IMPACTS 的原始数据体积很大,仓库通常不会直接塞进来,而是给你一个下载入口或子集样本。
我一般会先执行这几条命令摸清底细:
# 看目录层级,重点关注 data、notebooks、src 三个位置 find . -maxdepth 2 -type d | sort # 看依赖清单,判断是 conda 还是 pip 为主 cat requirements.txt 2>/dev/null || cat environment.yml 2>/dev/null # 看有没有数据下载脚本或说明 ls -la data/ 2>/dev/null; grep -ri "download" README* 2>/dev/null | head逻辑说明:find限制两层深度是为了不被深层缓存目录刷屏;requirements.txt和environment.yml二选一读取,是因为不同项目习惯不同,conda 环境更适合带地理空间库(如 cartopy、geopandas)的项目,pip 则轻量。参数上-maxdepth 2可根据实际调整,如果项目嵌套深就加到 3。这一步的产出是:你知道了数据在哪、依赖多重、入口是 notebook 还是脚本。
2.2 P-3 飞机数据的组织方式与关键变量
IMPACTS 的核心观测平台是 P-3 飞机,它上面挂载的仪器产生的数据大致分几类:导航与状态数据(经纬度、高度、时间戳)、遥感数据(机载雷达反射率、微波辐射计亮温)、原位微物理数据(粒子图像、粒径谱、液态水含量)。这些数据在原始分发时往往是按飞行架次(flight)分文件,每个文件里变量名带仪器前缀,时间戳可能是秒级偏移或 UTC 字符串。
nasa_impacts项目要解决的第一件事就是把这些异构变量对齐到统一的时间轴和坐标系。常见做法是用xarray把 NetCDF 文件读成 Dataset,用pandas处理时间索引,再用cartopy做地图投影。下面是一段我常用的读取与初步检查代码,你可以直接套:
import xarray as xr import pandas as pd # 打开一个架次的 NetCDF 文件,decode_times 让时间自动解析 ds = xr.open_dataset("data/IMPACTS_P3_flight01.nc", decode_times=True) # 打印变量清单和维度,先搞清楚有什么 print(ds.data_vars) print(ds.dims) # 把时间坐标转成 pandas 索引,方便后续重采样 time_index = pd.to_datetime(ds["time"].values) print("时间范围:", time_index.min(), "到", time_index.max()) # 检查经纬度是否在合理范围,暴风雪区域一般落在北美东海岸 print("纬度范围:", float(ds["lat"].min()), float(ds["lat"].max())) print("经度范围:", float(ds["lon"].min()), float(ds["lon"].max()))逻辑说明:decode_times=True是关键参数,很多 NetCDF 的时间是 "seconds since 1970-01-01" 这类编码,不解析的话后面画图时间轴会乱。data_vars和dims先看结构再动手,避免对着不存在的变量名写代码。经纬度范围检查是血泪经验——曾经有次数据里混进了地面测试架次,经纬度是实验室坐标,直接画图会得到一张莫名其妙的图。参数上,如果你的文件时间变量不叫time,用ds.attrs或ncdump -h先确认实际名字。
2.3 环境搭建:地理空间库的版本坑
这类项目最容易翻车的不是业务代码,而是环境。cartopy、geopandas、shapely这几个库对底层 GEOS、PROJ 版本敏感,pip 直接装经常报编译错误。我一般优先用 conda 建环境:
# 用 conda 建一个独立环境,指定 python 版本 conda create -n impacts python=3.10 -y conda activate impacts # 先装地理空间全家桶,走 conda-forge 渠道 conda install -c conda-forge cartopy geopandas xarray netcdf4 matplotlib -y # 再装项目自身的轻量依赖 pip install -r requirements.txt逻辑说明:conda-forge渠道的二进制包已经处理好 GEOS/PROJ 依赖,比 pip 源码编译省心得多。Python 版本选 3.10 是因为多数科研库对 3.11+ 的支持还在跟进,3.10 是当前稳妥区间。先 conda 装重库、再 pip 装项目依赖的顺序,能避免 pip 把 conda 装好的库覆盖成不兼容版本。如果你只有 pip,那就先pip install shapely --no-binary之类折腾,但我不推荐,时间成本太高。
3. 跑通可视化流程:从飞行轨迹到雷达剖面
3.1 飞行轨迹与时间-高度剖面图
把数据读进来之后,第一张该出的图是 P-3 的飞行轨迹叠加在地图上,再配一张时间-高度剖面。轨迹图告诉你飞机在哪飞,剖面图告诉你它在哪个高度采样,两张图对照着看,才能判断某段雷达数据对应的空间位置。
import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig = plt.figure(figsize=(12, 5)) # 左图:飞行轨迹 ax1 = fig.add_subplot(121, projection=ccrs.PlateCarree()) ax1.add_feature(cfeature.COASTLINE) ax1.add_feature(cfeature.STATES, linewidth=0.3) ax1.plot(ds["lon"], ds["lat"], color="red", linewidth=0.8, transform=ccrs.PlateCarree()) ax1.set_title("P-3 Flight Track") ax1.set_extent([-85, -60, 30, 48], crs=ccrs.PlateCarree()) # 右图:时间-高度剖面 ax2 = fig.add_subplot(122) ax2.plot(time_index, ds["altitude"], color="blue", linewidth=0.8) ax2.set_xlabel("Time (UTC)") ax2.set_ylabel("Altitude (m)") ax2.set_title("Altitude Profile") plt.xticks(rotation=30) plt.tight_layout() plt.savefig("flight_overview.png", dpi=150)逻辑说明:PlateCarree是最常用的等经纬投影,set_extent把范围锁在北美东海岸,避免全球地图上轨迹缩成一点。transform=ccrs.PlateCarree()这个参数不能漏,否则经纬度会被当成投影坐标画错位置。右图直接用 pandas 时间索引当 x 轴,matplotlib 会自动处理刻度。dpi=150是出图清晰度和文件大小的折中,投稿级可以提到 300。参数上,set_extent的四个值按你的实际架次区域改,别照抄。
3.2 雷达反射率剖面的绘制与色标选择
IMPACTS 的机载雷达数据是重头戏,反射率(reflectivity)随高度和距离的变化能反映降雪带结构。画这类剖面图,关键是色标(colormap)和坐标轴方向——雷达剖面通常是距离-高度,距离沿飞行路径累积。
import numpy as np # 假设 ds 里有 reflectivity 变量,维度是 (time, range) # 先构造距离轴,用飞行速度积分或直接用雷达距离门 distance = ds["range"].values # 单位通常是米或公里 altitude = ds["altitude"].values # 用 pcolormesh 画剖面,注意转置让距离在 x 轴 fig, ax = plt.subplots(figsize=(10, 5)) mesh = ax.pcolormesh(distance, altitude, ds["reflectivity"].T, cmap="viridis", vmin=-20, vmax=40, shading="auto") ax.set_xlabel("Distance (km)") ax.set_ylabel("Altitude (m)") ax.set_title("Radar Reflectivity Profile") cbar = plt.colorbar(mesh, ax=ax) cbar.set_label("Reflectivity (dBZ)") plt.tight_layout() plt.savefig("radar_profile.png", dpi=150)逻辑说明:pcolormesh比imshow更适合非均匀网格,shading="auto"让 matplotlib 自动判断网格对齐方式。vmin=-20, vmax=40是雷达反射率的常用范围,低于 -20 基本是晴空噪声,高于 40 是强降水或冰雹,暴风雪研究里重点看 0 到 25 dBZ 这段。.T转置是因为原始数据维度顺序可能是 (range, time),而画图需要 x 轴是距离。色标选viridis是因为它对色盲友好且打印成灰度仍有层次,如果你要突出弱回波,可以换jet但我不推荐,jet的视觉误导是老问题了。
3.3 微物理变量的联合展示
光看雷达还不够,IMPACTS 的价值在于它同时有原位微物理观测。把液态水含量(LWC)、冰水含量(IWC)、粒子浓度和雷达反射率放一张多面板图里,能看出微物理量和遥感量的对应关系。
fig, axes = plt.subplots(4, 1, figsize=(10, 10), sharex=True) variables = ["lwc", "iwc", "particle_conc", "reflectivity"] labels = ["LWC (g/m3)", "IWC (g/m3)", "Particle Conc (#/L)", "Reflectivity (dBZ)"] colors = ["green", "blue", "orange", "purple"] for ax, var, label, color in zip(axes, variables, labels, colors): if var in ds: ax.plot(time_index, ds[var], color=color, linewidth=0.7) ax.set_ylabel(label) else: ax.text(0.5, 0.5, f"{var} not found", transform=ax.transAxes, ha="center") axes[-1].set_xlabel("Time (UTC)") plt.xticks(rotation=30) plt.tight_layout() plt.savefig("microphysics_panel.png", dpi=150)逻辑说明:sharex=True让四个面板共用时间轴,方便上下对照。循环里加了if var in ds判断,是因为不同架次的仪器工作状态不同,某些变量可能缺失,硬写会直接报 KeyError。linewidth=0.7是时间序列的常用细线,太粗会糊成一片。这段代码的产出是一张能直接放进报告的多面板图,参数上你可以按需增删变量,但建议不超过 5 个面板,否则单张图信息过载。
4. 避坑与排查:那些让我重跑一整天的坑
4.1 时间戳时区不一致导致轨迹错位
现象:飞行轨迹图和雷达剖面对不上,飞机位置和回波区域偏移了几十公里。原因:导航数据的时间戳是 UTC,而某些仪器数据的时间戳是本地时或 GPS 秒,直接合并会错位。解决:统一转成 UTC 再对齐,用pd.to_datetime(..., utc=True)强制带时区,合并前先ds = ds.assign_coords(time=("time", unified_time))重写时间坐标。我现在的习惯是读进来第一件事就打印时间范围,和飞行日志核对。
4.2 NetCDF 变量名大小写与下划线不一致
现象:代码里写ds["Reflectivity"]报 KeyError,但ncdump看明明有。原因:不同仪器团队命名习惯不同,有的是reflectivity,有的是Reflectivity,还有的是radar_reflectivity。解决:别硬编码,先print(list(ds.data_vars))看实际名字,或者写一个模糊匹配函数[v for v in ds.data_vars if "reflect" in v.lower()]。这个坑我踩过不止一次,后来养成了先列变量再写代码的习惯。
4.3 cartopy 地图要素下载失败
现象:第一次运行绘图代码卡在Downloading: 100%然后超时,或者报URLError。原因:cartopy 的cfeature.COASTLINE等要素首次使用会从远程下载自然地球数据,网络不稳就失败。解决:提前手动下载或配置本地数据目录,设置环境变量CARTOPY_DATA_DIR指向已缓存目录。如果实在下不动,用cfeature.COASTLINE.with_scale("110m")降低精度,文件更小。这个坑的教训是:别在演示前才第一次跑绘图代码。
4.4 内存溢出:一次读入所有架次
现象:读第三个架次时进程被 kill,或者机器卡死。原因:IMPACTS 单架次数据可能几个 GB,一次性xr.open_mfdataset全部架次直接吃满内存。解决:用xr.open_dataset逐架次处理,或者用 dask 分块chunks={"time": 1000}惰性加载。我一般会先估算:单文件大小 × 架次数,超过内存一半就分批。参数上chunks的块大小按你的内存调,8GB 内存建议 time 维块不超过 2000。
4.5 色标范围不当导致弱回波看不见
现象:雷达剖面图上一片深色,看不出降雪带结构。原因:vmin/vmax设得太宽,弱回波被压缩到色标底部。解决:先画直方图看反射率分布,把vmin设在 5% 分位数、vmax设在 95% 分位数,或者针对暴风雪场景手动设vmin=-10, vmax=25。这个坑的本质是:默认色标范围是给通用场景的,科研图必须按数据实际分布调。
5. 进阶技巧:把可视化流程封装成可复用函数
跑通单架次之后,你肯定不想每换一个架次就复制粘贴改路径。我一般会把读取、对齐、绘图三步封装成函数,用一个配置字典控制参数。下面这个模式我用了很多次,改改就能套:
def load_flight(nc_path, time_var="time", lat_var="lat", lon_var="lon"): """读取单个架次并统一时间坐标""" ds = xr.open_dataset(nc_path, decode_times=True) ds = ds.assign_coords(time=pd.to_datetime(ds[time_var].values, utc=True)) ds = ds.rename({lat_var: "lat", lon_var: "lon"}) return ds def plot_overview(ds, extent=(-85, -60, 30, 48), out="overview.png"): """出轨迹+高度剖面图""" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), subplot_kw={"projection": ccrs.PlateCarree()}) # ... 绘图逻辑同上,略 plt.savefig(out, dpi=150) plt.close() # 批量处理 flights = ["flight01.nc", "flight02.nc", "flight03.nc"] for f in flights: ds = load_flight(f"data/{f}") plot_overview(ds, out=f"figs/{f.replace('.nc', '.png')}")逻辑说明:load_flight把时间对齐和变量重命名收进去,plot_overview只管画图,职责分离。plt.close()不能漏,批量循环里不关 figure 会累积内存直到崩。参数上extent给了默认值但可覆盖,out路径按架次名生成,避免覆盖。这套封装跑 20 个架次也就几分钟,比手动改代码快得多。
验证封装是否可靠,我有个笨办法但很管用:拿两个已知差异明显的架次(比如一个飞得高、一个飞得低)跑一遍,对比高度剖面图的 y 轴范围是否合理,再抽查一张轨迹图的经纬度是否落在预期区域。如果两架次图看起来一模一样,那多半是数据没读对或者坐标没对齐。
从那以后我每次拿到新的观测数据包,都强制先跑一遍「变量清单 + 时间范围 + 经纬度范围」三件套检查,确认数据本身没问题再写业务代码。这个习惯帮我省下了至少两次重跑一整天的后悔药。希望这份拆解能帮你把nasa_impacts这批暴风雪数据真正用起来,少走点我走过的弯路。
本文还有配套的精品资源,点击获取