数据处理 | python计算经典相当位温2.0

文摘   2024-12-11 18:11   北京  

相当位温

镜像:气象分析3.7
数据:wrfout文件

🎓 作者:酷炫用户名

📢 版权声明:公益性质转载需联系作者本人获取授权。转载本文时,请务必文字注明“来自:和鲸社区:酷炫用户名”,并附带本项目超链接。

前言

本项目旨在通过Python实现经典相当位温的计算方法,帮助大家更好理解位温概念。无论你是从事气象科研,还是从事天气预报,掌握相当位温,能帮助你更好地了解大气状态。

内容包括:相当位温的基本概念,计算方法,Python代码示例与简单可视化

相当位温(Equivalent Potential Temperature)是描述大气状态的一个重要指标。它是指将某一气块抬升到凝结高度,并使其水汽凝结释放所有潜热后得到的位温。换句话说,相当位温表示了气块在绝热抬升至相同压强下的稳定状态下的温度。

为什么相当位温如此重要呢?这是因为它具有以下几个特点:

包含了水汽的影响:相当位温考虑了水汽在临界抬升过程中释放的潜热,因此能够较准确地描述湿空气的状态。

反映了稳定性:相当位温是一个稳定性指标,稳定的大气层中相当位温变化较小,而不稳定的大气层中相当位温随高度增加而减小。

描述了气块的来源:相当位温还可以用来区分气块的不同来源,比如热带或极地地区。

通过计算和分析相当位温,我们可以更好地理解大气的垂直结构、判断不同气团的性质和运动趋势,对天气的形成和演变提供重要参考。因此,在气象科研、天气预报和气候研究等领域,相当位温是一项必不可少的指标。

在接下来的项目中,我们将详细介绍相当位温的计算方法和应用,帮助您更好地理解和应用这一重要概念。

Holton(1972)针对有凝结过程参与的饱和大气,提出的相当位温
公式为

其中,θ 为位温,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 * 1e6
    Cp = 1004
    exponent = (Lv0 * qs) / (Cp * T)
    theta_e = theta * np.exp(exponent)
    return theta_e

# 示例调用
theta = 300  # 位温
qs = 0.01  # 水汽饱和比湿
T = 300     # 温度

result = calculate_theta_e(theta, qs, T)
print(result)
326.5587371401846

2. 饱和比湿计算

那么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
qs

0.010921512911518452 dimensionless

此处的单位为kg/kg

3. 实际应用 : WRF后处理提取相关变量计算相当位温

还是从老伙计wrfout中提取需要的变量:位温 温度 气压

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/typhoon9537"
filename_prefix = "wrfout_d01_"

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]
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')
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[3,20,:,:].plot(cmap='bwr')
<matplotlib.collections.QuadMesh at 0x7f961dff7990>

4. metpy的相当位温

公式为

from  heremetpy

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[3,20,:,:].plot(cmap='bwr')
<matplotlib.collections.QuadMesh at 0x7f95b6ee5790>

4. 小结

  1. 本来想细化绘图,犯懒了直接plot
  2. 是不是觉得这俩位温差得也太大了。只能说同为位温,亦有差距。当然也可能我算错了,欢迎指正。
  3. 位温概念还有许多种,有兴趣可自行尝试比较。

气python风雨
主要发一些涉及大气科学的Python文章与个人学习备忘录
 最新文章