简介本资源面向气象学、大气科学及相关专业本科生与科研初学者提供一套基于Python的典型动力气象诊断量计算实战案例重点解决涡度、散度、涡度平流和温度平流等关键物理量的编程实现问题。压缩包共658个文件包含273个adfGRIB格式气象数据、117个nit与117个dat数值模式输出及中间结果、42个log运行日志、39个001分卷文件大体积数据切片以及少量ovr地理配准文件、xml元数据和doc/pdf文档总大小29.6MB结构完整覆盖数据读取、坐标处理、差分计算、可视化全流程。已有3364人学习下载资源内含可直接运行的Python核心代码、配套说明文档及多组实测数据集所有脚本均针对真实气象场设计注释详尽支持NCL/GrADS用户快速迁移至Python生态并附有常见报错解析与单位转换对照表显著降低气象诊断编程入门门槛。1. 项目缘起从气象数据到物理量诊断如果你处理过气象再分析数据比如ERA5或者GFS的格点资料大概率会和我有一样的经历下载下来的数据包变量名是u、v、t、z分别代表纬向风、经向风、温度和位势高度。这些是基础场但气象分析和预报中真正用来判断天气系统强度、识别锋区、预报降水落区的往往是它们的衍生物理量比如涡度、散度、平流项。我第一次需要计算涡度平流来做强对流天气分析时翻遍了常用的气象软件和库的文档。像Grads、NCL这类传统工具当然能算但脚本写起来不够灵活数据前后处理也麻烦。用Matlab吧矩阵运算方便但涉及到地球球坐标下的差分各种cos(lat)的因子一掺和代码就容易写乱而且对于批量处理大量时序数据效率也是个问题。直到我把目光彻底转向Python配合xarray、metpy和numpy才发现这条路径既清晰又高效。这次我就把自己在业务和研究中反复打磨的一套计算方案分享出来核心就是如何用Python准确、高效地从标准格点风场、温度场中计算出涡度、散度、涡度平流和温度平流。这几个量是天气动力学诊断的基石。简单来说涡度描述空气块的旋转程度正涡度对应气旋式旋转是低压、风暴发展的标志散度描述空气块的辐散辐合低层辐合高层辐散是垂直上升运动的动力条件。而平流描述物理量被风输送的过程涡度平流是预报槽脊移动的关键温度平流冷暖平流直接关联锋生、锋消和垂直运动。算对了这些一张天气图在你眼里就不再是静止的等高线和等温线而是流动的、相互作用的物理过程。本文将完全基于Python生态假设你已经有了一份netCDF或GRIB格式的格点数据例如从ERA5下载的。我们会从数据读取开始一步步推导公式在离散格点上的实现方法处理球坐标下的特殊问题比如经纬度网格并给出完整的、可复现的代码。过程中我会重点分享几个我踩过的坑比如为什么直接用numpy.gradient算出的涡度在极区和高纬度会“爆炸”如何处理地图投影和格距变化以及在计算平流时是选择“中心差分”还是“一次逆风差分”不同的选择对结果有什么影响这些细节才是从“能算”到“算得准、用得对”的关键。2. 环境准备与数据读取构建可复现的计算基础工欲善其事必先利其器。一个稳定、一致的计算环境是后续所有工作的前提。我强烈建议使用conda来管理你的Python环境这能完美解决气象包依赖复杂的问题。首先创建一个专属环境conda create -n meteorology python3.9 conda activate meteorology接下来安装核心依赖。这里有个关键点metpy是计算气象物理量的神器但它的一些高级依赖特别是cartopy在pip安装时容易出错。最稳妥的方式是通过conda的conda-forge通道来安装它能自动处理好所有二进制依赖。conda install -c conda-forge xarray netcdf4 cfgrib metpy cartopy numpy scipy matplotlib jupyter这条命令一次性安装了我们需要的大部分包xarray用于优雅地处理网格数据netcdf4和cfgrib是后端引擎用于读取不同格式的数据metpy提供气象计算函数cartopy用于绘图numpy和scipy是数值计算基石matplotlib和jupyter则是可视化和交互环境。数据方面我们以欧洲中期天气预报中心ECMWF的ERA5再分析数据为例这是目前最易获取且质量很高的全球数据。假设你已经通过CDS API下载好了单层或多层的数据通常是一个netCDF文件里面包含了u1010米纬向风、v1010米经向风、t温度、z位势高度等变量其维度通常是(time, latitude, longitude)。使用xarray打开数据非常直观import xarray as xr import numpy as np import metpy.calc as mpcalc from metpy.units import units # 打开数据集 ds xr.open_dataset(your_era5_data.nc) # 查看数据结构和变量 print(ds)这里会遇到第一个实际操作细节单位。ERA5数据的风场单位通常是m s**-1温度是K位势高度是m**2 s**-2这个需要除以9.80665才是位势米。metpy的强项之一就是带单位的计算。我们需要将xarray.DataArray转换为metpy能识别的带单位数组。但注意直接使用mpcalc的函数时它通常能自动从xarray属性中解析单位。为了保险起见可以这样处理# 选取一个时间点和一个层次如果是多层数据 # 例如我们计算850hPa上的物理量 time_idx 0 # 第一个时次 level 850 * units.hPa # 指定层次 # 提取变量并确保它们具有正确的单位属性 u ds[u].sel(levellevel, timeds.time[time_idx]).metpy.quantify() v ds[v].sel(levellevel, timeds.time[time_idx]).metpy.quantify() temp ds[t].sel(levellevel, timeds.time[time_idx]).metpy.quantify() # 对于高度场ERA5中通常是‘z’单位是 m^2/s^2 geopotential ds[z].sel(levellevel, timeds.time[time_idx]) # 将位势转换为位势高度位势米 height geopotential / 9.80665 height.attrs[units] m height height.metpy.quantify() # 获取经纬度坐标并附加单位 lats u.latitude.metpy.quantify() lons u.longitude.metpy.quantify()这段代码里的.metpy.quantify()是关键一步它尝试将数据转换为pint库管理的、带物理单位的数组。如果原始数据有units属性ERA5通常有这一步会自动完成。如果没有你需要手动赋值如u.attrs[units] m/s。注意数据插值与格点类型有时你下载的数据可能是高斯格点或其他投影坐标而我们的计算默认在规则的经纬度网格上进行。metpy的许多函数要求输入latitude和longitude坐标。如果你的数据坐标名不是这个或者需要插值到规则网格可以使用xarray的rename和interp方法。例如ERA5的纬度坐标有时叫latitude有时叫lat统一一下会更方便ds ds.rename({lat: latitude, lon: longitude})。3. 核心原理与公式离散化理解每个微分背后的物理在连续的大气运动中这些物理量有精确的数学定义。但我们的数据是离散的格点所以必须理解如何将连续的微分方程转化为离散的差分格式。这是保证计算精度的理论基础。3.1 相对涡度与散度的定义相对涡度ζ和散度D在球坐标下的表达式为ζ (1 / (a cosφ)) * (∂v/∂λ - ∂(u cosφ)/∂φ) D (1 / (a cosφ)) * (∂u/∂λ ∂(v cosφ)/∂φ)其中a是地球半径约6371kmφ是纬度λ是经度u是纬向风东正西负v是经向风北正南负。看到公式里的cosφ纬度余弦和分母中的a了吗这就是在球面上计算必须考虑的曲率效应和格距随纬度变化。在低纬度经度方向上的实际距离a cosφ Δλ很大而在高纬度同样的经度差对应的实际距离很小。如果直接用简单的Δv/Δx - Δu/Δy笛卡尔坐标近似在低纬度误差尚可接受但在中高纬度尤其是靠近极地计算结果会完全失真数值可能异常巨大这就是我前面提到的“爆炸”问题。3.2 平流项的计算平流是标量属性如涡度ζ、温度T被风场输送的过程。其表达式为涡度平流 Adv_ζ - (u * ∂ζ/∂x v * ∂ζ/∂y) 温度平流 Adv_T - (u * ∂T/∂x v * ∂T/∂y)注意前面的负号它意味着物理量沿着风的方向减少。例如北风v为负吹向温度梯度为正向北温度增加的区域将导致-v * ∂T/∂y为正即暖平流。这里有一个非常重要的算法选择如何计算导数∂ζ/∂x和∂ζ/∂y中心差分最常用精度二阶。例如∂ζ/∂x ≈ (ζ[i1] - ζ[i-1]) / (2Δx)。但它假设流场平滑在物理量梯度极大的区域如急流轴、锋区可能产生数值振荡甚至出现“伪极值”。一次逆风差分具有数值耗散性格式稳定。它总是用上游点的信息来计算下游点的变化。虽然精度为一阶但在实际天气图分析中有时反而能产生更平滑、更符合预报员直观的平流场因为它抑制了小尺度噪音。在业务中对平滑后的场进行计算时中心差分是主流但对原始分辨率数据需要谨慎。在我们的实现中我们将使用metpy的函数它默认采用中心差分并自动处理球坐标下的格距变化。这是最通用和推荐的做法。3.3 离散化实现的关键metpy.calc的幕后工作当我们调用mpcalc.vorticity(u, v)时metpy在内部做了什么它首先从输入的DataArray中提取经纬度坐标并转换为弧度。计算每个格点的Δx和Δy。Δx a * cosφ * ΔλΔy a * Δφ。这里Δλ和Δφ是经纬度的弧度差。使用中心差分公式计算风场的空间导数∂u/∂x,∂v/∂x,∂u/∂y,∂v/∂y。注意在计算∂/∂y时对于含有cosφ的项如u cosφ它会在差分前先乘以cosφ差分后再除以cosφ这正是球坐标公式的离散体现。最后组合出涡度ζ ∂v/∂x - ∂u/∂y和散度D ∂u/∂x ∂v/∂y。理解这个过程你就能明白为什么用numpy.gradient(u, lons, lats, axis(2,1))直接算会出问题——它默认是在笛卡尔坐标下进行等间距差分没有考虑球面几何。而metpy帮我们封装了这一切复杂性。4. 实战计算一步步生成诊断场理论清晰后我们开始动手计算。整个过程会非常简洁这得益于metpy的封装。4.1 计算相对涡度和散度# 计算相对涡度 (Relative Vorticity) vorticity mpcalc.vorticity(u, v) print(f涡度范围: {vorticity.min().values:.2e} 到 {vorticity.max().values:.2e} /s) # 通常转换为更常用的量级 10^-5 /s vorticity_per_5e5 vorticity.to(10^-5/s) # 计算散度 (Divergence) divergence mpcalc.divergence(u, v) print(f散度范围: {divergence.min().values:.2e} 到 {divergence.max().values:.2e} /s) divergence_per_5e5 divergence.to(10^-5/s)两行核心代码涡度和散度就计算完毕了。mpcalc.vorticity和mpcalc.divergence函数会自动识别数据的经纬度坐标和单位并应用正确的球面差分公式。输出结果的单位是/s气象上常用10^-5 /s来度量所以我们可以方便地转换一下。4.2 计算涡度平流和温度平流接下来计算平流。这里需要特别注意函数的参数顺序和单位。# 计算涡度平流 (Vorticity Advection) # 参数顺序风速u分量风速v分量涡度场经纬度 vorticity_advection mpcalc.advection(vorticity, [u, v]) print(f涡度平流范围: {vorticity_advection.min().values:.2e} 到 {vorticity_advection.max().values:.2e} /s^2) # 转换为常用单位 10^-9 /s^2 vort_adv_per_9e9 vorticity_advection.to(10^-9/s^2) # 计算温度平流 (Temperature Advection) # 注意温度平流是温度场被风场平流。metpy.advection函数第一个参数是被平流的标量场。 temperature_advection mpcalc.advection(temp, [u, v]) print(f温度平流范围: {temperature_advection.min().values:.2e} 到 {temperature_advection.max().values:.2e} K/s) # 转换为常用单位 10^-5 K/s temp_adv_per_5e5 temperature_advection.to(10^-5 K/s)mpcalc.advection函数同样封装了球坐标下的平流计算。它内部先计算标量场涡度或温度在经向和纬向上的梯度考虑球面再与风场点乘并加上负号。重要提示边界上的NaN值由于中心差分需要用到i1和i-1的格点因此计算出的涡度、散度和平流场在区域的最北、最南、最东、最西一圈格点上的值会是NaN无效值。这是差分方法的固有特性不是错误。在绘图或后续分析时你需要决定是保留这些NaN绘图时显示为空白还是进行填充例如用外推或简单复制邻近值。在大多数天气分析中我们关注的是区域内部所以通常可以直接用xarray的.where()方法或matplotlib绘图时忽略NaN。4.3 结果验证与快速可视化计算完成后快速画个图验证一下结果是好习惯。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 创建一个地图投影例如兰伯特投影 proj ccrs.LambertConformal(central_longitudelons.mean().values, central_latitudelats.mean().values) # 绘制850hPa相对涡度 fig, ax plt.subplots(figsize(12, 8), subplot_kw{projection: proj}) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5) # 绘制填色图 contourf ax.contourf(lons, lats, vorticity_per_5e5, levelsnp.linspace(-20, 20, 41), cmapRdBu_r, transformccrs.PlateCarree()) # 叠加等高线例如位势高度场 contour ax.contour(lons, lats, height, levelsnp.arange(1200, 1600, 40), colorsblack, linewidths1, transformccrs.PlateCarree()) ax.clabel(contour, inlineTrue, fontsize10, fmt%d) plt.colorbar(contourf, axax, orientationhorizontal, pad0.05, labelRelative Vorticity (10^-5 /s)) ax.set_title(f850hPa Relative Vorticity and Geopotential Height at {ds.time[time_idx].values}) plt.show()这段代码绘制了涡度填色和位势高度等值线。你应该能看到正涡度区暖色与位势高度低值区气旋有很好的对应关系负涡度区冷色对应高压脊。这是对计算结果最直观的物理验证。5. 高级话题精度、效率与常见陷阱当你成功跑出第一张图后我们深入聊聊那些影响结果可靠性和计算效率的细节。这些是我在长期使用中积累的经验很多在官方文档里不会强调。5.1 差分格式的选择与影响前面提到metpy默认用中心差分。但在某些场景下你可能需要自己实现差分。比如你的数据分辨率非常高如对流尺度模式格距小于5公里中心差分可能放大小尺度噪音。这时可以考虑平滑滤波后再计算或者使用更高阶的差分格式如谱方法或紧致差分但这超出了metpy的内置功能需要借助scipy或自己实现。一个简单的平滑方法是使用scipy.ndimage的高斯滤波from scipy import ndimage # 对风场进行轻微平滑sigma1格点 u_smooth ndimage.gaussian_filter(u.values, sigma1.0) v_smooth ndimage.gaussian_filter(v.values, sigma1.0) # 将平滑后的数组重新包装成带单位的DataArray u_smoothed xr.DataArray(u_smooth, dimsu.dims, coordsu.coords).metpy.quantify() v_smoothed xr.DataArray(v_smooth, dimsv.dims, coordsv.coords).metpy.quantify() # 再用平滑后的风场计算涡度 vort_smoothed mpcalc.vorticity(u_smoothed, v_smoothed)平滑会损失一些真实的微小尺度特征但能有效抑制由计算噪声产生的虚假小涡旋使大尺度特征更清晰。在天气尺度分析中适度的平滑通常是可接受的。5.2 极地和高纬度地区的特殊处理这是球坐标计算的最大挑战。在北极点纬度90°N所有经线汇合cosφ趋近于0导致公式分母为零。即使不在极点在很高纬度如80°N以上cosφ非常小使得经向格距Δx变得极小微小的风场差分误差会被极度放大产生无意义的巨大涡度/散度值。metpy的函数内部通过数值方法避免除以零但在极高纬度结果仍然不可信。实践中的黄金法则是通常只分析60°S到60°N之间的中低纬度区域。如果你必须分析极地应考虑将数据投影到极射赤面投影如ccrs.NorthPolarStereo下的笛卡尔坐标然后在投影后的直角坐标网格上计算涡度和散度。这涉及到数据重投影和插值是一个更复杂的话题可以使用cartopy的transform_points功能结合scipy的插值来实现。5.3 计算性能优化处理三维时空数据上面的例子是针对单个时次、单层的数据。实际研究中我们常需要处理包含几十个层次、上百个时次的数据集。循环调用mpcalc函数虽然简单但效率可能不高。优化策略是向量化和使用xarray的apply_ufunc。metpy的许多函数支持xarray的apply_ufunc这允许你将函数应用到多维数组的指定维度上而无需显式循环。import xarray as xr import metpy.calc as mpcalc # 假设ds是一个包含多个时间和层次的数据集 # 计算所有时间和层次上的涡度 vorticity_4d xr.apply_ufunc( mpcalc.vorticity, ds.u, ds.v, input_core_dims[[latitude, longitude], [latitude, longitude]], output_core_dims[[latitude, longitude]], vectorizeTrue, # 重要允许循环处理非核心维度如time, level daskparallelized, # 如果使用dask数组可以并行计算 output_dtypes[np.float64] )这段代码会一次性计算出整个数据立方体时间、层次、纬度、经度的涡度场。input_core_dims指定了函数实际处理的维度这里是经纬度apply_ufunc会自动在其他维度时间、层次上循环。vectorizeTrue是关键它告诉xarray对非核心维度进行循环调用。如果数据是分块的dask数组设置daskparallelized可以利用多核进行并行计算大幅提升处理大量数据的速度。5.4 单位制的陷阱与一致性metpy的单元管理非常强大但有时也会带来困惑。一个常见错误是计算出的物理量单位看起来很奇怪比如涡度平流的单位是meter / second ** 2而不是预期的/s^2。这是因为在计算过程中距离单位米被保留了。你需要确保在计算链的每一步单位都符合你的预期。使用.to()方法进行转换是最安全的。另一个陷阱是无量纲化。有些模式输出或再分析数据其涡度、散度可能已经做了尺度化例如除以科氏参数f。在计算平流时你必须使用原始的、有物理单位的量。始终用print(data_array.units)检查中间结果的单位。6. 结果解读与天气学应用让数据“说话”计算出漂亮的场只是第一步更重要的是理解它们在天气图上的意义并用于诊断分析。这里结合几个典型天气系统说说怎么看图。6.1 涡度平流与槽脊发展在500hPa等压面图上正涡度平流区暖色通常位于高空槽前负涡度平流区冷色位于脊前。根据准地转理论正涡度平流导致地面气压下降气旋发展负涡度平流导致地面气压上升反气旋发展。因此你可以在500hPa图上画出强正涡度平流中心。向下游方向大致沿气流方向寻找地面天气图。通常会发现该正涡度平流中心的下方或略偏下游对应着一个发展中的地面低压或低压槽。这就是利用涡度平流预报气旋发展的基本原理。6.2 温度平流与锋面、垂直运动温度平流直接关联到冷暖空气的输送。暖平流区Adv_T 0通常出现在高空槽前低层暖舌区域。暖平流导致局地增温根据热成风关系会产生涡度的垂直变化进而强迫出上升运动。因此强暖平流区是潜在的对流和降水区。冷平流区Adv_T 0通常出现在高空槽后脊前。冷平流导致局地降温强迫出下沉运动天气往往晴朗。在850hPa或700hPa层面分析温度平流尤其有用。你可以叠加风场和等温线风从冷区吹向暖区且与等温线有交角的地方就是明显的冷平流区反之则为暖平流区。计算出的温度平流场可以定量地验证你的定性判断并找出平流最强的核心区域。6.3 散度场与垂直运动的间接推断在摩擦层以上大尺度运动的散度通常很小10^-5 /s量级但其垂直分布至关重要。根据连续方程低层辐合D0对应高层辐散D0时其间必有强烈的上升运动。单独看某一层的散度意义不大需要看整层积分或垂直剖面。一个实用的近似是850hPa的辐合中心往往与降水区有较好的对应特别是当该处同时有暖平流和正涡度平流时上升运动的三重条件就具备了很可能发生强降水。你可以尝试将850hPa散度、700hPa温度平流、500hPa涡度平流叠加在一张综合诊断图上这样天气系统的动力、热力结构就一目了然。6.4 一个综合诊断案例识别气旋发展区假设你有一组完整的多层数据可以执行以下步骤进行综合诊断计算各层物理量对850hPa、700hPa、500hPa分别计算涡度、散度、温度平流、涡度平流。绘制叠加图用填色图显示500hPa涡度平流。用等值线显示500hPa位势高度场看槽脊位置。用矢量箭头显示850hPa风场看低层气流。用另一种填色或等值线显示850hPa温度平流。分析找到500hPa正涡度平流最大中心。看其下游方向850hPa是否有明显的暖平流和风场辐合。如果两者重叠或紧密相邻这个区域就是气旋发生发展的关键区PVA最大值下游、暖平流最大值上空。再结合比湿场或抬升凝结高度就能进一步判断降水潜势。通过这样的定量计算和叠加分析你对天气系统的理解就从“看图说话”进入了“物理诊断”的层次预报和研究的底气都会足很多。7. 代码封装与批量处理实战最后我们把散落的代码封装成一个可复用的函数并演示如何批量处理多个时次的数据输出为新的netCDF文件方便后续分析和绘图。import xarray as xr import numpy as np import metpy.calc as mpcalc from metpy.units import units import warnings warnings.filterwarnings(ignore) # 忽略一些metpy的单位警告 def calculate_diagnostic_fields(ds, level_hPa850): 计算单层数据集的涡度、散度、平流场。 参数: ds (xarray.Dataset): 包含u, v, t, z变量的数据集。 level_hPa (int/float): 需要计算的等压面层次单位hPa。 返回: xarray.Dataset: 包含原始变量和新增诊断变量的数据集。 # 1. 选择层次并附加单位 level level_hPa * units.hPa try: u ds[u].sel(levellevel).metpy.quantify() v ds[v].sel(levellevel).metpy.quantify() t ds[t].sel(levellevel).metpy.quantify() z ds[z].sel(levellevel).metpy.quantify() except KeyError as e: # 如果变量名不同尝试常见别名 var_map {u: U, v: V, t: T, z: Z, level: pressure} for old, new in var_map.items(): if old in ds.dims or old in ds.coords: ds ds.rename({old: new}) u ds[U].sel(pressurelevel).metpy.quantify() v ds[V].sel(pressurelevel).metpy.quantify() t ds[T].sel(pressurelevel).metpy.quantify() z ds[Z].sel(pressurelevel).metpy.quantify() # 2. 计算诊断量 print(f计算 {level_hPa}hPa 的相对涡度...) vort mpcalc.vorticity(u, v) print(f计算 {level_hPa}hPa 的散度...) div mpcalc.divergence(u, v) print(f计算 {level_hPa}hPa 的涡度平流...) vort_adv mpcalc.advection(vort, [u, v]) print(f计算 {level_hPa}hPa 的温度平流...) temp_adv mpcalc.advection(t, [u, v]) # 3. 转换为常用单位并剥离单位便于存储为netCDF vort_common vort.to(10^-5/s).metpy.dequantify() div_common div.to(10^-5/s).metpy.dequantify() vort_adv_common vort_adv.to(10^-9/s^2).metpy.dequantify() temp_adv_common temp_adv.to(10^-5 K/s).metpy.dequantify() # 4. 创建包含结果的新Dataset result_ds xr.Dataset({ u: u.metpy.dequantify(), v: v.metpy.dequantify(), t: t.metpy.dequantify(), z: z.metpy.dequantify(), relative_vorticity: vort_common, divergence: div_common, vorticity_advection: vort_adv_common, temperature_advection: temp_adv_common }) # 复制全局属性 result_ds.attrs ds.attrs.copy() result_ds.attrs[diagnostic_calculated] fusing MetPy on {level_hPa}hPa return result_ds def batch_process_era5(input_file_pattern, output_dir, levels[850, 700, 500]): 批量处理多个ERA5文件。 import glob import os # 找到所有输入文件 input_files sorted(glob.glob(input_file_pattern)) print(f找到 {len(input_files)} 个文件待处理。) for i, f in enumerate(input_files): print(f处理文件 ({i1}/{len(input_files)}): {os.path.basename(f)}) ds xr.open_dataset(f) # 对每个需要计算的层次进行处理 for lev in levels: print(f 处理 {lev}hPa 等压面...) result_ds calculate_diagnostic_fields(ds, lev) # 构建输出文件名 basename os.path.basename(f).replace(.nc, ) output_file os.path.join(output_dir, f{basename}_diagnostic_{lev}hPa.nc) # 保存为netCDF encoding {var: {zlib: True, complevel: 5} for var in result_ds.data_vars} result_ds.to_netcdf(output_file, encodingencoding) print(f 结果已保存至: {output_file}) ds.close() # 使用示例 if __name__ __main__: # 假设你的ERA5文件命名类似 era5_20230101.nc batch_process_era5( input_file_pattern./era5_data/era5_*.nc, output_dir./diagnostic_results/, levels[850, 500] # 只计算850和500hPa两层以节省时间 )这个脚本提供了两个核心函数。calculate_diagnostic_fields函数封装了单层计算的所有步骤包括单位处理和常见变量名的兼容。batch_process_era5函数则展示了如何批量处理多个文件和多层数据并将结果保存为新的netCDF文件压缩存储以节省空间。运行这个脚本你就可以将原始风温压场数据一键转化为包含涡度、散度、平流等诊断量的、可直接用于绘图和分析的数据集。这极大地提升了工作效率让你能更专注于物理过程的分析而不是重复的数据处理劳动。本文还有配套的精品资源点击获取