- 积分
- 1428
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-7-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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')
|
|