Python | ERA5任意经纬度风矢量剖面绘制

用户11172986

发布于 2025-01-20 20:29:56

9900

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

Python | ERA5任意经纬度风矢量剖面绘制

摘要

既然做了wrf的风剖面绘制，那么ERA5也一起做了当然还是基于metpy的函数

导入必要的库

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.interpolate import cross_section
from matplotlib.colors import Normalize

数据读取


# 加载 ERA5 数据
file_path = '/home/mw/input/era58091/ERA5-2023-08_pl.nc'
ds = xr.open_dataset(file_path)

###  定义剖面起点和终点

# 定义剖面起点和终点（经纬度）
start_point = (45, 100.0)  # (lat, lon)
end_point = (20, 130.0)    # (lat, lon)

### 5. 使用 MetPy 提取剖面数据

# 取第一个时刻，将数据转换为 MetPy 的 CF 标准格式
ds = ds.isel(time=0).metpy.parse_cf()

# 使用 MetPy 的 cross_section 函数提取剖面数据
cross = cross_section(ds, start_point, end_point)
print(cross)

<xarray.Dataset> Size: 210kB
Dimensions:    (level: 37, index: 100)
Coordinates:
  * level      (level) int32 148B 1 2 3 5 7 10 20 ... 875 900 925 950 975 1000
    time       datetime64[ns] 8B 2023-08-02
    metpy_crs  object 8B Projection: latitude_longitude
    longitude  (index) float64 800B 100.0 100.4 100.8 ... 129.5 129.8 130.0
    latitude   (index) float64 800B 45.0 44.79 44.59 44.38 ... 20.56 20.28 20.0
  * index      (index) int64 800B 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables:
    z          (level, index) float64 30kB 4.752e+05 4.751e+05 ... 445.0 507.9
    r          (level, index) float64 30kB -0.0005798 -0.0005798 ... 87.7 89.32
    t          (level, index) float64 30kB 259.9 259.9 260.1 ... 299.6 299.3
    q          (level, index) float64 30kB 4.048e-06 4.048e-06 ... 0.01905 0.019
    u          (level, index) float64 30kB -39.16 -38.9 -38.8 ... 8.306 7.847
    v          (level, index) float64 30kB 2.141 2.308 2.252 ... 13.12 13.05
    w          (level, index) float64 30kB -0.00102 -0.001254 ... 0.01719
Attributes:
    Conventions:  CF-1.6
    history:      2024-01-09 16:38:26 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

计算切向与法向风

import metpy.calc as mpcalc
### 6. 计算沿剖面的风场数据
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(
    cross['u'],
    cross['v']
)

计算W(m/s)

cross['vertical_velocity'] = mpcalc.vertical_velocity(
    cross['w'],
    cross['level'],
    cross['t']
)

绘图

import matplotlib.pyplot as plt
import numpy as np

# 设置LaTeX字体
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman']  # 使用Times New Roman字体

# 定义图形对象和主坐标轴
fig, ax = plt.subplots(1, 1, figsize=(12., 8.), dpi=200)

# 绘制相对湿度的等高线图
rh_contour = ax.contourf(cross['longitude'], cross['level'], cross['vertical_velocity'],
                         levels=10, cmap='YlGnBu')
rh_colorbar = fig.colorbar(rh_contour)
rh_colorbar.set_label(r'$\omega$ (m/s)')  # 使用LaTeX格式的标签

# 绘制风场箭头，使用自定义索引以减少箭头的密度
wind_slc_vert = list(range(0, 19, 2)) + list(range(19, 29))
wind_slc_horz = slice(5, 100, 5)
ax.quiver(cross['longitude'][::2], cross['level'][::2],
          cross['n_wind'][::2,::2],
          cross['vertical_velocity'][::2,::2]* 100, color='k')

# 设置y轴反向（气压递减）
ax.invert_yaxis()

# 设置标题和标签
ax.set_title(r'Vertical Velocity ($\omega$) and Wind Vectors', fontsize=16)
ax.set_ylabel(r'Pressure (hPa)', fontsize=14)
ax.set_xlabel(r'Longitude (degrees east)', fontsize=14)

# 设置刻度字体大小
ax.tick_params(axis='both', which='major', labelsize=12)

# 增加网格线
ax.grid(True, linestyle='--', alpha=0.6)

# 显示图形
plt.show()