保姆级教程:用Python+Cartopy绘制专业气象图(以ERA5 500hPa位势高度场为例)
Python气象可视化实战:从ERA5数据到出版级位势高度场
当气象数据遇上Python的Cartopy库,枯燥的数字矩阵就能蜕变为直观的专业图表。对于气象、地理或环境科学领域的研究者而言,掌握这套工具链意味着能够自主完成从原始数据到发表级图形的全流程工作。本文将以ERA5再分析数据中的500hPa位势高度场为例,手把手带你构建一个可复用的气象绘图模板。
1. 环境配置与数据准备
工欲善其事,必先利其器。在开始绘图前,我们需要搭建一个稳定的Python环境并获取正确的数据集。推荐使用Anaconda创建专属的气象分析环境:
conda create -n weather_plot python=3.9 conda activate weather_plot conda install -c conda-forge cartopy matplotlib xarray numpyERA5数据可以通过Copernicus Climate Data Store获取,这里我们选择2023年夏季(6-8月)北半球500hPa位势高度场月平均数据。下载时需注意:
- 时间范围:2023年6月1日至8月31日
- 空间范围:经度0-360°,纬度90°N-0°
- 变量选择:Geopotential (z)
数据下载完成后,建议按照以下结构组织项目目录:
/project /data ERA5_500hPa_2023_JJA.nc /scripts plot_geopotential.py /output2. 数据处理与夏季平均计算
ERA5数据通常以NetCDF格式存储,xarray库为此提供了完美的解决方案。与传统的numpy数组不同,xarray保留了维度信息,使得数据操作更加直观。
import xarray as xr # 读取数据并选择东亚区域(60-140°E, 10-70°N) ds = xr.open_dataset('data/ERA5_500hPa_2023_JJA.nc') hgt = ds['z'].sel(longitude=slice(60,140), latitude=slice(70,10)) # 计算夏季平均(注意单位转换) hgt_summer_avg = hgt.mean(dim='time') / 9.8 # 转换为位势米处理过程中常见的坑点:
- 经度范围表示:ERA5使用0-360°,而Cartopy默认-180-180°
- 单位一致性:ERA5的位势高度单位为m²/s²,需除以重力加速度9.8
- 内存管理:大数据集处理时建议使用chunk分块读取
3. Cartopy绘图核心框架搭建
Cartopy的强大之处在于其丰富的地图投影系统和地理特征支持。我们先构建一个基础绘图框架:
import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 创建图形和坐标轴 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # 设置地图范围 ax.set_extent([60, 140, 10, 70], crs=ccrs.PlateCarree()) # 添加基本地理要素 ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8) ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='lightgray') ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='lightblue')专业气象图需要精细的经纬网格标注,Cartopy提供了灵活的网格线控制:
# 配置经纬网格 gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') gl.top_labels = False gl.right_labels = False gl.xlabel_style = {'size': 10} gl.ylabel_style = {'size': 10}4. 位势高度场等值线绘制技巧
等值线是位势高度场的核心表达方式,matplotlib的contour函数提供了丰富的控制参数:
import numpy as np # 生成等值线 levels = np.arange(5600, 6000, 40) cs = ax.contour(hgt_summer_avg.longitude, hgt_summer_avg.latitude, hgt_summer_avg, levels=levels, colors='black', linewidths=1.2) # 添加等值线标签 ax.clabel(cs, fmt='%1.0f', fontsize=10) # 特别突出588线 cs_588 = ax.contour(hgt_summer_avg.longitude, hgt_summer_avg.latitude, hgt_summer_avg, levels=[5880], colors='red', linewidths=2.5) ax.clabel(cs_588, fmt='5880', fontsize=12)专业气象图的视觉优化要点:
- 等值线间隔:通常选择40位势米间隔
- 关键线突出:588线用红色加粗表示副热带高压
- 标签清晰度:避免标签重叠,适当调整inline_spacing
5. 出版级图形优化与输出
最后的视觉调整往往决定图形的专业程度:
# 添加标题和注释 plt.title('2023年夏季500hPa位势高度场(单位:位势米)', fontsize=14, pad=20) # 添加比例尺和指北针 ax.add_artist(ScaleBar(ax)) ax.add_artist(NorthArrow(ax)) # 保存高分辨率图片 plt.savefig('output/500hPa_geopotential.png', dpi=600, bbox_inches='tight', facecolor='white')输出格式选择建议:
- 期刊投稿:TIFF或PDF格式,分辨率≥300dpi
- 网页展示:PNG格式,适当降低分辨率
- 矢量图形:SVG或PDF格式,方便后期编辑
6. 进阶技巧与模块化封装
将绘图流程封装成函数可以提高代码复用率:
def plot_geopotential(data, region, levels, highlight_lines=None): """绘制位势高度场 参数: data -- xarray DataArray格式的位势高度数据 region -- 绘图区域[lon_min, lon_max, lat_min, lat_max] levels -- 等值线层级数组 highlight_lines -- 需要突出的等值线值列表 """ fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # 绘图核心逻辑 # ... return fig, ax其他实用进阶技巧:
- 多子图布局:比较不同季节或年份的形势
- 叠加风场:用quiver叠加风矢量
- 地形效应:添加地形阴影relief效果
- 交互式探索:结合ipywidgets创建交互界面
掌握这些技术后,你可以轻松应对各种气象可视化需求,从日常分析到学术发表,Python+Cartopy的组合将成为你的得力助手。在实际项目中,我习惯将常用配置保存为模板文件,新项目只需调整数据源和少量参数即可快速生成专业图表。
