metpy绘制锋生与冷锋

用户11172986

发布于 2024-06-20 17:37:39

810

发布于 2024-06-20 17:37:39

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

以下全文代码和数据均已发布至和鲸社区，阅读原文前往，可一键fork跑通

https://www.heywhale.com/mw/project/65485a22d74b63fed5f03f49

本项目旨在利用Python中的MetPy库来进行气象数据分析，特别是关于锋生函数和冷锋的研究与绘制。

镜像：气象分析3.9

锋生函数公式如下

这是11月7日的中央气象台的地面天气图

导入库

import xarray as xr
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from cartopy.io.shapereader import BasicReader

读取数据

data = xr.open_dataset('/home/mw/input/1107125177/2023110720.nc')
provinces = BasicReader('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp')

计算位温与锋生

计算位温与锋生
# 从ERA5数据集中读取500hPa高度和温度数据
lon=data.longitude
lat=data.latitude
#heights_500hPa = data['p'].sel(level=500,time=data.time[0])
temperatures_500hPa = data['t'].sel(level=1000,time=data.time[0])
u=  data['u'].sel(level=1000 ,time=data.time[0])
v=  data['v'].sel(level=1000,time=data.time[0] )
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# 计算位温
p =850. *units.mbars
potential_temperature = mpcalc.potential_temperature(p, temperatures_500hPa* units.kelvin)


# 计算500hPa上的锋生函数
frontogenesis = mpcalc.frontogenesis(potential_temperature, u, v, dx=dx, dy=dy)

‍

绘图部分

# 绘制位温等值线图
convert_to_per_100km_3h = 1000*100*3600*3
fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())
cs = ax.contour(data.longitude, data.latitude, potential_temperature, colors='k',levels=np.arange(270, 310, 2))
cf = ax.contourf(data.longitude, data.latitude, frontogenesis*convert_to_per_100km_3h, cmap='coolwarm',levels=np.arange(-8,12, 0.5))
plt.colorbar(cf, ax=ax, label='K/100km/3h')
ax.set_title('1000 hPa Frontogenesis')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.add_geometries(provinces.geometries(), linewidth=.5, edgecolor='black', crs=ccrs.PlateCarree(),
                  facecolor='none')


ax.set_extent( [100, 130, 20,60])
ax.set_xticks(range(100, 130, 5))  # 设置经度刻度
ax.set_yticks(range(20,60, 5))  # 设置纬度刻度
plt.show()

绘图过程较慢，可缩小绘制范围或对程序进行优化

冷锋绘制

from metpy.plots import ColdFront, WarmFront
对曲线需要多次调整，有更好的曲线绘制方法可评论区讨论


low_lon, low_lat = 120, 50  # 起点经纬度
h_lon,h_lat = 98, 49
# 生成冷锋曲线坐标
cold_lat = np.linspace(low_lat , low_lat - 10, 100)
slope_multiplier = 0.1 # 坡度增大因子
cold_lon = (low_lon + 0.25) - slope_multiplier * (cold_lat - (low_lat - 0.5))**2
convert_to_per_100km_3h = 1000*100*3600*3
# 创建曲线路径
vertices = np.column_stack((cold_lon, cold_lat))
path = mpath.Path(vertices)


# 创建路径补丁
patch = mpatches.PathPatch(path, fc="none")
##建立图像
fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())


ax.text(low_lon, low_lat, 'L', color='red', size=30,
        horizontalalignment='center', verticalalignment='center')
ax.text(h_lon, h_lat, 'H', color='blue', size=30,
        horizontalalignment='center', verticalalignment='center')
ax.add_patch(patch)
ax.plot(cold_lon, cold_lat , 'blue', path_effects=[ColdFront(size=8, spacing=1.5)])
cs = ax.contour(data.longitude, data.latitude, potential_temperature, colors='k',levels=np.arange(270, 310, 5))
mesh = ax.pcolormesh(data.longitude, data.latitude, frontogenesis*convert_to_per_100km_3h, cmap='coolwarm', vmin=-8, vmax=8,
                     transform=ccrs.PlateCarree())
# cf = ax.contourf(data.longitude, data.latitude, frontogenesis*convert_to_per_100km_3h, cmap='coolwarm',levels=np.arange(-8,12, 0.5))
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.add_geometries(provinces.geometries(), linewidth=.5, edgecolor='black', crs=ccrs.PlateCarree(),
                  facecolor='none')
ax.set_extent( [90, 130, 20,60])
ax.set_xticks(range(100, 130, 5))  # 设置经度刻度
ax.set_yticks(range(20,60, 5))  # 设置纬度刻度  
plt.colorbar(mesh)  # 添加颜色条              
plt.show()