-
Notifications
You must be signed in to change notification settings - Fork 0
/
mwe1_metpy.py
93 lines (64 loc) · 3.29 KB
/
mwe1_metpy.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import os.path
import xarray as xr
import numpy as np
import boto3
import metpy.calc as mpcalc
from botocore import UNSIGNED
from botocore.config import Config
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
#from matplotlib.colormaps import get_cmap
import matplotlib.ticker as mticker
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import sys
import time
if (not os.path.isfile('gfs.t12z.pgrb2.0p25.f000')):
client = boto3.client('s3', config=Config(signature_version=UNSIGNED))
client.download_file('noaa-gfs-bdp-pds', 'gfs.20230809/12/atmos/gfs.t12z.pgrb2.0p25.f000', 'gfs.t12z.pgrb2.0p25.f000')
u850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'u', 'level': 850}})
u = u850.u
print(u.shape)
v850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib', backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'v', 'level': 850}})
v = v850.v
#print(v.shape)
# Compute the 850 hPa relative vorticity.
vort850 = mpcalc.vorticity(u, v)
# Compute the 850 hPa divergence.
div850 = mpcalc.divergence(u, v)
vortmask = mpcalc.bounding_box_mask(vort850,5.,13.5,191.,202.)
divmask = mpcalc.bounding_box_mask(div850,5.,13.5,191.,202.)
i_bb_indices = mpcalc.find_bounding_box_indices(vortmask,5.,13.5,191.,202.)
o_bb_indices = mpcalc.find_bounding_box_indices(vortmask,0.,30,180.,220)
dx, dy = mpcalc.lat_lon_grid_deltas(vortmask.longitude, vortmask.latitude)
upsi,vpsi = mpcalc.rotational_wind_from_inversion(vortmask,dx,dy,o_bb_indices,i_bb_indices)
# Create a set of axes for the figure and set
# its map projection to that of the input data.
ax = plt.axes(projection=crs.PlateCarree())
# Add country borders and coastlines.
countries = NaturalEarthFeature(category="cultural", scale="50m",
facecolor="none",
name="admin_0_countries")
ax.add_feature(countries, linewidth=.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)
# Compute the magnitude of the non-divergent component of the 850 hPa wind.
nd_spd = np.sqrt(upsi**2 + vpsi**2)
# Plot this using xarray's plot functionality.
plot = nd_spd.plot(levels=np.arange(0., 13., 1.), cmap=get_cmap('YlGnBu'), transform=crs.PlateCarree(), cbar_kwargs={'label':'non-divergent wind ($m s^{-1}$)', 'shrink': 0.98})
# Set the map's extent to match that over which we computed the non-divergent wind.
ax.set_extent([-180.,-140.,0.,30.],crs=crs.PlateCarree())
# Add latitude/longitude gridlines.
gridlines = ax.gridlines(color="grey", linestyle="dotted", draw_labels=True)
gridlines.top_labels = False
gridlines.right_labels = False
gridlines.xlocator = mticker.FixedLocator(np.arange(-180.,139.,5.))
gridlines.ylocator = mticker.FixedLocator(np.arange(0.,31.,5.))
gridlines.xlabel_style = {'size':12, 'color':'black'}
gridlines.ylabel_style = {'size':12, 'color':'black'}
gridlines.xformatter = LONGITUDE_FORMATTER
gridlines.yformatter = LATITUDE_FORMATTER
# Add a plot title, then show the image.
plt.title("GFS 0-h 850 hPa non-divergent wind magnitude ($m s^{-1}$) due to Dora at 1200 UTC 9 August 2023")
plt.savefig('vectorized_version')
plt.show()