经典相当位温的Python 实现
前言
本项目旨在通过Python实现经典相当位温的计算方法，帮助大家更好理解位温概念。无论你是从事气象科研，还是从事天气预报，掌握相当位温，能帮助你更好地了解大气状态。
内容包括：相当位温的基本概念，计算方法，Python代码示例与简单可视化
相当位温（Equivalent Potential Temperature）是描述大气状态的一个重要指标。它是指将某一气块抬升到凝结高度，并使其水汽凝结释放所有潜热后得到的位温。换句话说，相当位温表示了气块在绝热抬升至相同压强下的稳定状态下的温度。
为什么相当位温如此重要呢？这是因为它具有以下几个特点：
包含了水汽的影响：相当位温考虑了水汽在临界抬升过程中释放的潜热，因此能够较准确地描述湿空气的状态。
反映了稳定性：相当位温是一个稳定性指标，稳定的大气层中相当位温变化较小，而不稳定的大气层中相当位温随高度增加而减小。
描述了气块的来源：相当位温还可以用来区分气块的不同来源，比如热带或极地地区。
通过计算和分析相当位温，我们可以更好地理解大气的垂直结构、判断不同气团的性质和运动趋势，对天气的形成和演变提供重要参考。因此，在气象科研、天气预报和气候研究等领域，相当位温是一项必不可少的指标。
在接下来的项目中，我们将详细介绍相当位温的计算方法和应用，帮助您更好地理解和应用这一重要概念。
Holton（1972）针对有凝结过程参与的饱和大气，提出的相当位温 公式为 ��=�exp(��0�����)
其中，θ 为位温，Lv0＝2.555×106 J kg-1 代表凝结潜热，qs 为水汽饱和比湿，cp＝1004 J kg-1 K-1 代表干空气的定压比热
参考文献：周括, 冉令坤, 齐彦斌, 等. 2020. 包含冻结过程的广义位温及位涡特征分析 [J]. 大气科学, 44(4): 816−834. ZHOU Kuo, RAN Lingkun, QI Yanbin, et al. 2020. Characteristic Analysis of Generalized Potential Temperature and Potential Vorticity during Freezing [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 44(4): 816−834. doi:10.3878/j.issn.1006-9895.1908.19154
小知识
html中上标的输入格式是数字< sup > 数字 < /sup > latex中分号是frac{{a}}{{b}} 下标格式为 a_b ,有复数字母时需要用{}括起来
1. 定义函数
import numpy as np
def calculate_theta_e(theta, qs, T):
Lv0 = 2.555 * 1e3
Cp = 1004
exponent = (Lv0 * qs) / (Cp * T)
theta_e = theta * np.exp(exponent)
return theta_e
# 示例调用
theta = 300 # 位温
qs = 20 # 水汽饱和比湿
T = 300 # 温度result = calculate_theta_e(theta, qs, T)
print(result)355.46869600864052. 饱和比湿计算
那么qs（饱和比湿）哪里来呢，检索一下，懂了，从es（饱和水汽压）计算
那么es从哪里来，答案是metpy函数metpy.calc import saturation_vapor_pressure
当然，metpy并没有直接计算饱和比湿的函数，倒是有饱和混合比
from metpy.calc import saturation_vapor_pressure
from metpy.units import units
p = 500 * units.hPa
es = saturation_vapor_pressure(5 * units.degC).to('hPa')
qs = (622 *es)/(p - 0.378*es)/1000
qs0.010921512911518452 dimensionless
3. 实际应用 : WRF后处理提取相关变量计算相当位温
还是从老伙计wrfout中提取需要的变量：位温 温度 气压
In [22]:
from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,ALL_TIMES
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants
import os
# 定义 WRF 文件夹路径和文件名前缀
wrfout_path = "/home/mw/input/wrfout3385/"
filename_prefix = "wrfout_d02_"
wrf_files = sorted([os.path.join(wrfout_path, f) for f in os.listdir(wrfout_path) if f.startswith(filename_prefix)])
wrf_list = [Dataset(f) for f in wrf_files]In [23]:
theta_wrf = getvar(wrf_list, 'theta',timeidx=ALL_TIMES, method='cat')
tk = getvar(wrf_list, 'tk', timeidx=ALL_TIMES, method='cat')
p_wrf = getvar(wrf_list, 'pressure', timeidx=ALL_TIMES, method='cat')
tc = getvar(wrf_list, 'tc', timeidx=ALL_TIMES, method='cat')
td = getvar(wrf_list, 'td', timeidx=ALL_TIMES, method='cat')In [24]:
es_wrf = saturation_vapor_pressure(to_np(tc) * units.degC).to('hPa')
qs_wrf = (622 *es_wrf)/(p_wrf* units.hPa - 0.378*es_wrf)/1000
result_wrf = calculate_theta_e(theta_wrf,qs_wrf,tk)
result_wrf.shape(9, 49, 90, 90)In [25]:
result_wrf[5,20,:,:].plot(cmap='bwr')
<matplotlib.collections.QuadMesh at 0x7f946ebb4b10>
4. metpy的相当位温
公式为
from here metpy
In [26]:
from metpy.calc import equivalent_potential_temperature
from metpy.units import units
result2 = equivalent_potential_temperature(p_wrf*units.hPa, tc *units.degC, td*units.degC)
result2[5,20,:,:].plot(cmap='bwr')<matplotlib.collections.QuadMesh at 0x7f946eb0fad0>
4. 小结
- 本来想细化绘图，犯懒了直接plot
- 是不是觉得这俩位温差得也太大了。只能说同为位温，亦有差距。当然也可能我算错了，欢迎指正。
- 位温概念还有许多种，有兴趣可自行尝试比较。