天气学诊断与应用实习3

xiaoxiao啊

可以参考实习2的函数将*.000格式文件转化为*.txt文件
涡度平流和温度平流同样可以使用metpy函数计算，需要灵活地使用，在这里留给大家自己思考！


import numpy as np
from scipy.interpolate import interp2d
from metpy.units import units
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps

plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False


lon = np.arange(30, 160 + 2.5, 2.5)
lat = np.arange(80, 10 - 2.5, -2.5)
new_lon = np.arange(30, 160 + 0.5, 0.5)
new_lat = np.arange(10, 80 + 0.5, 0.5)


def my_interploate(*, var):
    func = interp2d(lon, lat, var, kind='cubic')
    return func(new_lon, new_lat)


def my_plot(*, level, time, var):
    u = np.loadtxt(f'u/{level}/{time}.txt')
    v = np.loadtxt(f'v/{level}/{time}.txt')
    t = np.loadtxt(f't/{level}/{time}.txt')

    u, v, t = u * units('m/s'), v * units('m/s'), t * units('degC')
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    vor = mpcalc.vorticity(u=u, v=v, dx=dx, dy=dy)
    div = mpcalc.divergence(u=u, v=v, dx=dx, dy=dy)
    new_div = my_interploate(var=div)
    new_vor = my_interploate(var=vor)

    fig = plt.figure(figsize=(6, 7), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    day = str(time)[4:6]
    ax.set_title(f'2013-05-{day}_20:00:00', loc='right')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.3, edgecolor='dimgrey')
    ax.add_feature(cfeature.LAKES.with_scale('110m'), linewidth=0.3, edgecolor='dimgrey', facecolor='none')
    ax.set_extent([90, 130, 20, 50], ccrs.PlateCarree())
    ax.set_xticks(np.arange(90, 140, 10), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(20, 60, 10), crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.tick_params(which='major', width=0.5, length=5)
    ax.spines['geo'].set_linewidth(0.5)

    if var == 'div':
        clevels = np.arange(-7, 8, 1)
        ax.set_title(f'{level}hPa Divergence', loc='left')
        ax.contour(new_lon, new_lat, new_div * 1e+5, colors='k', linewidths=0.2, levels=clevels, transform=ccrs.PlateCarree())
        cf = ax.contourf(new_lon, new_lat, new_div * 1e+5, cmap=cmaps.NCV_blu_red, levels=clevels, transform=ccrs.PlateCarree())

    if var == 'vor':
        clevels = np.arange(-12, 14, 2)
        ax.set_title(f'{level}hPa Vorticity', loc='left')
        ax.contour(new_lon, new_lat, new_vor * 1e+5, colors='k', linewidths=0.2, levels=clevels, transform=ccrs.PlateCarree())
        cf = ax.contourf(new_lon, new_lat, new_vor * 1e+5, cmap=cmaps.NCV_blu_red, levels=clevels, transform=ccrs.PlateCarree())

    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.95, pad=0.1, aspect=20, ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    plt.savefig(f'{var}_{level}_{day}.png', dpi=300, bbox_inches='tight')


my_plot(level=500, time=13052520, var='div')
my_plot(level=500, time=13052620, var='div')
my_plot(level=850, time=13052520, var='div')
my_plot(level=850, time=13052620, var='div')
my_plot(level=500, time=13052520, var='vor')
my_plot(level=500, time=13052620, var='vor')
my_plot(level=850, time=13052520, var='vor')
my_plot(level=850, time=13052620, var='vor')

雪域小丹

这个是NCL的图片吧？不太像python画的，楼主上传错了吧？

来自苹果APP客户端

xiaoxiao啊

雪域小丹 发表于 2023-5-8 03:31
这个是NCL的图片吧？不太像python画的，楼主上传错了吧？

python画的，NCL画的比我好看

雪域小丹

xiaoxiao啊 发表于 2023-6-2 16:02
python画的，NCL画的比我好看

平滑用的是什么呀？效果怎么这么好？