请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 5985|回复: 2

[经验总结] 天气学诊断与应用实习2

[复制链接]

新浪微博达人勋

发表于 2022-7-2 21:57:44 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 xiaoxiao啊 于 2022-7-2 21:59 编辑

首先,可以先将老师所给的*****.000格式的micaps数据处理为便于numpy读取的*****.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


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


def micaps_to_txt(*, var, level, time, skip_rows):
    file_micaps = f'{var}/{level}/{time}.000'
    file_txt = f'{var}/{level}/{time}.txt'
    with open(file_micaps, mode='r', encoding='gbk') as f1:
        with open(file_txt, mode='w') as f2:
            for i in range(skip_rows):
                f1.readline()  # 过滤无用行
            for i in range(29):
                string = ''
                for j in range(6):
                    string += f1.readline().split('\n')[0]  # 去除换行符
                f2.write(string)
                f2.write('\n')


# 500hPa
micaps_to_txt(var='t', level=500, time=13052520, skip_rows=4)
micaps_to_txt(var='t', level=500, time=13052620, skip_rows=4)
micaps_to_txt(var='t-td', level=500, time=13052520, skip_rows=4)
micaps_to_txt(var='t-td', level=500, time=13052620, skip_rows=4)
# 850hPa
micaps_to_txt(var='t', level=850, time=13052520, skip_rows=4)
micaps_to_txt(var='t', level=850, time=13052620, skip_rows=4)
micaps_to_txt(var='t-td', level=850, time=13052520, skip_rows=4)
micaps_to_txt(var='t-td', level=850, time=13052620, skip_rows=4)


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_plot(*, level, time, var):
    t = np.loadtxt(f't/{level}/{time}.txt')
    delta_t = np.loadtxt(f't-td/{level}/{time}.txt')
    td = t - delta_t
    t, td = t * units('degC'), td * units('degC')
    p = 500 * units('hPa')
    theta = mpcalc.potential_temperature(p, t)
    theta_e = mpcalc.equivalent_potential_temperature(p, t, td)
   
    func = interp2d(lon, lat, theta, kind='cubic')
    new_theta = func(new_lon, new_lat)
    func = interp2d(lon, lat, theta_e, kind='cubic')
    new_theta_e = func(new_lon, new_lat)
   
    fig = plt.figure(figsize=(6, 6), 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([30, 150, 10, 70], ccrs.PlateCarree())
    ax.set_xticks(np.arange(30, 160, 30), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(10, 80, 20), 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)
    ax.set_aspect(1.3)
   
    if var == 'theta':
        if level == 500:
            clevels=np.arange(290, 344, 4)
        elif level == 850:
            clevels = np.arange(315, 380, 5)
        ax.set_title(f'{level}hPa Potential Temperature', loc='left')
        cs = ax.contour(new_lon, new_lat, new_theta, colors='k', levels=clevels, linewidths=0.2, transform=ccrs.PlateCarree())
        cf = ax.contourf(new_lon, new_lat, new_theta, cmap='jet', levels=clevels, transform=ccrs.PlateCarree())
   
    elif var == 'theta_e':
        if level == 500:
            clevels=np.arange(290, 360, 5)
        elif level == 850:
            clevels = np.arange(300, 520, 20)
        ax.set_title(f'{level}hPa Equivalent Potential Temperature', loc='left')
        cs = ax.contour(new_lon, new_lat, new_theta_e, colors='k', levels=clevels, linewidths=0.3, transform=ccrs.PlateCarree())
        cf = ax.contourf(new_lon, new_lat, new_theta_e, cmap='jet', levels=clevels, transform=ccrs.PlateCarree())
   
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.9, pad=0.1)
    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='theta')
my_plot(level=500, time=13052620, var='theta')
my_plot(level=850, time=13052520, var='theta')
my_plot(level=850, time=13052620, var='theta')
my_plot(level=500, time=13052520, var='theta_e')
my_plot(level=500, time=13052620, var='theta_e')
my_plot(level=850, time=13052520, var='theta_e')
my_plot(level=850, time=13052620, var='theta_e')

theta_500_25.png
theta_500_26.png
theta_850_25.png
theta_850_26.png
theta_e_500_25.png
theta_e_500_26.png
theta_e_850_25.png
theta_e_850_26.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-7-5 23:26:39 | 显示全部楼层
请问有数据格式可以看一下吗
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-7-6 18:04:55 | 显示全部楼层
子非鱼。 发表于 2022-7-5 23:26
请问有数据格式可以看一下吗

micaps文件的格式吗?我们老师给的是第11类文件
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表