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 metpy.plots import ctables
from matplotlib.colors import Normalize
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, BoundaryNorm
import scipy
import cfgrib
import matplotlib
# configure backend here
matplotlib.use('Agg')
file_dir = '/data/icond2'
# Changing the directory
os.chdir(file_dir)
#data = xr.open_dataset('ec_pl3_small.nc')
data_dbz = xr.open_dataset('icond2_dbz.nc').sel(longitude=slice(-1,15,1),latitude=slice(45,55,1))
#data_gust = xr.open_dataset('icond2_vmax.nc').sel(longitude=slice(-1,15,1),latitude=slice(45,55,1))
#data_uh = xr.open_dataset('icond2_uh_max.nc').sel(longitude=slice(-1,15,1),latitude=slice(45,55,1))
data_sat = xr.open_dataset('icond2_ir108.nc').sel(longitude=slice(-1,15,1),latitude=slice(45,55,1))
# View a summary of the Dataset
print(data_dbz)
print(data_sat)
<xarray.Dataset> Dimensions: (longitude: 801, latitude: 501, time: 48) Coordinates: * longitude (longitude) float32 -1.0 -0.98 -0.96 -0.94 ... 14.96 14.98 15.0 * latitude (latitude) float32 45.0 45.02 45.04 45.06 ... 54.96 54.98 55.0 * time (time) datetime64[ns] 2023-07-13T04:00:00 ... 2023-07-15T03:00:00 Data variables: bref (time, latitude, longitude) float32 ... Attributes: Conventions: CF-1.6 history: 2023-07-13 04:43:10 GMT by grib_to_netcdf-2.6.0: grib_to_ne... <xarray.Dataset> Dimensions: (longitude: 801, latitude: 501, time: 48) Coordinates: * longitude (longitude) float32 -1.0 -0.98 -0.96 -0.94 ... 14.96 14.98 15.0 * latitude (latitude) float32 45.0 45.02 45.04 45.06 ... 54.96 54.98 55.0 * time (time) datetime64[ns] 2023-07-13T04:00:00 ... 2023-07-15T03:00:00 Data variables: clbt (time, latitude, longitude) float32 ... Attributes: Conventions: CF-1.6 history: 2023-07-13 04:43:12 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_dbz = data_dbz.metpy.parse_cf()
#data_uh = data_uh.metpy.parse_cf()
data_sat = data_sat.metpy.parse_cf()
x, y = data_sat['clbt'].metpy.coordinates('x', 'y')
time = data_sat['clbt'].metpy.time
time2 = data_dbz['bref'].metpy.time
timeinit = time[0]
timeinit = datetime.utcfromtimestamp(timeinit.item()/1e9)-timedelta(hours=1)
print(timeinit)
timeinit2 = time2[0]
timeinit2 = datetime.utcfromtimestamp(timeinit2.item()/1e9)-timedelta(hours=1)
print(timeinit2)
#uh = data_uh['UH_MAX']
dbz = data_dbz['bref']
clbt = data_sat['clbt']
dbz.data = np.nan_to_num(dbz.data, copy=True, nan=-1)
#print(dbz.data)
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 03:00:00 2023-07-13 03: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)
cmap = ctables.colortables.get_colortable('NWSStormClearReflectivity')
newcmap = ListedColormap(cmap(np.linspace(0.25, 0.92, 28)))
cmap2 = ctables.colortables.get_colortable('NWSReflectivity')
newcmap2 = ListedColormap(cmap2(np.linspace(0.2, 0.96, 28)))
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.388235, 0.462745, 0.658824), (0.372549, 0.45098, 0.654902), (0.372549, 0.45098, 0.654902),
(0.356863, 0.439216, 0.65098), (0.341176, 0.427451, 0.643137), (0.32549, 0.415686, 0.639216),
(0.309804, 0.403922, 0.635294), (0.294118, 0.392157, 0.631373), (0.278431, 0.380392, 0.627451),
(0.262745, 0.368627, 0.623529), (0.262745, 0.380392, 0.635294), (0.270588, 0.407843, 0.65098),
(0.282353, 0.435294, 0.666667), (0.290196, 0.462745, 0.682353), (0.301961, 0.490196, 0.698039),
(0.309804, 0.517647, 0.713725), (0.317647, 0.545098, 0.733333), (0.329412, 0.572549, 0.74902),
(0.337255, 0.6, 0.764706), (0.34902, 0.623529, 0.780392), (0.356863, 0.65098, 0.796078),
(0.368627, 0.678431, 0.811765), (0.376471, 0.705882, 0.831373), (0.384314, 0.733333, 0.847059),
(0.396078, 0.760784, 0.862745), (0.403922, 0.788235, 0.878431), (0.415686, 0.815686, 0.894118),
(0.435294, 0.839216, 0.909804), (0.407843, 0.839216, 0.843137), (0.380392, 0.839216, 0.772549),
(0.34902, 0.839216, 0.701961), (0.321569, 0.839216, 0.635294), (0.294118, 0.839216, 0.564706),
(0.262745, 0.839216, 0.494118), (0.235294, 0.839216, 0.427451), (0.207843, 0.839216, 0.356863),
(0.066667, 0.835294, 0.094118), (0.066667, 0.819608, 0.090196), (0.062745, 0.803922, 0.090196),
(0.062745, 0.784314, 0.086275), (0.062745, 0.768627, 0.086275), (0.058824, 0.752941, 0.082353),
(0.058824, 0.737255, 0.082353), (0.058824, 0.717647, 0.078431), (0.054902, 0.701961, 0.078431),
(0.054902, 0.686275, 0.07451), (0.054902, 0.670588, 0.07451), (0.05098, 0.65098, 0.070588),
(0.05098, 0.635294, 0.070588), (0.05098, 0.619608, 0.066667), (0.047059, 0.6, 0.066667),
(0.047059, 0.584314, 0.062745), (0.047059, 0.568627, 0.062745), (0.043137, 0.552941, 0.058824),
(0.043137, 0.533333, 0.058824), (0.043137, 0.517647, 0.054902), (0.039216, 0.501961, 0.054902),
(0.039216, 0.486275, 0.05098), (0.039216, 0.466667, 0.05098), (0.035294, 0.45098, 0.047059),
(0.035294, 0.435294, 0.047059), (0.035294, 0.419608, 0.043137), (0.031373, 0.4, 0.043137),
(0.031373, 0.384314, 0.039216), (0.035294, 0.368627, 0.035294), (0.113725, 0.407843, 0.035294),
(0.196078, 0.45098, 0.031373), (0.27451, 0.490196, 0.031373), (0.356863, 0.533333, 0.027451),
(0.435294, 0.572549, 0.027451), (0.517647, 0.615686, 0.023529), (0.596078, 0.658824, 0.023529),
(0.678431, 0.698039, 0.019608), (0.756863, 0.741176, 0.019608), (0.839216, 0.780392, 0.015686),
(0.917647, 0.823529, 0.015686), (1.0, 0.886275, 0.0), (1.0, 0.847059, 0.0), (1.0, 0.827451, 0.0),
(1.0, 0.788235, 0.0), (1.0, 0.768627, 0.0), (1.0, 0.733333, 0.0), (1.0, 0.713725, 0.0),
(1.0, 0.693725, 0.0), (1.0, 0.67451, 0.0), (1.0, 0.654902, 0.0), (1.0, 0.619608, 0.0),
(1.0, 0.6, 0.0), (1.0, 0.580392, 0.0), (1.0, 0.541176, 0.0), (1.0, 0.521569, 0.0),
(1.0, 0.501569, 0.0), (0.945098, 0.0, 0.0),
(0.917647, 0.0, 0.0), (0.890196, 0.0, 0.0), (0.862745, 0.0, 0.0), (0.835294, 0.0, 0.0),
(0.803922, 0.0, 0.0), (0.776471, 0.0, 0.0), (0.74902, 0.0, 0.0), (0.721569, 0.0, 0.0),
(0.694118, 0.0, 0.0), (0.666667, 0.0, 0.0), (0.639216, 0.0, 0.0), (0.607843, 0.0, 0.0),
(0.580392, 0.0, 0.0), (0.552941, 0.0, 0.0), (0.52549, 0.0, 0.0), (0.498039, 0.0, 0.0),
(0.470588, 0.0, 0.0), (0.443137, 0.0, 0.0), (1.0, 0.960784, 1.0),
(1.0, 0.917647, 1.0), (1.0, 0.87451, 1.0), (1.0, 0.831373, 1.0), (1.0, 0.788235, 1.0),
(1.0, 0.745098, 1.0), (1.0, 0.701961, 1.0), (1.0, 0.658824, 1.0), (1.0, 0.615686, 1.0),
(1.0, 0.572549, 1.0), (1.0, 0.458824, 1.0), (0.988235, 0.419608, 0.992157),
(0.976471, 0.376471, 0.980392), (0.964706, 0.337255, 0.968627), (0.952941, 0.294118, 0.956863),
(0.941176, 0.25098, 0.945098), (0.929412, 0.211765, 0.937255), (0.917647, 0.168627, 0.92549),
(0.905882, 0.12549, 0.913725), (0.894118, 0.086275, 0.901961), (0.882353, 0.043137, 0.890196),
(0.698039, 0.0, 1.0), (0.67451, 0.0, 0.988235), (0.643137, 0.0, 0.968627), (0.607843, 0.0, 0.956863),
(0.576471, 0.0, 0.937255), (0.533333, 0.0, 0.917647), (0.513725, 0.0, 0.909804),
(0.47451, 0.0, 0.886275), (0.447059, 0.0, 0.866667), (0.411765, 0.0, 0.858824),
(0.388235, 0.0, 0.839216)]
cmap1 = ListedColormap(colors)
newcmap1 = ListedColormap(cmap1(np.linspace(0, 1, 17)))
cmap2 = ListedColormap(colors2)
newcmap2 = ListedColormap(cmap2(np.linspace(0, 0.9, 29)))
colors3=[(35/255,35/255,35/255),(66/255,21/255,21/255),(126/255,1/255,1/255),(186/255,0/255,0/255),(255/255,45/255,0/255),
(255/255,155/255,0/255),(253/255,255/255,2/255),
(192/255,255/255,63/255),
(96/255,255/255,159/255),(0,249/255,255/255),(0,151/255,255/255),(0/255, 59/255, 255/255),
(0/255, 0/255, 190/255),(239/255, 239/255, 239/255),
(210/255, 210/255, 210/255),(167/255, 167/255, 167/255),(135/255, 135/255, 135/255),
(110/255, 110/255, 110/255),
(89/255, 89/255, 89/255),(65/255,65/255,65/255)]
cmap3 = ListedColormap(colors3)
newcmap3 = ListedColormap(cmap3(np.linspace(0, 1, 20)))
# Create the figure and plot background on different axes
clevs_sat = np.arange(-70,35,5)
clevs_dbz = np.arange(-2.5,72.5,2.5)
# cmap = plt.get_cmap('gist_ncar')
# newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
time = data_dbz['bref'].metpy.time
data_crs = data_dbz['bref'].metpy.time
crs = ccrs.Mercator()
for i in range(48):
fig, axarr = plt.subplots(nrows=1, ncols=2, figsize=(25, 10), constrained_layout=False,
subplot_kw={'projection': crs})
# Set height padding for plots
fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
axlist = axarr.flatten()
for ax in axlist:
plot_background(ax)
# Upper left plot
cf1 = axlist[0].contourf(dbz.longitude, dbz.latitude, dbz.metpy.loc[{'time': time[i]}],
clevs_dbz, cmap=newcmap2, transform=ccrs.PlateCarree())
#ccf1= axlist[0].contour(gust.longitude, gust.latitude, gust.metpy.loc[{'time': time[i]}]*3.6,
#clevs_gust, colors='darkviolet', transform=ccrs.PlateCarree())
#axlist[0].clabel(ccf1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[0].set_title('850 hPa Reflectivity', fontsize=16)
cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical', ticks=(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70),
shrink=0.81, fraction=0.1, pad=0)
cb1.set_label('dBZ', size='x-large')
cf2 = axlist[1].contourf(clbt.longitude, clbt.latitude, clbt.metpy.loc[{'time': time[i]}]-273.15,
clevs_sat, cmap=newcmap3, transform=ccrs.PlateCarree())
#ccf2= axlist[1].contour(gust.longitude, gust.latitude, gust.metpy.loc[{'time': time[i]}]*3.6,
#clevs_gust, colors='black', transform=ccrs.PlateCarree())
#axlist[0].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
axlist[1].set_title('Synthetic MSG SEVIRI Image at 10.8 μm', fontsize=16)
cb2= fig.colorbar(cf2, ax=axlist[1], orientation='vertical',ticks=(-70,-60,-50,-40,-30,-20,-10,0,10,20,30),
shrink=0.81, fraction=0.1, pad=0)
cb2.set_label('°C', size='x-large')
# Set figure title
plt.gcf().text(0.125, 0.88, 'Model: ICON-D2 0.02° | ' + timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ') + data_dbz['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values,fontsize=20)
#plt.gcf().text(0.55, 0.09, 'Note: Updraft helicity is vertically averaged between 500 and 6000 m AGL.', fontsize=10)
# Display the plot
time2 = str(i*1+1)
base_filename='icond2_comp2_'
suffix='.jpeg'
my_file = base_filename+time2+suffix
print(my_file)
plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
plt.close(fig)
icond2_comp2_1.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,
icond2_comp2_2.jpeg icond2_comp2_3.jpeg icond2_comp2_4.jpeg icond2_comp2_5.jpeg icond2_comp2_6.jpeg icond2_comp2_7.jpeg icond2_comp2_8.jpeg icond2_comp2_9.jpeg icond2_comp2_10.jpeg icond2_comp2_11.jpeg icond2_comp2_12.jpeg icond2_comp2_13.jpeg icond2_comp2_14.jpeg icond2_comp2_15.jpeg icond2_comp2_16.jpeg icond2_comp2_17.jpeg icond2_comp2_18.jpeg icond2_comp2_19.jpeg icond2_comp2_20.jpeg icond2_comp2_21.jpeg icond2_comp2_22.jpeg icond2_comp2_23.jpeg icond2_comp2_24.jpeg icond2_comp2_25.jpeg icond2_comp2_26.jpeg icond2_comp2_27.jpeg icond2_comp2_28.jpeg icond2_comp2_29.jpeg icond2_comp2_30.jpeg icond2_comp2_31.jpeg icond2_comp2_32.jpeg icond2_comp2_33.jpeg icond2_comp2_34.jpeg icond2_comp2_35.jpeg icond2_comp2_36.jpeg icond2_comp2_37.jpeg icond2_comp2_38.jpeg icond2_comp2_39.jpeg icond2_comp2_40.jpeg icond2_comp2_41.jpeg icond2_comp2_42.jpeg icond2_comp2_43.jpeg icond2_comp2_44.jpeg icond2_comp2_45.jpeg icond2_comp2_46.jpeg icond2_comp2_47.jpeg icond2_comp2_48.jpeg