In [1]:
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 wrf
import scipy
#import xcape
In [2]:
os.chdir('/data/ecmwf')
#data = xr.open_dataset('ec_pl3_small.nc')
#data_ml = xr.open_dataset('ec_ml.nc')
data_sfc = xr.open_dataset('ec_sfc2.nc')

# View a summary of the Dataset
#print(data_ml)
print(data_sfc)

<xarray.Dataset>
Dimensions:    (longitude: 151, latitude: 101, time: 48)
Coordinates:
  * longitude  (longitude) float32 0.0 0.1 0.2 0.3 0.4 ... 14.7 14.8 14.9 15.0
  * latitude   (latitude) float32 55.0 54.9 54.8 54.7 ... 45.3 45.2 45.1 45.0
  * time       (time) datetime64[ns] 2023-07-12T13:00:00 ... 2023-07-14T12:00:00
Data variables: (12/13)
    z          (time, latitude, longitude) float32 ...
    t2m        (time, latitude, longitude) float32 ...
    d2m        (time, latitude, longitude) float32 ...
    cin        (time, latitude, longitude) float32 ...
    tcwv       (time, latitude, longitude) float32 ...
    mld        (time, latitude, longitude) float32 ...
    ...         ...
    u10        (time, latitude, longitude) float32 ...
    v10        (time, latitude, longitude) float32 ...
    mlcape50   (time, latitude, longitude) float32 ...
    mlcin50    (time, latitude, longitude) float32 ...
    mucape     (time, latitude, longitude) float32 ...
    litota1    (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2023-07-12 18:07:18 GMT by grib_to_netcdf-2.30.2: grib_to_n...
In [3]:
# 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()

x, y = data_sfc['u10'].metpy.coordinates('x', 'y')

time = data_sfc['u10'].metpy.time

timeinit = time[0]
timeinit = datetime.utcfromtimestamp(timeinit.item()/1e9)-timedelta(hours=1)
print(timeinit)

sfcuwind =  data_sfc['u10']
sfcvwind =  data_sfc['v10']
mslp = data_sfc['msl']
mslp = mslp/100

#qv_bl = data_ml['q'].where(data_ml['q'].level>123, drop=True)
#qv_mean = np.mean(qv_bl, axis=1)
#qv_mean = qv_mean*1000

usfc = data_sfc['u10']
vsfc = data_sfc['v10']
cape = data_sfc['mucape']
mlcape = data_sfc['mlcape50']
mlcin = data_sfc['mlcin50']
cin = data_sfc['cin']
li = data_sfc['litota1']*0.041666  #convert from /(km2*day) to /(km2*h)
tcwv = data_sfc['tcwv']

cin_new = data_sfc['cin'].where(data_sfc['cin']<0.01,-0.01)
mlcin_new = data_sfc['mlcin50'].where(data_sfc['mlcin50']<0.01,-0.01)
#li_new = data_sfc['litota1'].where(data_sfc['litota1']==0.0,-0.1)*0.041666  #convert from /(km2*day) to /(km2*h)
np.max(li)

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
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-12 12:00:00
Out[3]:
<xarray.DataArray 'litota1' ()>
array(2.3796413, dtype=float32)
Coordinates:
    metpy_crs  object Projection: latitude_longitude
xarray.DataArray
'litota1'
  • 2.38
    array(2.3796413, dtype=float32)
    • metpy_crs
      ()
      object
      Projection: latitude_longitude
      array(<metpy.plots.mapping.CFProjection object at 0x7f0c8af7e8e0>,
      dtype=object)
In [4]:
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)
In [5]:
for i in range(48):
    
    crs = ccrs.Mercator()
    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)
        
    clevs_cape = np.arange(100,5100,400)
    # 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['mucape'].metpy.time
    data_crs = data_sfc['cin'].metpy.cartopy_crs
        
    # Upper left plot
    cf1 = axlist[0].contourf(data_sfc.longitude, data_sfc.latitude, cape.metpy.loc[{'time': time[i]}], 
                      clevs_cape, cmap=newcmap, extend='max', transform=data_crs)
    #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('Most-Unstable CAPE', 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=0.81, fraction=0.1, pad=0)
    cb1.set_label('J/kg', size='x-large')

    cf2 = axlist[1].contourf(data_sfc.longitude, data_sfc.latitude, mlcape.metpy.loc[{'time': time[i]}], 
                      clevs_cape, cmap=newcmap, extend='max', transform=data_crs)
    #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('50-hPa Mixed-Layer CAPE', fontsize=16)
    cb2= fig.colorbar(cf2, ax=axlist[1], orientation='vertical', ticks=(100,500,900,1300,1700,2100,2500,2900,3300,3700,4100,4500),
                              shrink=0.81, fraction=0.1, pad=0)
    cb2.set_label('J/kg²', size='x-large')


    # Set figure title
    plt.gcf().text(0.125, 0.88, 'Model: ECMWF IFS 0.1° | ' + 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.125, 0.13, 'Note: In the IFS "MU CAPE" refers to the most unstable parcel (the parcel with the largest CAPE) found in the atmosphere from the surface up to 350 hPa. For all the model levels in the lowest 60 hPa, 30-hPa mixed-layer parameters are used for the calculation.', fontsize=10)

    # Display the plot
    time2 = str(i+1)
    base_filename='ecmwf_cape_'
    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)

ecmwf_cape_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,

ecmwf_cape_2.jpeg
ecmwf_cape_3.jpeg
ecmwf_cape_4.jpeg
ecmwf_cape_5.jpeg
ecmwf_cape_6.jpeg
ecmwf_cape_7.jpeg
ecmwf_cape_8.jpeg
ecmwf_cape_9.jpeg
ecmwf_cape_10.jpeg
ecmwf_cape_11.jpeg
ecmwf_cape_12.jpeg
ecmwf_cape_13.jpeg
ecmwf_cape_14.jpeg
ecmwf_cape_15.jpeg
ecmwf_cape_16.jpeg
ecmwf_cape_17.jpeg
ecmwf_cape_18.jpeg
ecmwf_cape_19.jpeg
ecmwf_cape_20.jpeg
ecmwf_cape_21.jpeg
ecmwf_cape_22.jpeg
ecmwf_cape_23.jpeg
ecmwf_cape_24.jpeg
ecmwf_cape_25.jpeg
ecmwf_cape_26.jpeg
ecmwf_cape_27.jpeg
ecmwf_cape_28.jpeg
ecmwf_cape_29.jpeg
ecmwf_cape_30.jpeg
ecmwf_cape_31.jpeg
ecmwf_cape_32.jpeg
ecmwf_cape_33.jpeg
ecmwf_cape_34.jpeg
ecmwf_cape_35.jpeg
ecmwf_cape_36.jpeg
ecmwf_cape_37.jpeg
ecmwf_cape_38.jpeg
ecmwf_cape_39.jpeg
ecmwf_cape_40.jpeg
ecmwf_cape_41.jpeg
ecmwf_cape_42.jpeg
ecmwf_cape_43.jpeg
ecmwf_cape_44.jpeg
ecmwf_cape_45.jpeg
ecmwf_cape_46.jpeg
ecmwf_cape_47.jpeg
ecmwf_cape_48.jpeg
In [6]:
for i in range(48):

    fig = plt.figure(1, figsize=(11, 9))
    ax = fig.add_subplot(111, projection=ccrs.Mercator())
    plot_background(ax)
    

    clevs_li = np.arange(0.01,0.51,0.02)
    # cmap = plt.get_cmap('gist_ncar')
    # newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))

    cmap1 = ListedColormap(colors)
    cmap2 = ListedColormap(colors2)
    newcmap = ListedColormap(cmap1(np.linspace(0, 1, 19)))

    time = data_sfc['mucape'].metpy.time
    data_crs = data_sfc['cin'].metpy.cartopy_crs
    
    ccf2 = ax.contourf(data_sfc.longitude, data_sfc.latitude, li.metpy.loc[{'time': time[i]}], 
                      clevs_li, cmap='gnuplot_r', extend='max', transform=data_crs)
    #axlist[0].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    ax.set_title('Hourly Averaged Total Lightning Density', fontsize=14)
    cb2= fig.colorbar(ccf2, ax=ax, orientation='vertical',
                          shrink=0.825, pad=0)
    cb2.set_label('km⁻² h⁻¹', size='x-large')
    
    # Set figure title
    plt.gcf().text(0.125, 0.92, 'Model: ECMWF IFS 0.1°', fontsize=16)
    #plt.gcf().text(0.4, 0.94, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC'), fontsize=20)
    plt.gcf().text(0.125, 0.89, timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+data_sfc['time'][i].dt.strftime('Valid: %d.%m.%Y %H:%M UTC').values, fontsize=16)
    #plt.gcf().text(0.13, 0.14, 'Note: Updraft helicity is vertically averaged between 1500 and 6000 m AMSL.', fontsize=10)

    # Display the plot
    time2 = str(i+1)
    base_filename='ecmwf_ldensity_'
    suffix='.jpeg'
    latest='latest'
    my_file = base_filename+time2+suffix
    print(my_file)
    plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=100)
    plt.close(fig)

ecmwf_ldensity_1.jpeg
ecmwf_ldensity_2.jpeg
ecmwf_ldensity_3.jpeg
ecmwf_ldensity_4.jpeg
ecmwf_ldensity_5.jpeg
ecmwf_ldensity_6.jpeg
ecmwf_ldensity_7.jpeg
ecmwf_ldensity_8.jpeg
ecmwf_ldensity_9.jpeg
ecmwf_ldensity_10.jpeg
ecmwf_ldensity_11.jpeg
ecmwf_ldensity_12.jpeg
ecmwf_ldensity_13.jpeg
ecmwf_ldensity_14.jpeg
ecmwf_ldensity_15.jpeg
ecmwf_ldensity_16.jpeg
ecmwf_ldensity_17.jpeg
ecmwf_ldensity_18.jpeg
ecmwf_ldensity_19.jpeg
ecmwf_ldensity_20.jpeg
ecmwf_ldensity_21.jpeg
ecmwf_ldensity_22.jpeg
ecmwf_ldensity_23.jpeg
ecmwf_ldensity_24.jpeg
ecmwf_ldensity_25.jpeg
ecmwf_ldensity_26.jpeg
ecmwf_ldensity_27.jpeg
ecmwf_ldensity_28.jpeg
ecmwf_ldensity_29.jpeg
ecmwf_ldensity_30.jpeg
ecmwf_ldensity_31.jpeg
ecmwf_ldensity_32.jpeg
ecmwf_ldensity_33.jpeg
ecmwf_ldensity_34.jpeg
ecmwf_ldensity_35.jpeg
ecmwf_ldensity_36.jpeg
ecmwf_ldensity_37.jpeg
ecmwf_ldensity_38.jpeg
ecmwf_ldensity_39.jpeg
ecmwf_ldensity_40.jpeg
ecmwf_ldensity_41.jpeg
ecmwf_ldensity_42.jpeg
ecmwf_ldensity_43.jpeg
ecmwf_ldensity_44.jpeg
ecmwf_ldensity_45.jpeg
ecmwf_ldensity_46.jpeg
ecmwf_ldensity_47.jpeg
ecmwf_ldensity_48.jpeg
In [ ]: