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

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_sm.nc').sel(longitude=slice(3,9,1),latitude=slice(51,49,1))
data_effis = xr.open_dataset('ec_effis.nc').sel(longitude=slice(3,9,1),latitude=slice(51,49,1))

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

<xarray.Dataset>
Dimensions:    (longitude: 61, latitude: 21, time: 6)
Coordinates:
  * longitude  (longitude) float32 3.0 3.1 3.2 3.3 3.4 ... 8.6 8.7 8.8 8.9 9.0
  * latitude   (latitude) float32 51.0 50.9 50.8 50.7 ... 49.3 49.2 49.1 49.0
  * time       (time) datetime64[ns] 2023-07-12T12:00:00 ... 2023-07-17T12:00:00
Data variables:
    fwinx      (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2023-07-12 06:56:19 GMT by grib_to_netcdf-2.30.2: grib_to_n...

# 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['swvl1'].metpy.coordinates('x', 'y')

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

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

sm =  data_sfc['swvl1']
fw =  data_effis['fwinx']

Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable

2023-07-12 00:00:00

def plot_background(ax):
    ax.set_extent([5.25, 7, 49.2, 50.4])
    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)

for i in range(6):
    
    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_sm = np.arange(0,100,5)
    clevs_fw = [0,5.2,11.2,21.3,38,50,70,100]
    # 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.4, 0.9, 20)))
    bounds = [0,5.2,11.2,21.3,38,50,70,100]
    norm = BoundaryNorm(bounds, newcmap1.N)

    time = data_sfc['swvl1'].metpy.time
    data_crs = data_sfc['swvl1'].metpy.cartopy_crs
        
    # Upper left plot
    cf1 = axlist[0].contourf(data_sfc.longitude, data_sfc.latitude, fw.metpy.loc[{'time': time[i]}], 
                       clevs_fw, cmap=newcmap1, transform=data_crs)
    #ccf1= axlist[0].contour(data_sfc.longitude, data_sfc.latitude, fw.metpy.loc[{'time': time[i]}],
                      #clevs_fw, colors='grey', transform=data_crs)
    #axlist[0].clabel(ccf1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[0].set_title('Forest Fire Weather Index', fontsize=18)
    cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical',
                              shrink=0.999, ticks=(2.6,8.2,16.25,29.6,44,60,85), fraction=0.1, pad=0)
    cb1.set_label('', size='x-large')
    cb1.ax.set_yticklabels(['Very Low','Low', 'Moderate', 'High', 'Very High', 'Extreme', 'Very Extreme'])
    

    cf2 = axlist[1].contourf(data_sfc.longitude, data_sfc.latitude, sm.metpy.loc[{'time': time[i]}]*100, 
                      clevs_sm, cmap='BrBG', 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('Volumetric Soil Water Content 0-7 cm', fontsize=18)
    cb2= fig.colorbar(cf2, ax=axlist[1], orientation='vertical', 
                              shrink=0.999, fraction=0.1, pad=0)
    cb2.set_label('%', size='x-large')


    # Set figure title
    plt.gcf().text(0.125, 0.95, '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.05, 'Fire Danger Classes: Very Low (< 5.2) | Low (5.2 - 11.2) | Moderate (11.2 - 21.3) | High (21.3 - 38.0) | Very High (38.0 - 50.0) | Extreme (50.0 - 70.0) |  Very Extreme (> 70.0) ', fontsize=14)

    # Display the plot
    time2 = str(i+1)
    base_filename='ecmwf_effis_'
    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_effis_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_effis_2.jpeg
ecmwf_effis_3.jpeg
ecmwf_effis_4.jpeg
ecmwf_effis_5.jpeg
ecmwf_effis_6.jpeg

for i in range(6): fig = plt.figure(1, figsize=(11, 9)) ax = fig.add_subplot(111, projection=ccrs.Mercator()) plot_background(ax) clevs_sm = np.arange(0,100,5) # cmap = plt.get_cmap('gist_ncar') # newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30))) time = data_sfc['swvl1'].metpy.time data_crs = data_sfc['swvl1'].metpy.cartopy_crs ccf2 = ax.contourf(data_sfc.longitude, data_sfc.latitude, sm.metpy.loc[{'time': time[i]}]*100, clevs_sm, cmap='BrBG', extend="max", transform=data_crs) #axlist[0].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True) ax.set_title('Volumetric Soil Water Content 0-7 cm', fontsize=16) cb2= fig.colorbar(ccf2, ax=ax, orientation='vertical', shrink=0.825, pad=0) cb2.set_label('%', size='x-large') # Set figure title plt.gcf().text(0.125, 0.96, '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.93, 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_sm_' 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)