学习笔记 | 用位势计算海平面气压的两种方法

用户11172986

发布于 2024-11-26 15:52:44

6700

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

学习笔记 | 用位势计算海平面气压的两种方法

温馨提示

由于可视化代码过长隐藏，可点击运行Fork查看若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

前言

在气象学的世界里，每一朵云彩背后都藏着一个故事，每一个数字背后都隐藏着大自然的秘密。今天，我们将踏上一段新的旅程，探索如何使用位势高度数据来计算海平面气压。

项目目标

理解位势高度和海平面气压的关系：揭开位势高度的神秘面纱，了解它是如何与海平面气压相互作用的。
数据准备：处理ERA5的NetCDF格式的位势高度数据，为我们的计算做好准备。
计算海平面气压：编写Python代码，将位势高度数据转化为海平面气压。
结果可视化：将计算结果可视化，用图表展现海平面气压的分布，让数据说话。

项目内容

1. 理解位势高度和海平面气压的关系

想象一下，你站在山顶上，俯瞰着下面的山谷。位势高度就像你站在山顶的高度，而海平面气压则是将你在山顶感受到的气压“拉”到海平面的高度。这个过程可以通过静力学方程来实现：

这里：

是气压的变化量，就像你从山顶下山时感受到的气压变化。
是空气的密度，就像你呼吸的空气有多重。
是重力加速度，就像地球对你的吸引力。
是高度的变化量，就像你下山的路程。

2.数据处理

现在我们需要将ERA5数据中的位势高度（geopotential height）进行单位转换

“在重力场中任一高度上，单位质量空气相对于海平面所具有的位能所表征的高度称为重力位势高度，简称位势高度。它通常以位势米为单位，并等于上述位能的1/(9.8m/s2)，比如NCEP用的位势高度就是以这个为主的。但是ECMWF直接是将单位质量空气相对于海平面所具有的位能（也就是从海平面移动到空中某一点重力所做的功）表示为位势高度，并没有除以9.8m/s2，所以单位会变为m2 s-2。用NCEP的位势高度值乘以9.8m/s**2就会变为ECMWF的值。本质上是一样的

import xarray as xr
# 打开样例数据
ds = xr.open_dataset("/home/mw/input/ERA5_Lekima4742/ERA5_Lekima.nc")
ds

<xarray.Dataset>
Dimensions:    (longitude: 401, latitude: 321, level: 7, time: 72)
Coordinates:


longitude  (longitude) float32 60.0 60.25 60.5 60.75 ... 159.5 159.8 160.0latitude   (latitude) float32 80.0 79.75 79.5 79.25 ... 0.75 0.5 0.25 0.0level      (level) int32 250 500 600 700 850 925 1000time       (time) datetime64[ns] 2019-08-08 ... 2019-08-10T23:00:00
Data variables:
z          (time, level, latitude, longitude) float32 ...
r          (time, level, latitude, longitude) float32 ...
q          (time, level, latitude, longitude) float32 ...
t          (time, level, latitude, longitude) float32 ...
u          (time, level, latitude, longitude) float32 ...
v          (time, level, latitude, longitude) float32 ...
w          (time, level, latitude, longitude) float32 ...
Attributes:
Conventions:  CF-1.6
history:      2023-05-29 04:40:37 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

# 提取位势高度数据
geop = ds['z'][0]
g = 9.81  # 重力加速度
Z = geop / g  # 位势高度
Z

<xarray.DataArray 'z' (level: 7, latitude: 321, longitude: 401)>
array([[[10222.244   , 10222.94    , 10223.81    , ..., 10281.713   ,
10281.191   , 10280.67    ],
[10228.853   , 10229.722   , 10230.417   , ..., 10291.625   ,
10291.103   , 10290.582   ],
[10235.459   , 10236.329   , 10237.198   , ..., 10301.362   ,
10300.841   , 10300.145   ],
...,
[10925.436   , 10924.913   , 10924.566   , ..., 10990.815   ,
10991.164   , 10991.338   ],
[10924.913   , 10924.566   , 10924.392   , ..., 10990.815   ,
10991.338   , 10991.859   ],
[10924.566   , 10924.392   , 10924.044   , ..., 10991.338   ,
10992.208   , 10992.902   ]],
   [[ 5474.315   ,  5474.6626  ,  5475.184   , ...,  5536.3916  ,
      5536.3916  ,  5536.3916  ],
    [ 5477.9663  ,  5478.488   ,  5478.836   , ...,  5543.695   ,
      5543.521   ,  5543.521   ],
    [ 5481.444   ,  5481.966   ,  5482.3135  , ...,  5550.4766  ,
      5550.3022  ,  5550.129   ],


...
[  774.2032  ,   773.8556  ,   773.6816  , ...,   776.1161  ,
776.2901  ,   775.7685  ],
[  774.3772  ,   774.0296  ,   773.6816  , ...,   776.46375 ,
776.1161  ,   775.7685  ],
[  774.5512  ,   774.2032  ,   773.6816  , ...,   776.46375 ,
776.1161  ,   775.9421  ]],

   [[  124.22113 ,   124.22113 ,   124.22113 , ...,   149.60857 ,
       149.95619 ,   150.30382 ],
    [  120.91734 ,   121.091354,   121.265366, ...,   152.39073 ,
       152.73836 ,   153.08597 ],
    [  117.43955 ,   117.43955 ,   117.613556, ...,   154.82487 ,
       155.3465  ,   155.69452 ],
    ...,
    [   88.574715,    87.531456,    87.18383 , ...,    93.965416,
        93.965416,    93.44379 ],
    [   88.74872 ,    87.70546 ,    87.18383 , ...,    93.791405,
        93.6174  ,    93.269775],
    [   88.92273 ,    87.70546 ,    87.18383 , ...,    93.965416,
        93.6174  ,    93.269775]]], dtype=float32)


Coordinates:

longitude (longitude) float32 60.0 60.25 60.5 60.75 ... 159.5 159.8 160.0
latitude (latitude) float32 80.0 79.75 79.5 79.25 ... 0.75 0.5 0.25 0.0
level (level) int32 250 500 600 700 850 925 1000 time datetime64[ns] 2019-08-08

计算海平面气压

# 选择1000 hPa层的位势高度数据
Zl = Z[-1, :, :]

# 假设空气密度为常数
rho = 1.255  # kg/m³
g = 9.81  # m/s²

# 计算压力变化量 dp
dp = -1 * (rho * g * Zl)  # 单位为帕斯卡 (Pa)

# 将压力变化量转换为百帕 (hPa)
dp_hPa = dp / 100

# 计算海平面气压
SLP = 1013.25 + dp_hPa  # 1000 hPa是1000 hPa层的标准气压

简单可视化

import matplotlib.pyplot as plt
import numpy as np
# 绘制海平面气压图
plt.figure(figsize=(10, 6))
plt.contourf(SLP.longitude,SLP.latitude,SLP, levels=np.linspace(980, 1040, 100), cmap='coolwarm')
plt.colorbar(label='Sea Level Pressure (hPa)')
plt.title('Sea Level Pressure from Geopotential Height')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

3. 为什么计算出的海平面气压在台风中心是高值？

位势高度的正负问题

位势高度是相对于海平面的高度。如果是正值，表示该点高于海平面；如果是负值，表示该点低于海平面。

在台风中心，由于气压非常低，位势高度通常是负值。这意味着会是一个正值（因为），从而导致计算出的海平面气压增加。

密度假设问题

我们假设空气密度为常数，但在实际情况中，空气密度会随着高度和温度的变化而变化。特别是在台风中心，空气密度可能会显著降低，这会影响的计算。

静力学方程的局限性

静力学方程是一个简化的模型，零级近似，适用于稳定的大气条件。在台风等极端天气条件下，大气动力学更为复杂，静态假设可能不再适用。

使用metpy函数验证计算结果

import metpy.calc
height = metpy.calc.geopotential_to_height(geop[-1])
p = metpy.calc.height_to_pressure_std(height)
# 绘制海平面气压图
plt.figure(figsize=(10, 6))
plt.contourf(p.longitude,p.latitude,p, levels=np.linspace(980, 1040, 100), cmap='coolwarm')
plt.colorbar(label='Sea Level Pressure (hPa)')
plt.title('Sea Level Pressure from Geopotential Height')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

np.max(p.values-SLP.values)

1.2125244

那么我们可以判定我们的计算是正确的

小结

在这次过程中我们学会了什么？

通过使用 MetPy 库，我们可以更方便地处理气象数据，特别是涉及单位转换和物理公式的计算。这样不仅可以减少出错的可能性，还可以提高代码的可读性和可维护性
大气的简易假设下物理模型面对极端天气会失去作用，那么需要进行进一步的设置与计算
位势高度与位势的概念
位势到气压的简易转换

如果你喜欢这篇文章，不妨点赞、分享给更多的朋友。关注我们的公众号，获取更多气象学和数据科学的实用技巧和知识。感谢你的支持！

参考链接 https://earthscience.stackexchange.com/questions/22998/how-to-calculate-sea-level-pressure-from-geopotential-height https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.relative_humidity_from_specific_humidity.html

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2024-11-25，如有侵权请联系 cloudcommunity@tencent.com 删除

可视化