import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import metpy
import metpy.calc as mpcalc
from metpy.plots import ctables
from metpy.cbook import get_test_data
from metpy.units import units
import os
import scipy.integrate as integrate
import datetime as dt
import glob
import json
from datetime import datetime
from datetime import timedelta
from matplotlib.colors import Normalize
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, BoundaryNorm
import scipy
file_dir = '/data/iconeu/'
# Changing the directory
os.chdir(file_dir)
data = xr.open_dataset('icon_uv10.nc').sel(longitude=slice(-5,15,1),latitude=slice(40,55,1))
data_ml = xr.open_dataset('icon_ml.nc').sel(longitude=slice(-5,15,1),latitude=slice(40,55,1))
data_sfc = xr.open_dataset('icon_pmsl.nc').sel(longitude=slice(-5,15,1),latitude=slice(40,55,1))
data_sfc2 = xr.open_dataset('icon_cape.nc').sel(longitude=slice(-5,15,1),latitude=slice(40,55,1))
# View a summary of the Dataset
print(data_ml)
print(data)
<xarray.Dataset>
Dimensions: (longitude: 321, latitude: 241, level: 40, time: 50)
Coordinates:
* longitude (longitude) float32 -5.0 -4.938 -4.875 ... 14.88 14.94 15.0
* latitude (latitude) float32 40.0 40.06 40.12 40.19 ... 54.88 54.94 55.0
* level (level) int32 30 31 32 33 34 35 36 37 ... 62 63 64 65 66 67 68 69
* time (time) datetime64[ns] 2023-07-13 ... 2023-07-15T01:00:00
Data variables:
u (time, level, latitude, longitude) float32 ...
v (time, level, latitude, longitude) float32 ...
q (time, level, latitude, longitude) float32 ...
Attributes:
Conventions: CF-1.6
history: 2023-07-13 04:44:51 GMT by grib_to_netcdf-2.6.0: grib_to_ne...
<xarray.Dataset>
Dimensions: (longitude: 321, latitude: 241, time: 50)
Coordinates:
* longitude (longitude) float32 -5.0 -4.938 -4.875 ... 14.88 14.94 15.0
* latitude (latitude) float32 40.0 40.06 40.12 40.19 ... 54.88 54.94 55.0
* time (time) datetime64[ns] 2023-07-13 ... 2023-07-15T01:00:00
Data variables:
u10 (time, latitude, longitude) float32 ...
v10 (time, latitude, longitude) float32 ...
Attributes:
Conventions: CF-1.6
history: 2023-07-13 04:47:32 GMT by grib_to_netcdf-2.6.0: grib_to_ne...
# To parse the full dataset, we can call parse_cf without an argument, and assign the returned Dataset.
data = data.metpy.parse_cf()
data_ml = data_ml.metpy.parse_cf()
data_sfc = data_sfc.metpy.parse_cf()
data_sfc2 = data_sfc2.metpy.parse_cf()
x, y = data_ml['u'].metpy.coordinates('x', 'y')
time = data_ml['u'].metpy.time
timeinit = time[0]
timeinit = datetime.utcfromtimestamp(timeinit.item()/1e9)-timedelta(hours=0)
print(timeinit)
mslp = data_sfc['prmsl']
mslp = mslp/100
usfc = data['u10']
vsfc = data['v10']
cape = data_sfc2['CAPE_ML']
u250m = data_ml['u'].metpy.loc[{'level': 57}]
v250m = data_ml['v'].metpy.loc[{'level': 57}]
u500m = data_ml['u'].metpy.loc[{'level': 55}]
v500m = data_ml['v'].metpy.loc[{'level': 55}]
u750m = data_ml['u'].metpy.loc[{'level': 53}]
v750m = data_ml['v'].metpy.loc[{'level': 53}]
u1km = data_ml['u'].metpy.loc[{'level': 51}]
v1km = data_ml['v'].metpy.loc[{'level': 51}]
u1_5km = data_ml['u'].metpy.loc[{'level': 49}]
v1_5km = data_ml['v'].metpy.loc[{'level': 49}]
u2km = data_ml['u'].metpy.loc[{'level': 47}]
v2km = data_ml['v'].metpy.loc[{'level': 47}]
u2_5km = data_ml['u'].metpy.loc[{'level': 45}]
v2_5km = data_ml['v'].metpy.loc[{'level': 45}]
u3km = data_ml['u'].metpy.loc[{'level': 43}]
v3km = data_ml['v'].metpy.loc[{'level': 43}]
u3_5km = data_ml['u'].metpy.loc[{'level': 41}]
v3_5km = data_ml['v'].metpy.loc[{'level': 41}]
u4km = data_ml['u'].metpy.loc[{'level': 39}]
v4km = data_ml['v'].metpy.loc[{'level': 39}]
u6km = data_ml['u'].metpy.loc[{'level': 34}]
v6km = data_ml['v'].metpy.loc[{'level': 34}]
q_ll = data_ml['q'].where(data_ml['q'].level>55, drop=True)
qmean = np.mean(q_ll, axis=1)
ulls = u1km - usfc
vlls = v1km - vsfc
umls = u3km - usfc
vmls = v3km - vsfc
udls = u6km - usfc
vdls = v6km - vsfc
lls = np.sqrt(ulls**2 + vlls**2)
mls = np.sqrt(umls**2 + vmls**2)
dls = np.sqrt(udls**2 + vdls**2)
u_motion = data_ml['u'].where(data_ml['u'].level>34, drop=True)
v_motion = data_ml['v'].where(data_ml['v'].level>34, drop=True)
umean = np.mean(u_motion, axis=1)
vmean = np.mean(v_motion, axis=1)
smotion = np.sqrt(umean**2 + vmean**2)
ubunkers = umean+(7.5/dls)*vdls
vbunkers = vmean-(7.5/dls)*udls
sr1=np.sqrt((usfc-ubunkers)**2+(vsfc-vbunkers)**2)
sr2=np.sqrt((u500m-ubunkers)**2+(v500m-vbunkers)**2)
sr3=np.sqrt((u1km-ubunkers)**2+(v1km-vbunkers)**2)
sr4=np.sqrt((u1_5km-ubunkers)**2+(v1_5km-vbunkers)**2)
sr5=np.sqrt((u2km-ubunkers)**2+(v2km-vbunkers)**2)
srmean_ll = (sr1+sr2+sr3+sr4+sr5)/5
sr6=np.sqrt((u2km-ubunkers)**2+(v2km-vbunkers)**2)
sr7=np.sqrt((u2_5km-ubunkers)**2+(v2_5km-vbunkers)**2)
sr8=np.sqrt((u3km-ubunkers)**2+(v3km-vbunkers)**2)
sr9=np.sqrt((u3_5km-ubunkers)**2+(v3_5km-vbunkers)**2)
sr10=np.sqrt((u4km-ubunkers)**2+(v4km-vbunkers)**2)
srmean_ul = (sr6+sr7+sr8+sr9+sr10)/5
sr1m=np.sqrt((usfc-umean)**2+(vsfc-vmean)**2)
sr2m=np.sqrt((u500m-umean)**2+(v500m-vmean)**2)
sr3m=np.sqrt((u1km-umean)**2+(v1km-vmean)**2)
sr4m=np.sqrt((u1_5km-umean)**2+(v1_5km-vmean)**2)
sr5m=np.sqrt((u2km-umean)**2+(v2km-vmean)**2)
srmean_llm = (sr1m+sr2m+sr3m+sr4m+sr5m)/5
sr6m=np.sqrt((u2km-umean)**2+(v2km-vmean)**2)
sr7m=np.sqrt((u2_5km-umean)**2+(v2_5km-vmean)**2)
sr8m=np.sqrt((u3km-umean)**2+(v3km-vmean)**2)
sr9m=np.sqrt((u3_5km-umean)**2+(v3_5km-vmean)**2)
sr10m=np.sqrt((u4km-umean)**2+(v4km-vmean)**2)
srmean_ulm = (sr6m+sr7m+sr8m+sr9m+sr10m)/5
srh2=((u500m-ubunkers)*(vsfc-vbunkers)-(usfc-ubunkers)*(v500m-vbunkers))
srh3=((u1km-ubunkers)*(v500m-vbunkers)-(u500m-ubunkers)*(v1km-vbunkers))
srh4=((u1_5km-ubunkers)*(v1km-vbunkers)-(u1km-ubunkers)*(v1_5km-vbunkers))
srh5=((u2km-ubunkers)*(v1_5km-vbunkers)-(u1_5km-ubunkers)*(v2km-vbunkers))
srh6=((u2_5km-ubunkers)*(v2km-vbunkers)-(u2km-ubunkers)*(v2_5km-vbunkers))
srh7=((u3km-ubunkers)*(v2_5km-vbunkers)-(u2_5km-ubunkers)*(v3km-vbunkers))
srh1a=((u250m-ubunkers)*(vsfc-vbunkers)-(usfc-ubunkers)*(v250m-vbunkers))
srh1b=((u500m-ubunkers)*(v250m-vbunkers)-(u250m-ubunkers)*(v500m-vbunkers))
#srh8=((u250m-ubunkers)*(vsfc-vbunkers)-(usfc-ubunkers)*(v250m-vbunkers))
#srh9=((u500m-ubunkers)*(v250m-vbunkers)-(u250m-ubunkers)*(v500m-vbunkers))
#srh10=((u750m-ubunkers)*(v500m-vbunkers)-(u500m-ubunkers)*(v750m-vbunkers))
#srh11=((u1km-ubunkers)*(v750m-vbunkers)-(u750m-ubunkers)*(v1km-vbunkers))
srh3km=srh2+srh3+srh4+srh5+srh6+srh7
srh1km=srh2+srh3
srh500m=srh1a+srh1b
wmaxshear = dls*np.sqrt(2*cape)
ehi3 = (cape*srh3km)/160000
ehi1 = (cape*srh1km)/160000
dx, dy = metpy.xarray.grid_deltas_from_dataarray(u500m)
div500 = mpcalc.divergence(u500m, v500m)
conv500 = np.ma.masked_where(div500>0,div500)
#cin_new = data_sfc['cin'].where(data_sfc['cin']<0.01,-0.01)
#print(conv500)
#np.max(cin)
Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable
2023-07-13 00:00:00
def plot_background(ax):
ax.set_extent([2, 9, 48, 51.5])
ax.add_feature(cfeature.COASTLINE.with_scale('10m'), LineWidth=2)
ax.add_feature(cfeature.BORDERS.with_scale('10m'),LineWidth=2)
gl = ax.gridlines(draw_labels=True,linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 12, 'color': 'black'}
gl.ylabel_style = {'size': 12, 'color': 'black'}
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return ax
#import matplotlib
#cmap = matplotlib.cm.get_cmap('cubehelix_r')
#for i in range(20):
#rgba = cmap(i)
# rgb2hex accepts rgb or rgba
#print(rgba)
# Create the figure and plot background on different axes
for i in range(49):
crs = ccrs.Mercator()
fig, axarr = plt.subplots(nrows=2, ncols=2, figsize=(20, 15), constrained_layout=False,
subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
plot_background(ax)
clevs_cape = np.arange(100,5100,400)
clevs_srh3 = np.arange(50,650,50)
clevs_lls = np.arange(4,32,2)
clevs_cin = [-25,0,25,50,75,100,125,150,175,200,225,250]
clevs_conv = np.arange(-100,0,10)
# cmap = plt.get_cmap('gist_ncar')
# newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
cmap = plt.get_cmap('gist_ncar')
newcmap1 = ListedColormap(cmap(np.linspace(0.15, 0.8, 31)))
colors=[(0.0, 0.9254901960784314, 0.9254901960784314),
(0.00392156862745098, 0.6274509803921569, 0.9647058823529412),
(0.0, 0.0, 0.9647058823529412),
(0.0, 1.0, 0.0),
(0.0, 0.7843137254901961, 0.0),
(0.0, 0.5647058823529412, 0.0),
(1.0, 1.0, 0.0),
(0.9058823529411765, 0.7529411764705882, 0.0),
(1.0, 0.5647058823529412, 0.0),
(1.0, 0.16078431372, 0.16078431372),
(0.7529411764705882, 0.0, 0.0),
(0.59765625, 0.0, 0.0),
(1.0, 0.0, 1.0),
(0.6, 0.3333333333333333, 0.788235294117647),
(0.27,0,0.4)]
colors2 = [(1,1,1),
(0.953,0.882,0.882),
(0.933,0.784,0.792),
(0.906,0.686,0.694),
(0.871,0.584,0.592),
(0.831,0.478,0.494),
(0.784,0.369,0.388),
(0.729,0.251,0.278),
(0.643,0.141,0.184),
(0.525,0.075,0.122),
(0.412,0.000,0.047)]
cmap1 = ListedColormap(colors)
cmap2 = ListedColormap(colors2)
newcmap = ListedColormap(cmap1(np.linspace(0, 1, 19)))
time = data_sfc['prmsl'].metpy.time
data_crs = data_sfc['prmsl'].metpy.cartopy_crs
# Upper left plot
cf1 = axlist[0].contourf(data_ml.longitude, data_ml.latitude, cape.metpy.loc[{'time': time[i]}],
clevs_cape, cmap=newcmap, extend='max', transform=data_crs)
wind_slice = (slice(None, None, 7), slice(None, None, 7))
c1=axlist[0].barbs(x[wind_slice[0]], y[wind_slice[1]],
udls.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
vdls.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
pivot='middle', color='red',transform=data_crs, length=7,
zorder=5)
#axlist[0].clabel(c1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[0].set_title('ML CAPE & 0-6 km Shear Vector', fontsize=16)
cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical',
ticks=(100,500,900,1300,1700,2100,2500,2900,3300,3700,4100,4500),
shrink=1, fraction=0.05, pad=0)
cb1.set_label('J/kg', size='x-large')
# Upper left plot
cf2 = axlist[1].contourf(data_ml.longitude, data_ml.latitude, lls.metpy.loc[{'time': time[i]}],
clevs_lls, cmap='BuPu', extend='max', transform=data_crs)
wind_slice = (slice(None, None,7), slice(None, None, 7))
c2=axlist[1].barbs(x[wind_slice[0]], y[wind_slice[1]],
umls.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
vmls.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
pivot='middle', color='red',transform=data_crs, length=7,
zorder=5)
ccf2 = axlist[1].contour(data_ml.longitude, data_ml.latitude, lls.metpy.loc[{'time': time[i]}],
[10,15], colors='blue' , transform=data_crs)
#axlist[0].clabel(c1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[1].set_title('0-1 km Bulk Shear & 0-3 km Shear Vector', fontsize=16)
cb2 = fig.colorbar(cf2, ax=axlist[1], orientation='vertical',
ticks=(4,6,8,10,12,14,16,18,20,22,24,26,28,30),shrink=1, fraction=0.05, pad=0)
cb2.set_label('m/s', size='x-large')
axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
# Upper left plot
cf3 = axlist[2].contourf(data_ml.longitude, data_ml.latitude, div500[i,:,:]*100000,
clevs_conv, cmap='Reds_r', extend='min', transform=data_crs)
ccf3= axlist[2].contour(data_ml.longitude, data_ml.latitude, div500[i,:,:]*100000,
[-100,-50,-30,-10], colors='blue', linestyles='dashed', transform=data_crs)
axlist[2].clabel(ccf3, fontsize=10, fmt='%i', rightside_up=True)
axlist[2].set_title('500 m Convergence (≤ -10x10⁻⁵ s⁻¹)', fontsize=16)
cb3 = fig.colorbar(cf3, ax=axlist[2], orientation='vertical'
,shrink=1, fraction=0.05, pad=0)
cb3.set_label('x10⁻⁵ s⁻¹', size='x-large')
# Upper left plot
cf4 = axlist[3].contourf(data_ml.longitude, data_ml.latitude, srh3km.metpy.loc[{'time': time[i]}],
clevs_srh3, cmap='plasma_r', extend='max', transform=data_crs)
ccf4 = axlist[3].contour(data_ml.longitude, data_ml.latitude, srh1km.metpy.loc[{'time': time[i]}],
[100,150,200,300], colors='black' , linestyles='dotted', transform=data_crs)
axlist[3].clabel(ccf4, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[3].set_title('0-3 km SRH & 0-1 km SRH (dotted lines)', fontsize=16)
cb4 = fig.colorbar(cf4, ax=axlist[3], orientation='vertical',
ticks=(50,100,150,200,250,300,350,400,450,500,550,600),shrink=1, fraction=0.05, pad=0)
cb4.set_label('m²/s²', size='x-large')
# Set height padding for plots
fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
# Set figure title
plt.gcf().text(0.13, 0.94, 'Model: ICON-EU 0.0625° | '+timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+data_ml['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
#plt.gcf().text(0.4, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC'), fontsize=20)
#plt.gcf().text(0.5, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC / ')+data_ml['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
plt.gcf().text(0.55, 0.09, 'Note: SRH is calculated for a right-moving storm cell.', fontsize=10)
# Display the plot
time2 = str(i)
base_filename='icon_comp_'
suffix='.jpeg'
latest='latest'
my_file = base_filename+time2+suffix
print(my_file)
plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
plt.close(fig)
icon_comp_0.jpeg
/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later c = matplotlib.collections.PathCollection(paths,
icon_comp_1.jpeg icon_comp_2.jpeg icon_comp_3.jpeg icon_comp_4.jpeg icon_comp_5.jpeg icon_comp_6.jpeg icon_comp_7.jpeg icon_comp_8.jpeg icon_comp_9.jpeg icon_comp_10.jpeg icon_comp_11.jpeg icon_comp_12.jpeg icon_comp_13.jpeg icon_comp_14.jpeg icon_comp_15.jpeg icon_comp_16.jpeg icon_comp_17.jpeg icon_comp_18.jpeg icon_comp_19.jpeg icon_comp_20.jpeg icon_comp_21.jpeg icon_comp_22.jpeg icon_comp_23.jpeg icon_comp_24.jpeg icon_comp_25.jpeg icon_comp_26.jpeg icon_comp_27.jpeg icon_comp_28.jpeg icon_comp_29.jpeg icon_comp_30.jpeg icon_comp_31.jpeg icon_comp_32.jpeg icon_comp_33.jpeg icon_comp_34.jpeg icon_comp_35.jpeg icon_comp_36.jpeg icon_comp_37.jpeg icon_comp_38.jpeg icon_comp_39.jpeg icon_comp_40.jpeg icon_comp_41.jpeg icon_comp_42.jpeg icon_comp_43.jpeg icon_comp_44.jpeg icon_comp_45.jpeg icon_comp_46.jpeg icon_comp_47.jpeg icon_comp_48.jpeg
# Create the figure and plot background on different axes
for i in range(49):
crs = ccrs.Mercator()
fig, axarr = plt.subplots(nrows=2, ncols=2, figsize=(20, 15), constrained_layout=False,
subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
plot_background(ax)
clevs_wmax = np.arange(-100,2500,200)
clevs_tcwv = np.arange(0,55,3)
clevs_sm = np.arange(0,39,3)
#clevs_cin = [-50,0,50,100,150,200,250,300,350,400,450,500]
clevs_qv = np.arange(0,20,2)
clevs_mslp = np.arange(950, 1050, 1)
#clevs_ehi = np.arange(-0.5, 6.5,0.5)
cmap3 = plt.get_cmap('cubehelix_r')
newcmap3 = ListedColormap(cmap3(np.linspace(0, 0.9, 20)))
bounds = [0,0.1,0.2,0.5,1,2,3,4,5]
norm = BoundaryNorm(bounds, newcmap3.N)
# cmap = plt.get_cmap('gist_ncar')
# newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
cmap = plt.get_cmap('gist_ncar')
newcmap1 = ListedColormap(cmap(np.linspace(0.15, 0.8, 30)))
colors=[(1,1,1),(0.0, 0.9254901960784314, 0.9254901960784314),
(0.00392156862745098, 0.6274509803921569, 0.9647058823529412),
(0.0, 0.0, 0.9647058823529412),
(0.0, 1.0, 0.0),
(0.0, 0.7843137254901961, 0.0),
(0.0, 0.5647058823529412, 0.0),
(1.0, 1.0, 0.0),
(0.9058823529411765, 0.7529411764705882, 0.0),
(1.0, 0.5647058823529412, 0.0),
(1.0, 0.16078431372, 0.16078431372),
(0.7529411764705882, 0.0, 0.0),
(0.59765625, 0.0, 0.0),
(1.0, 0.0, 1.0),
(0.6, 0.3333333333333333, 0.788235294117647),
(0.27,0,0.4)]
colors2 = [(1,1,1),
(0.953,0.882,0.882),
(0.933,0.784,0.792),
(0.906,0.686,0.694),
(0.871,0.584,0.592),
(0.831,0.478,0.494),
(0.784,0.369,0.388),
(0.729,0.251,0.278),
(0.643,0.141,0.184),
(0.525,0.075,0.122),
(0.412,0.000,0.047)]
cmap1 = ListedColormap(colors)
cmap2 = ListedColormap(colors2)
newcmap = ListedColormap(cmap1(np.linspace(0, 1, 17)))
time = data_ml['u'].metpy.time
data_crs = data_ml['u'].metpy.cartopy_crs
# Upper left plot
cf1 = axlist[0].contourf(data_ml.longitude, data_ml.latitude, wmaxshear.metpy.loc[{'time': time[i]}],
clevs_wmax, cmap=newcmap, extend='max', transform=data_crs)
#axlist[0].clabel(c1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[0].set_title('WMAXSHEAR', fontsize=16)
cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical',
ticks=(100,300,500,700,900,1100,1300,1500,1700,1900,2100),
shrink=1, fraction=0.05, pad=0)
cb1.set_label('m²/s²', size='x-large')
plt.gcf().text(0.13, 0.50, 'Note: WMAXSHEAR is calculated as follows: sqrt(2 x ML CAPE) x (0-6 km bulk shear)', fontsize=10)
# Upper left plot
cf2 = axlist[1].contourf(data_ml.longitude, data_ml.latitude, ehi3.metpy.loc[{'time': time[i]}],
[0,0.1,0.2,0.5,1,2,3,4,5], cmap=newcmap3, norm=norm, transform=data_crs)
#ccf2 = axlist[1].contour(data_ml.longitude, data_ml.latitude, ehi1.metpy.loc[{'time': time[i]}],
#[0.1,1], colors='blue' , transform=data_crs)
#axlist[0].clabel(c1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[1].set_title('0-3 km EHI', fontsize=16)
cb2 = fig.colorbar(cf2, ax=axlist[1], orientation='vertical', shrink=1, fraction=0.05, pad=0)
#cb2.set_label('', size='x-large')
#cb2.set_ticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
#axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
plt.gcf().text(0.55, 0.50, 'Note: EHI is calculated as follows: (ML CAPE x SRH)/(1.6 x 10⁵)', fontsize=10)
# Upper left plot
cf3 = axlist[2].contourf(data_ml.longitude, data_ml.latitude, qmean.metpy.loc[{'time': time[i]}]*1000,
clevs_qv, cmap='YlGnBu', extend='max', transform=data_crs)
ccf3= axlist[2].contour(data_ml.longitude, data_ml.latitude, mslp.metpy.loc[{'time': time[i]}],
clevs_mslp, colors='black', transform=data_crs)
cccf3 = axlist[2].contour(data_ml.longitude, data_ml.latitude, qmean.metpy.loc[{'time': time[i]}]*1000,
[6,10,14], colors='black' , linestyles='dotted', transform=data_crs)
axlist[2].clabel(ccf3, fontsize=10, fmt='%i', rightside_up=True)
axlist[2].clabel(cccf3, fontsize=10, fmt='%i', rightside_up=True)
axlist[2].set_title('MSLP & 0-500 m Mean Specific Humidity', fontsize=16)
cb3 = fig.colorbar(cf3, ax=axlist[2], orientation='vertical', shrink=1, fraction=0.05, pad=0)
cb3.set_label('g/kg', size='x-large')
# Upper left plot
cf4 = axlist[3].contourf(data_ml.longitude, data_ml.latitude, smotion.metpy.loc[{'time': time[i]}],
clevs_sm, cmap='pink_r', extend='max', transform=data_crs)
#ccf4 = axlist[3].contour(data_ml.longitude, data_ml.latitude, tcwv.metpy.loc[{'time': time[i]}],
#[10,20,30,40], colors='black' , linestyles='dotted', transform=data_crs)
#axlist[3].clabel(ccf4, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[3].set_title('Storm Motion Vector', fontsize=16)
wind_slice = (slice(None, None, 7), slice(None, None, 7))
c4=axlist[3].barbs(x[wind_slice[0]], y[wind_slice[1]],
umean.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
vmean.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
pivot='middle', color='blue',transform=data_crs, length=7,
zorder=3)
cb4 = fig.colorbar(cf4, ax=axlist[3], orientation='vertical',shrink=1, fraction=0.05, pad=0)
cb4.set_label('m/s', size='x-large')
# Set height padding for plots
fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
# Set figure title
plt.gcf().text(0.13, 0.94, 'Model: ICON-EU 0.0625° | '+timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+data_ml['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
#plt.gcf().text(0.4, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC'), fontsize=20)
#plt.gcf().text(0.5, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC / ')+data_sfc['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
plt.gcf().text(0.55, 0.09, 'Note: Motion of a non-deviant storm cell.', fontsize=10)
# Display the plot
time2 = str(i)
base_filename='icon_comp2_'
suffix='.jpeg'
#latest='latest'
my_file = base_filename+time2+suffix
print(my_file)
plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
plt.close(fig)
icon_comp2_0.jpeg icon_comp2_1.jpeg icon_comp2_2.jpeg icon_comp2_3.jpeg icon_comp2_4.jpeg icon_comp2_5.jpeg icon_comp2_6.jpeg icon_comp2_7.jpeg icon_comp2_8.jpeg icon_comp2_9.jpeg icon_comp2_10.jpeg icon_comp2_11.jpeg icon_comp2_12.jpeg icon_comp2_13.jpeg icon_comp2_14.jpeg icon_comp2_15.jpeg icon_comp2_16.jpeg icon_comp2_17.jpeg icon_comp2_18.jpeg icon_comp2_19.jpeg icon_comp2_20.jpeg icon_comp2_21.jpeg icon_comp2_22.jpeg icon_comp2_23.jpeg icon_comp2_24.jpeg icon_comp2_25.jpeg icon_comp2_26.jpeg icon_comp2_27.jpeg icon_comp2_28.jpeg icon_comp2_29.jpeg icon_comp2_30.jpeg icon_comp2_31.jpeg icon_comp2_32.jpeg icon_comp2_33.jpeg icon_comp2_34.jpeg icon_comp2_35.jpeg icon_comp2_36.jpeg icon_comp2_37.jpeg icon_comp2_38.jpeg icon_comp2_39.jpeg icon_comp2_40.jpeg icon_comp2_41.jpeg icon_comp2_42.jpeg icon_comp2_43.jpeg icon_comp2_44.jpeg icon_comp2_45.jpeg icon_comp2_46.jpeg icon_comp2_47.jpeg icon_comp2_48.jpeg
# Create the figure and plot background on different axes
for i in range(49):
crs = ccrs.Mercator()
fig, axarr = plt.subplots(nrows=2, ncols=2, figsize=(20, 15), constrained_layout=False,
subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
plot_background(ax)
clevs_wmax = np.arange(-100,2500,200)
clevs_tcwv = np.arange(0,55,3)
clevs_sm = np.arange(0,39,3)
clevs_lls = np.arange(0,24,2)
#clevs_cin = [-50,0,50,100,150,200,250,300,350,400,450,500]
clevs_qv = np.arange(0,20,2)
clevs_mslp = np.arange(950, 1050, 1)
#clevs_ehi = np.arange(-0.5, 6.5,0.5)
cmap3 = plt.get_cmap('cubehelix_r')
newcmap3 = ListedColormap(cmap3(np.linspace(0, 0.9, 20)))
bounds = [0,0.1,0.2,0.5,1,2,3,4,5]
norm = BoundaryNorm(bounds, newcmap3.N)
# cmap = plt.get_cmap('gist_ncar')
# newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
cmap = plt.get_cmap('gist_ncar')
newcmap1 = ListedColormap(cmap(np.linspace(0.15, 0.8, 30)))
colors=[(1,1,1),(0.0, 0.9254901960784314, 0.9254901960784314),
(0.00392156862745098, 0.6274509803921569, 0.9647058823529412),
(0.0, 0.0, 0.9647058823529412),
(0.0, 1.0, 0.0),
(0.0, 0.7843137254901961, 0.0),
(0.0, 0.5647058823529412, 0.0),
(1.0, 1.0, 0.0),
(0.9058823529411765, 0.7529411764705882, 0.0),
(1.0, 0.5647058823529412, 0.0),
(1.0, 0.16078431372, 0.16078431372),
(0.7529411764705882, 0.0, 0.0),
(0.59765625, 0.0, 0.0),
(1.0, 0.0, 1.0),
(0.6, 0.3333333333333333, 0.788235294117647),
(0.27,0,0.4)]
colors2 = [(1,1,1),
(0.953,0.882,0.882),
(0.933,0.784,0.792),
(0.906,0.686,0.694),
(0.871,0.584,0.592),
(0.831,0.478,0.494),
(0.784,0.369,0.388),
(0.729,0.251,0.278),
(0.643,0.141,0.184),
(0.525,0.075,0.122),
(0.412,0.000,0.047)]
cmap1 = ListedColormap(colors)
cmap2 = ListedColormap(colors2)
newcmap = ListedColormap(cmap1(np.linspace(0, 1, 17)))
time = data_ml['u'].metpy.time
data_crs = data_ml['u'].metpy.cartopy_crs
# Upper left plot
cf1 = axlist[0].contourf(data_ml.longitude, data_ml.latitude, srmean_ll.metpy.loc[{'time': time[i]}],
clevs_lls, cmap='BuPu', extend='max',transform=data_crs)
ccf1 = axlist[0].contour(data_ml.longitude, data_ml.latitude, srmean_ll.metpy.loc[{'time': time[i]}],
[5,10,15,20], colors='black' , linestyles='dotted', transform=data_crs)
axlist[0].clabel(ccf1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[0].set_title('0-2 km Mean Storm-Relative Wind (RM)', fontsize=16)
cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical',
shrink=1, fraction=0.05, pad=0)
cb1.set_label('m/s', size='x-large')
#plt.gcf().text(0.13, 0.50, 'Note: WMAXSHEAR is calculated as follows: sqrt(2 x MU CAPE) x (0-6 km bulk shear)', fontsize=10)
# Upper left plot
cf2 = axlist[1].contourf(data_ml.longitude, data_ml.latitude, srmean_llm.metpy.loc[{'time': time[i]}],
clevs_lls, cmap='BuPu', extend='max', transform=data_crs)
ccf2 = axlist[1].contour(data_ml.longitude, data_ml.latitude, srmean_llm.metpy.loc[{'time': time[i]}],
[5,10,15,20], colors='black' , linestyles='dotted', transform=data_crs)
axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[1].set_title('0-2 km Mean Storm-Relative Wind (non-deviant motion)', fontsize=16)
cb2 = fig.colorbar(cf2, ax=axlist[1], orientation='vertical', shrink=1, fraction=0.05, pad=0)
cb2.set_label('m/s', size='x-large')
#cb2.set_ticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
#axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
#plt.gcf().text(0.55, 0.50, 'Note: EHI is calculated as follows: (MU CAPE x SRH)/(1.6 x 10⁵)', fontsize=10)
# Upper left plot
cf3 = axlist[2].contourf(data_ml.longitude, data_ml.latitude, srmean_ul.metpy.loc[{'time': time[i]}],
clevs_lls, cmap='BuPu', extend='max', transform=data_crs)
ccf3= axlist[2].contour(data_ml.longitude, data_ml.latitude, srmean_ul.metpy.loc[{'time': time[i]}],
[5,10,15,20], colors='black', linestyles='dotted', transform=data_crs)
#cccf3 = axlist[2].contour(data_ml.longitude, data_ml.latitude, qmean.metpy.loc[{'time': time[i]}]*1000,
#[6,10,14], colors='black' , linestyles='dotted', transform=data_crs)
axlist[2].clabel(ccf3, fontsize=10, fmt='%i', rightside_up=True)
#axlist[2].clabel(cccf3, fontsize=10, fmt='%i', rightside_up=True)
axlist[2].set_title('2-4 km Mean Storm-Relative Wind (RM)', fontsize=16)
cb3 = fig.colorbar(cf3, ax=axlist[2], orientation='vertical',
shrink=1, fraction=0.05, pad=0)
cb3.set_label('m/s', size='x-large')
# Upper left plot
cf4 = axlist[3].contourf(data_ml.longitude, data_ml.latitude, srmean_ulm.metpy.loc[{'time': time[i]}],
clevs_lls, cmap='BuPu', extend='max',transform=data_crs)
ccf4 = axlist[3].contour(data_ml.longitude, data_ml.latitude, srmean_ulm.metpy.loc[{'time': time[i]}],
[5,10,15,20], colors='black' , linestyles='dotted', transform=data_crs)
axlist[3].clabel(ccf4, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[3].set_title('2-4 km Mean Storm-Relative Wind (non-deviant motion)', fontsize=16)
#wind_slice = (slice(None, None, 4), slice(None, None, 4))
#c4=axlist[3].barbs(x[wind_slice[0]], y[wind_slice[1]],
#umean.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
#vmean.metpy.loc[{'time': time[i]}].values[wind_slice[0], wind_slice[1]]*1.94384,
#pivot='middle', color='red',transform=data_crs, length=7,
#zorder=3)
cb4 = fig.colorbar(cf4, ax=axlist[3],
orientation='vertical',shrink=1, fraction=0.05, pad=0)
cb4.set_label('m/s', size='x-large')
# Set height padding for plots
fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
# Set figure title
plt.gcf().text(0.13, 0.94, 'Model: ICON-EU 0.0625° | '+timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+data_ml['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
#plt.gcf().text(0.4, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC'), fontsize=20)
#plt.gcf().text(0.5, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC / ')+data_sfc['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=20)
#plt.gcf().text(0.55, 0.09, 'Note: Motion of a non-deviant storm cell.', fontsize=10)
# Display the plot
time2 = str(i)
base_filename='icon_srwind_'
suffix='.jpeg'
#latest='latest'
my_file = base_filename+time2+suffix
print(my_file)
plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
plt.close(fig)
icon_srwind_0.jpeg icon_srwind_1.jpeg icon_srwind_2.jpeg icon_srwind_3.jpeg icon_srwind_4.jpeg icon_srwind_5.jpeg icon_srwind_6.jpeg icon_srwind_7.jpeg icon_srwind_8.jpeg icon_srwind_9.jpeg icon_srwind_10.jpeg icon_srwind_11.jpeg icon_srwind_12.jpeg icon_srwind_13.jpeg icon_srwind_14.jpeg icon_srwind_15.jpeg icon_srwind_16.jpeg icon_srwind_17.jpeg icon_srwind_18.jpeg icon_srwind_19.jpeg icon_srwind_20.jpeg icon_srwind_21.jpeg icon_srwind_22.jpeg icon_srwind_23.jpeg icon_srwind_24.jpeg icon_srwind_25.jpeg icon_srwind_26.jpeg icon_srwind_27.jpeg icon_srwind_28.jpeg icon_srwind_29.jpeg icon_srwind_30.jpeg icon_srwind_31.jpeg icon_srwind_32.jpeg icon_srwind_33.jpeg icon_srwind_34.jpeg icon_srwind_35.jpeg icon_srwind_36.jpeg icon_srwind_37.jpeg icon_srwind_38.jpeg icon_srwind_39.jpeg icon_srwind_40.jpeg icon_srwind_41.jpeg icon_srwind_42.jpeg icon_srwind_43.jpeg icon_srwind_44.jpeg icon_srwind_45.jpeg icon_srwind_46.jpeg icon_srwind_47.jpeg icon_srwind_48.jpeg