首页 > Python资料 博客日记
Python纬向垂直剖面图绘制
2025-01-16 09:30:06Python资料围观30次
纬向垂直剖面图绘制
文章目录
以上风羽为U-V风； - -为大气温度；等值面填色图为 ω \omega ω。
一、数据内容
根据上一期切割的NC数据，内容如下：
- 垂直速度omega
- 大气温度air
- u方向风uwnd
- v方向风vwnd（新增v方向风的文件）
在xarray中得到的各要素内容为：
>>>data
<xarray.Dataset>
Dimensions: (level: 12, lon: 13)
Coordinates:
* level (level) float32 1e+03 925.0 850.0 700.0 ... 250.0 200.0 150.0 100.0
* lat float32 ...
* time datetime64[ns] ...
* lon (lon) float32 110.0 112.5 115.0 117.5 ... 132.5 135.0 137.5 140.0
Data variables:
omega (level, lon) float32 ...
air (level, lon) float32 ...
uwnd (level, lon) float32 ...
>>>vdata
<xarray.Dataset>
Dimensions: (valid_time: 1, pressure_level: 12, latitude: 1, longitude: 121)
Coordinates:
* valid_time (valid_time) datetime64[ns] 2018-09-15
* pressure_level (pressure_level) float64 1e+03 925.0 850.0 ... 150.0 100.0
* latitude (latitude) float64 17.5
* longitude (longitude) float64 110.0 110.2 110.5 ... 139.5 139.8 140.0
Data variables:
v (valid_time, pressure_level, latitude, longitude) float32 ...
Attributes:
GRIB_centre: ecmf
GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts
GRIB_subCentre: 0
Conventions: CF-1.7
institution: European Centre for Medium-Range Weather Forecasts
history: 2024-10-23T00:56 GRIB to CDM+CF via cfgrib-0.9.1...
二、所用到的库与数据读取
import datetime
import numpy as np
import matplotlib.pyplot as plt
from metpy.units import units
import metpy.constants as constants
import metpy.calc as mpcalc
import xarray as xr
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
file = 'F:/data/cut_example_data.nc'
data = xr.open_dataset(file)
vfile = 'F:/data/vwnd_example.nc'
vdata = xr.open_dataset(vfile)
# 对数据数组的维度进行重排（transpose）。希望所有变量的维度按照，时间-高度-纬度-经度 排序。
# 如果对自己读取数据的维度顺序满意，就不需要这条。
data = data.transpose('time','level','lat','lon')
time = data.time
# 等压面，赋予其对应单位
level = data.level * units.hPa
lon = data.lon
lat = data.lat
# P坐标系下垂直速度，赋予对应单位
omega = data.omega * units('Pa/s')
# 将开尔文温标转换成摄氏度，并赋予对应单位
air = (data.air - 273.15) * units.degC
uwnd = data.uwnd
vwnd = vdata.v
三、绘图
fig = plt.figure(figsize=(12, 6))
ax = plt.subplot(111)
（1）垂直速度填色图：
contour_color = ['white','#b5c9ff','#7f96ff','#009e1e','#b9f96e','#ffa309','#810000'][::-1]
contour_clevs = [0,-1,-2,-3,-4,-5,-6,-7][::-1]
# omega填色图
ct = ax.contourf(lon, level, omega[0,:,0,:] * 10, colors=contour_color, levels=contour_clevs)
cb = plt.colorbar(ct,shrink=0.9)
# 重设色标刻度显示
cb.ax.set_yticklabels([''] + [str(tick) for tick in cb.get_ticks()[1:-1]] + [''])
# 色标标题
cb.ax.set_title(r'$\omega}$', fontsize=15)
（2）温度等值线图：
clevs_t = np.arange(-40, 20, 8)
air_contour = ax.contour(lon, level, air[0,:,0,:], colors='r', levels=clevs_t, linestyles='--', alpha=0.7)
ax.clabel(air_contour, clevs_t, inline=True, fmt='%r ℃')
（3）u-v风风羽图
## 设置风羽，半杠2m/s，一杠4m/s，箭头20m/s
ax.barbs(lon, level, uwnd[0,:,0,:], vwnd.loc[time,:,lat,lon][0,:,0,:], alpha=0.7,
barb_increments={'half': 2, 'full': 4, 'flag': 20})
（4）经纬标题等设置
# 采用latex方法设置标题
t = ax.set_title(r'Vertical Velocity($10^{-1}Pa/s$) & Horizontal Wind$(m/s)$ & Temp$(\degree C)$',
loc='right', fontsize=15)
# y轴的标题
ax.set_ylabel('Pressure(hPa)',fontsize=15)
# 将坐标轴倒置，使1000hPa从底下开始
plt.gca().invert_yaxis()
（5）完整代码
import datetime
import numpy as np
import matplotlib.pyplot as plt
from metpy.units import units
import metpy.constants as constants
import metpy.calc as mpcalc
import xarray as xr
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
# 参考资料
# https://blog.csdn.net/weixin_44237337/article/details/124759780
file = 'F:/data/cut_example_data.nc'
data = xr.open_dataset(file)
vfile = 'F:/data/vwnd_example.nc'
vdata = xr.open_dataset(vfile)
# 对数据数组的维度进行重排（transpose）。希望所有变量的维度按照，时间-高度-纬度-经度 排序。
# 如果对自己读取数据的维度顺序满意，就不需要这条。
data = data.transpose('time','level','lat','lon')
time = data.time
# 等压面，赋予其对应单位
level = data.level * units.hPa
lon = data.lon
lat = data.lat
# P坐标系下垂直速度，赋予对应单位
omega = data.omega * units('Pa/s')
# 将开尔文温标转换成摄氏度，并赋予对应单位
air = (data.air - 273.15) * units.degC
uwnd = data.uwnd
vwnd = vdata.v
fig = plt.figure(figsize=(12, 6))
ax = plt.subplot(111)
# ………………………………………………………垂直剖面图………………………………………………………………………………………………………………#
# 我是对着参考图的色标从上到下抄的，但是绘制填色图要求要递增。懒得重抄，加了[::-1]倒叙。
contour_color = ['white','#b5c9ff','#7f96ff','#009e1e','#b9f96e','#ffa309','#810000'][::-1]
contour_clevs = [0,-1,-2,-3,-4,-5,-6,-7][::-1]
# omega填色图
ct = ax.contourf(lon, level, omega[0,:,0,:] * 10, colors=contour_color, levels=contour_clevs)
cb = plt.colorbar(ct,shrink=0.9)
cb.ax.set_yticklabels([''] + [str(tick) for tick in cb.get_ticks()[1:-1]] + [''])
cb.ax.set_title(r'$\omega}$', fontsize=15)
# 等温线图
clevs_t = np.arange(-40, 20, 8)
air_contour = ax.contour(lon, level, air[0,:,0,:], colors='r', levels=clevs_t, linestyles='--', alpha=0.7)
ax.clabel(air_contour, clevs_t, inline=True, fmt='%r ℃')
# u-v风羽图
## 索引到的v风建议还需要压缩维度，所以直接选取time第0个和lat第0个
## 设置风羽，半杠2m/s，一杠4m/s，箭头20m/s
ax.barbs(lon, level, uwnd[0,:,0,:], vwnd.loc[time,:,lat,lon][0,:,0,:], alpha=0.7,
barb_increments={'half': 2, 'full': 4, 'flag': 20})
t = ax.set_title(r'Vertical Velocity($10^{-1}Pa/s$) & Horizontal Wind$(m/s)$ & Temp$(\degree C)$',
loc='right', fontsize=15)
ax.set_ylabel('Pressure(hPa)',fontsize=15)
plt.gca().invert_yaxis()
plt.show()
探讨：u-w合成风绘制的垂直风场图
首先在使用 w w w 风前，我们需要将 P P P坐标系下的 ω \omega ω转换。 P a / s — > m / s Pa/s—>m/s Pa/s—>m/s。
（1）ω 转 w 详解。
metpy 库配备了可以灵活处理常量和变量的计算公式，采用公式：
metpy.calc.vertical_velocity(omega,pressure, temperature, mixing_ratio=0)
可以将 ω 转 w 。参考自python 将垂直速度从压力坐标系转为高度坐标系（pa/s转为m/s）_垂直速度单位转换-CSDN博客
根据官方文档
metpy.calc.vertical_velocity(omega,pressure, temperature, mixing_ratio=0)
Calculate
wfromωassuming hydrostatic conditions.假定流体静力学的条件，根据
ω计算w。
This function converts vertical velocity with respect to pressure ( ω = D p D t ) \left(\omega = \frac{Dp}{Dt}\right) (ω=DtDp) to that with respect to height ( w = D z D t ) \left(w = \frac{Dz}{Dt}\right) (w=DtDz)assuming hydrostatic conditions on the synoptic scale. By Equation 7.33 in [Hobbs2006],假定流体静力学的条件，该函数将垂直速度与压力 ( ω = D p D t ) \left(\omega = \frac{Dp}{Dt}\right) (ω=DtDp)的关系转换为垂直速度与高度 ( w = D z D t ) \left(w = \frac{Dz}{Dt}\right) (w=DtDz)的关系
根据 ω ≃ − ρ g w \omega \simeq -\rho g w ω≃−ρgw ，得 w ≃ − ω ρ g w \simeq \frac{- \omega}{\rho g} w≃ρg−ω 。
Density ρ \rho ρ is calculated using thedensity()function, from the given pressure and temperature. Ifmixing_ratiois given, the virtual temperature correction is used, otherwise, dry air is assumed.在给定气压和温度下，计算密度
ρ可以用density()函数。如果mixing_ratio给出，会借助虚温订正，否则，将认为是干空气。
引用自metpy 库中记录的函数 vertical_velocity() ：
def vertical_velocity(omega, pressure, temperature, mixing_ratio=0):
"""
balabalabala……
mixing_ratio: `pint.Quantity`, optional
Mixing ratio of air 空气混合比例
"""
rho = density(pressure, temperature, mixing_ratio)
return (omega / (- mpconsts.g * rho)).to('m/s')
def density(pressure, temperature, mixing_ratio, molecular_weight_ratio=mpconsts.epsilon):
"""
balabalabala……
mixing_ratio : `pint.Quantity`
Mass mixing ratio (dimensionless) 无量纲质量混合比
"""
virttemp = virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio)
return (pressure / (mpconsts.Rd * virttemp)).to('kg/m**3')
def virtual_temperature(
temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon
):
"""
balabalabala……
mixing_ratio : `pint.Quantity`
Mass mixing ratio (dimensionless) 无量纲质量混合比
molecular_weight_ratio : `pint.Quantity` or float, optional
The ratio of the molecular weight of the constituent gas to that assumed
for air. Defaults to the ratio for water vapor to dry air.
组成气体的分子量与假设的空气分子量之比。默认为水蒸气与干燥空气的比率。
(:math:`\epsilon\approx0.622`)
"""
return temperature * ((mixing_ratio + molecular_weight_ratio)
/ (molecular_weight_ratio * (1 + mixing_ratio)))
根据以上函数vertical_velocity() 内容及展开，我们可以将公式解析为：
w
=
−
ω
g
∗
ρ
w=-\frac{\omega}{g*\rho}
w=−g∗ρω
其中密度
ρ
=
P
R
d
∗
T
v
\rho = \frac{P}{R_d*T_v}
ρ=Rd∗TvP；
T
v
T_v
Tv为虚温，
T
v
=
T
∗
w
m
+
ϵ
ϵ
∗
(
1
+
w
m
)
T_v = T*\frac{w_m+\epsilon}{\epsilon*(1+w_m)}
Tv=T∗ϵ∗(1+wm)wm+ϵ;g为重力加速度常数：
9.80665
m
/
s
2
9.80665m/s^2
9.80665m/s2；
R
d
R_d
Rd——mpconsts.Rd，dry_air_gas_constant，地表干空气气体常数：
287.047491
J
/
(
K
⋅
k
g
)
287.047491J/(K·kg)
287.047491J/(K⋅kg) 。
在虚温的公式中，
w
m
w_m
wm为质量混合比，
w
m
=
m
i
m
d
w_m=\frac{m_i}{m_d}
wm=mdmi，这里默认值为 0 ；
ϵ
\epsilon
ϵ——mpconsts.nounit.epsilon，molecular_weight_ratio，水的分子量与干燥空气的分子量之比：
0.62195691
0.62195691
0.62195691。
综上，笛卡尔坐标系的垂直速度，其运算公式可以推导为：
w
=
−
ω
∗
T
P
∗
R
d
∗
(
w
m
+
ϵ
)
g
∗
ϵ
∗
(
1
+
w
m
)
w=-\frac{\omega*T}{P}*\frac{R_d*(w_m+\epsilon)}{g*\epsilon*(1+w_m)}
w=−Pω∗T∗g∗ϵ∗(1+wm)Rd∗(wm+ϵ)
w = − ω ∗ T P ∗ R d g = − ω ∗ T P ∗ 29.27069 w=-\frac{\omega*T}{P}*\frac{R_d}{g}=-\frac{\omega*T}{P}*29.27069 w=−Pω∗T∗gRd=−Pω∗T∗29.27069
根据官方文档，其算法参考自
Hobbs, P. V., and J. M. Wallace, 2006: Atmospheric Science: An Introductory Survey. 2nd ed. Academic Press, 504 pp.
感兴趣的可以看一下，全英的。
那么在代码中的引用则是（需要使用的变量记得要先赋予其对应单位）：
time = data.time
# 等压面，赋予其对应单位
level = data.level * units.hPa
lon = data.lon
lat = data.lat
# P坐标系下垂直速度，赋予对应单位
omega = data.omega * units('Pa/s')
# 将开尔文温标转换成摄氏度，并赋予对应单位
air = (data.air - 273.15) * units.degC
uwnd = data.uwnd
# 笛卡尔坐标系下的垂直速度
w = mpcalc.vertical_velocity(omega, level, air, mixing_ratio=0)
而对应的，其纬向垂直速度剖面图如下：
# omega填色图
contour_color = ['white','#b5c9ff','#7f96ff','#009e1e','#b9f96e','#ffa309','#810000'][::-1]
contour_clevs = [0,-1,-2,-3,-4,-5,-6,-7][::-1]
# omega填色图
ct = ax.contourf(lon, level, omega[0,:,0,:] * 10, colors=contour_color, levels=contour_clevs)
cb = plt.colorbar(ct,shrink=0.9)
cb.ax.set_yticklabels([''] + [str(tick) for tick in cb.get_ticks()[1:-1]] + [''])
cb.ax.set_title(r'$\omega}$', fontsize=15)
# w等值线图
# 过滤小于0的值
w_filtered = w[0,:,0,:].where(w[0,:,0,:]>=0,0)
w_contour = ax.contour(lon, level, w_filtered, colors='black', linestyles='--')
ax.clabel(w_contour, inline=True, fmt='%2.2f')
（2）垂直风速尺度与水平风尺度对比
根据《天气学原理和方法》第四版, Page629。
ω = 1 P a / s \omega=1Pa/s ω=1Pa/s时 w w w 对应的值：
| 1000hPa | -0.086m/s |
|---|---|
| 850hPa | -0.094m/s |
| 700hPa | -0.108m/s |
| 500hPa | -0.137m/s |
从量级上看，水平方向上的风远大于垂直方向上的风，量级至少是十倍以上。
所以一般用于绘制垂直风速风场的剖面图时，会将垂直风速的量级放大，放大到100倍左右。
因此常用的绘制方法是，
w
=
−
ω
∗
100
w=-\omega*100
w=−ω∗100。则合成如下：
而为了单位的统一性，用metpy库将
ω
\omega
ω转换
w
w
w，再对量级放大100倍后绘制，得到的结果效果反而并不明显。
这里需要注意的是：
垂直速度 w = w ∗ 100 ∗ 1 0 − 2 m / s w=w*100*10^{-2}m/s w=w∗100∗10−2m/s ，即单位： 1 0 − 2 m / s 10^{-2}m/s 10−2m/s
合成风U-W= ( w ∗ 1 0 − 2 m / s ) 2 + ( u ∗ m / s ) 2 \sqrt{(w*10^{-2}m/s)^2+(u*m/s)^2} (w∗10−2m/s)2+(u∗m/s)2，即单位： 1 0 − 4 + 1 ∗ m / s ≃ 1.00005 ∗ m / s ≈ m / s \sqrt{10^{-4}+1}*m/s\simeq1.00005*m/s\approx m/s 10−4+1∗m/s≃1.00005∗m/s≈m/s
因此，按照垂直速度放大100倍进行绘图，对最后的结果并不会有太大的影响。但是，如果有需要绘制该类型的风场图，更建议直接采用 ω \omega ω放大100倍来进行绘制。
标签：
相关文章
最新发布
- 光流法结合深度学习神经网络的原理及应用（完整代码都有Python opencv）
- Python 图像处理进阶：特征提取与图像分类
- 大数据可视化分析-基于python的电影数据分析及可视化系统_9532dr50
- 【Python】入门（运算、输出、数据类型）
- 【Python】第一弹---解锁编程新世界：深入理解计算机基础与Python入门指南
- 华为OD机试E卷 --第k个排列 --24年OD统一考试（Java & JS & Python & C & C++）
- Python已安装包在import时报错未找到的解决方法
- 【Python】自动化神器PyAutoGUI —告别手动操作，一键模拟鼠标键盘，玩转微信及各种软件自动化
- Pycharm连接SQL Sever（详细教程）
- Python编程练习题及解析（49题）
点击排行
- 版本匹配指南：Numpy版本和Python版本的对应关系
- 版本匹配指南：PyTorch版本、torchvision 版本和Python版本的对应关系
- Python 可视化 web 神器：streamlit、Gradio、dash、nicegui；低代码 Python Web 框架：PyWebIO
- 相关性分析——Pearson相关系数+热力图(附data和Python完整代码)
- Anaconda版本和Python版本对应关系（持续更新...）
- Python与PyTorch的版本对应
- Windows上安装 Python 环境并配置环境变量 （超详细教程）
- Python pyinstaller打包exe最完整教程
本站推荐
