python的气象应用

修—养

已于 2024-06-14 20:25:18 修改

阅读量908

点赞数 14

文章标签： python 开发语言 python3.11

于 2024-06-14 20:20:08 首次发布

本文链接：https://blog.csdn.net/sudo_xiuyang/article/details/139689422

版权

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from wrf import getvar, interplevel, CoordPair, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords, geo_bounds, to_np
import netCDF4 as nc

print('import OK !')

####
ncfile = nc.Dataset('wrfout_d01_2023-08-01_18_00_00.wrfout')
print('ncfile OK !')
# print(ncfile)

ntimes = ncfile.dimensions['Time'].size
nlons = ncfile.dimensions['west_east'].size
nlats = ncfile.dimensions['south_north'].size
nlevs = ncfile.dimensions['bottom_top'].size
XLONG = ncfile.variables['XLONG'][0, :, :]
XLAT = ncfile.variables['XLAT'][0, :, :]
leftlon = np.amin(XLONG)
leftlon = round(leftlon - 0.5)
rightlon = np.amax(XLONG)
rightlon = round(rightlon + 0.5)
lowerlat = np.amin(XLAT)
lowerlat = round(lowerlat - 0.5)
upperlat = np.amax(XLAT)
upperlat = round(upperlat + 0.5)
print("ntimes=" + str(ntimes))
print("nlevs=" + str(nlevs))
print("nlats=" + str(nlats))
print("nlons=" + str(nlons))
print("leftlon =" + str(leftlon))
print("rightlon=" + str(rightlon))
print("lowerlat=" + str(lowerlat))
print("upperlat=" + str(upperlat))
Region = [leftlon, rightlon, lowerlat, upperlat] # Ҫ��Ƶķ�Χ lon1,lon2,lat1,lat2
proj = ccrs.PlateCarree()
proj_data = ccrs.PlateCarree()

for i in range(0, ntimes, 6):
print("Time" + str(i))

times = getvar(ncfile, "Times", i)
lats = getvar(ncfile, "XLAT", i)
lons = getvar(ncfile, "XLONG", i)

# ��θ߶�
hgt = getvar(ncfile, "HGT", i)
# ��ƽ��ѹ��
slp = getvar(ncfile, "slp", i)
# ��2m�� λת��Ϊ��϶�
t0 = getvar(ncfile, "T2", i)
t2 = t0 - 273.15
# �ۼƽ�ˮ ��ˮ+��ƽ�ˮ
rainc = getvar(ncfile, "RAINC", i)
rainnc = getvar(ncfile, "RAINNC", i)
rain = rainc + rainnc
# ��߶��
p = getvar(ncfile, "pressure", i)
# print(p.shape)
# �ر��ѹ hPa
Ps = p[0, :, :]
# print(Ps.shape)
# print(Ps)
z = getvar(ncfile, "z", i)
ua = getvar(ncfile, "ua", i)
va = getvar(ncfile, "va", i)
wspd = getvar(ncfile, "wspd_wdir", i)[0, :]
# wspd = getvar(ncfile, "wspd_wdir",i)
# print(wspd.shape)
# print(wspd)
rh = getvar(ncfile, "rh", i)
qv = getvar(ncfile, "QVAPOR", i)
# print(qv)
# Interpolate geopotential height, u, and v winds to 500 hPa
# 500hpaλ�Ƹ߶ȳ�
ht_500mb = interplevel(z, p, 500.)
ht_500mb = ht_500mb / 10.
# print(np.min(ht_500mb))
# print(np.max(ht_500mb))
# lon = ht_500mb.XLONG
# lat = ht_500mb.XLAT
# print(ht_500mb.shape)
# 500hpa��
ht_500 = interplevel(z, p, 500)
u_500 = interplevel(ua, p, 500)
v_500 = interplevel(va, p, 500)
wspd_500 = interplevel(wspd, p, 500)
# 7 00hpa��+��ʪ��
rh_700 = interplevel(rh, p, 700)
u_700 = interplevel(ua, p, 700)
v_700 = interplevel(va, p, 700)
wspd_700 = interplevel(wspd, p, 700)
# 850hpa��+ˮ��ϱ�
qv_850 = interplevel(qv, p, 850)
u_850 = interplevel(ua, p, 850)
v_850 = interplevel(va, p, 850)
wspd_850 = interplevel(wspd, p, 850)
qv_850 = interplevel(rh, p, 850)

# ��HGT ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
levels = np.arange(300, 2400, 300) ##
ac = ax1.contourf(lons, lats, hgt, cmap='RdYlBu_r', extend='both', levels=levels, transform=proj_data)
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout HGT time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() # bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/hgt_time%d.png' % (i + 1))
plt.close()

# ��Ƶر��ѹ ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
levels = np.arange(700, 1050, 50)
ac = ax1.contourf(lons, lats, Ps, cmap='RdYlBu_r', extend='both', levels=levels, transform=proj_data)
ax1.set_xticks(np.arange(leftlon, rightlon, 5))
ax1.set_yticks(np.arange(lowerlat, upperlat, 5))
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout Surface Pressure time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() # bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/Ps_time%d.png' % (i + 1))
plt.close()

# ��SLP ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
levels = np.arange(996, 1024, 4)
ac = ax1.contourf(lons, lats, slp, cmap='RdYlBu_r', levels=levels, extend='both', transform=proj_data)
ax1.set_xticks(np.arange(leftlon, rightlon, 5))
ax1.set_yticks(np.arange(lowerlat, upperlat, 5))
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout slp time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries()
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/slp_time%d.png' % (i + 1))
plt.close()

# ��T2 ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, t2, levels=np.arange(12, 36, 3), cmap='Reds', extend='both') ##
ax1.set_xticks(np.arange(leftlon, rightlon, 5))
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout T2 time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/T2_time%d.png' % (i + 1))
plt.close()

# ��rain ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, rain, levels=np.arange(1.5, 10.5, 1.5), cmap='RdYlBu_r', extend='both') ##
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout rain time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/rain_time%d.png' % (i + 1))
plt.close()

#��500hpaλ�Ƹ߶ȳ� ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, ht_500mb, levels=np.arange(540, 600, 10), cmap='RdYlBu_r', extend='both') ##
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout ht_500mb time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/ht500mb_time%d.png' % (i + 1))
plt.close()

#��500hpa�糡 ��Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, wspd_500, levels=np.arange(0, 30, 3), cmap='RdYlBu_r', extend='both') ##
ax1.quiver(to_np(lons[::3, ::3]), to_np(lats[::3, ::3]), to_np(u_500[::3, ::3]), to_np(v_500[::3, ::3]),
to_np(wspd_500[::3, ::3]), transform=ccrs.PlateCarree(), scale=500, cmap='jet')
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout 500hpa wind' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/wind500_time%d.png' % (i + 1))
plt.close()

#��700�糡+��ʪ�� Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, rh_700, levels=np.arange(15, 105, 15), cmap='Greens', extend='both') ##
#ax1.quiver(lons[::3],lats[::3],u_700[::3,::3],v_700[::3,::3],wspd_700[::3,::3],transform=ccrs.PlateCarree(),scale=500,cmap=plt.cm.jet)
ax1.barbs(lons[::3, ::3], lats[::3, ::3], u_700[::3, ::3], v_700[::3, ::3], wspd_700[::3, ::3], length=4,
transform=ccrs.PlateCarree())
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout 700hpa wind+RH time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/wind&RH700_time%d.png' % (i + 1))
plt.close()
#help(plt.barbs)

#��850hpa�糡+��ʪ�� Ҫ�޸�levels�Ȳ��
fig = plt.figure(figsize=(12, 6), dpi=150)
ax1 = plt.subplot(121, projection=proj)
ax1.set_extent(Region, crs=proj_data)
ac = ax1.contourf(lons, lats, qv_850, levels=np.arange(15, 105, 15), cmap='Greens', extend='both') ##
#ax1.quiver(lons[::3],lats[::3],u_700[::3,::3],v_700[::3,::3],wspd_700[::3,::3],transform=ccrs.PlateCarree(),scale=500,cmap=plt.cm.jet)
ax1.barbs(lons[::3, ::3], lats[::3, ::3], u_850[::3, ::3], v_850[::3, ::3], wspd_850[::3, ::3], length=4,
transform=ccrs.PlateCarree())
ax1.set_xticks(np.arange(leftlon, rightlon, 5)) ##
ax1.set_yticks(np.arange(lowerlat, upperlat, 5)) ##
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_title('WRFout 850hpa wind+RH time' + str(i + 1))
ax1.add_feature(cfeat.COASTLINE, linewidth=0.6, zorder=1)
SHP = Reader('./shp/bou2_4p.shp').geometries() #bou2_4p.shp ʡ��
ax1.add_geometries(SHP, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.5)
c1bar = fig.colorbar(ac, ax=ax1, fraction=0.02, pad=0.025, orientation='vertical')
plt.show()
plt.savefig('wrfplots/wind&RH850_time%d.png' % (i + 1))
plt.close()