<jupyter_start><jupyter_text>Using Copernicus data to investigate Marine Heatwaves     Version: 2.0    Date:    29/06/2020    Author:  <NAME> (EUMETSAT) and <NAME> (Innoflair, Plymouth Marine Laboratory)    Credit:  This code was developed for EUMETSAT under contracts for the Copernicus              programme.    License: This code is offered as open source and free-to-use in the public domain,              with no warranty, under the MIT license associated with this code repository. **What is this notebook for?**This notebook will download EUMETSAT Sentinel-3 SLSTR data for composite plotting, as well as CMEMS time series data of SST, from both reprocessed and NRT data stream. The notebook covers the application of some simple plotting techniques and application of basic analysis to investigate both historical and current potential for the occurrence of marine heat waves.**What specific tools does this notebook use?**Beyond general Python modules, this notebook imports some functions for managing the harmonised data access api (harmonised_data_access_api_tools.py) which can be found in the wekeo-hda folder on JupyterLab, and additional libraries for marine heatwave analysis provided openly and based on the work of  Hobday et al..**What are marine heatwaves and how can Copernicus data be used to observe them?**Like heatwaves on land, marine heatwaves are extended periods of higher than normal temperatures. They have occurred in different areas around the world and can be caused by different oceanographic driving forces. You can find out all about them here. The variable drivers, and historical context for defining heatwaves regionally, mean that we need the ability to measure the range of ocean temperatures that occur all over the world at any given time, and also a long-term base-line understanding of what ‘normal’ temperatures look like.Sea surface temperature measurements from satellites can support the identification of marine heatwaves, both through the daily measurements they make, and their contributions to long-term data records. The Sentinel-3 satellites are the Copernicus programme's contribution to climate-scale monitoring of sea surface temperatures, with the Sea and Land Surface Temperature Radiometer (SLSTR) on Sentinel-3A now able to function as a reference sensor, when combining multiple sea surface temperature data sets, such as those available from the Copernicus Marine Service. *** Let's get started! Python is divided into a series of modules that each contain a series of methods for specific tasks. The box below imports all of the moduls we need to complete our plotting task<jupyter_code># standard tools
import os, sys, json
from zipfile import ZipFile
import shutil
import sys
import datetime
import numpy as np
from IPython.core.display import HTML
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib import gridspec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import xmltodict
import glob
import warnings
warnings.filterwarnings("ignore")

# HDA API tools
sys.path.append(os.path.join(os.path.dirname(os.path.dirname(os.getcwd())),'wekeo-hda'))
import hda_api_functions as hapi # this is the PYTHON version
        
# specific tools (which can be found here ../../tools/)
sys.path.append(os.path.join(os.getcwd(),'tools'))
sys.path.append(os.path.join(os.getcwd(),'tools','mhw_master'))
import image_tools as img
import marineHeatWaves as mhw<jupyter_output><empty_output><jupyter_text>WEkEO provides access to a huge number of datasets through its **'harmonised-data-access'** API. This allows us to query the full data catalogue and download data quickly and directly onto the Jupyter Lab. You can search for what data is available hereIn order to use the HDA-API we need to provide some authentication credentials, which comes in the form of an API key and API token. In this notebook we have provided functions so you can retrieve the API key and token you need directly. You can find out more about this process in the notebook on HDA access (wekeo_harmonized_data_access_api.ipynb) that can be found in the **wekeo-hda** folder on your Jupyterlab.We will also define a few other parameters including where to download the data to, and if we want the HDA-API functions to be verbose. **Lastly, we will tell the notebook where to find the query we will use to find the data.** These 'JSON' queries are what we use to ask WEkEO for data. They have a very specific form, but allow us quite fine grained control over what data to get. You can find the ones that we will use here: **JSON_templates/mhw/..**<jupyter_code># your WEkEO API username and password (needs to be in '  ')
user_name = 'Hek17'
password = '<PASSWORD>!'

# Generate an API key
api_key = hapi.generate_api_key(user_name, password)
display(HTML('Your API key is: <b>'+api_key+'</b>'))
# where the data should be downloaded to:
download_dir_path = os.path.join('C:/','Users','Evers','Desktop','products')
# where we can find our data query form:
JSON_query_dir = os.path.join(os.getcwd(),'JSON_templates','mhw')
# HDA-API loud and noisy?
verbose = False

# make the output directory if required
if not os.path.exists(download_dir_path):
    os.makedirs(download_dir_path)<jupyter_output><empty_output><jupyter_text>Now we are ready to get our data.<jupyter_code># SLSTR L2 SST KEY
dataset_id = "EO:EUM:DAT:SENTINEL-3:SL_2_WST___"<jupyter_output><empty_output><jupyter_text>We use our dataset ID to load the correct JSON query file from ../JSON_templates/mhw/You can edit this query if you want to get different data, but be aware of asking for too much data - you could be here a while and might run out of space to use this data in the JupyterLab. The box below gets the correct query file.<jupyter_code># find query file
JSON_query_file = os.path.join(JSON_query_dir,dataset_id.replace(':','_')+".json")
if not os.path.exists(JSON_query_file):
    print('Query file ' + JSON_query_file + ' does not exist')
else:
    print('Found JSON query file for '+dataset_id)<jupyter_output><empty_output><jupyter_text>Now we have a query, we need to launch it to WEkEO to get our data. The box below takes care of this through the following steps:    1. initialise our HDA-API    2. get an access token for our data    3. accepts the WEkEO terms and conditions    4. loads our JSON query into memory    5. launches our search    6. waits for our search results    7. gets our result list    8. downloads our dataThis is quite a complex process, so much of the functionality has been buried 'behind the scenes'. If you want more information, you can check out the Harmonised Data Access API  tutorials. The code below will report some information as it runs. At the end, it should tell you that one product has been downloaded.<jupyter_code>HAPI_dict = hapi.init(dataset_id, api_key, download_dir_path)
HAPI_dict = hapi.get_access_token(HAPI_dict)
HAPI_dict = hapi.acceptTandC(HAPI_dict)

# load the query
with open(JSON_query_file, 'r') as f:
    query = json.load(f)

# launch job
print('Launching job...')
HAPI_dict = hapi.get_job_id(HAPI_dict, query)

# check results
print('Getting results...')
HAPI_dict = hapi.get_results_list(HAPI_dict)
HAPI_dict = hapi.get_order_ids(HAPI_dict)

# download data
print('Downloading data...')
HAPI_dict = hapi.download_data(HAPI_dict, file_extension='.zip')<jupyter_output><empty_output><jupyter_text>Sentinel data is usually distributed as a zip file, which contains the SAFE format data within. To use this, we must unzip the file. The bow below handles this.<jupyter_code># unzip file
for filename in HAPI_dict['filenames']:
    if os.path.splitext(filename)[-1] == '.zip':
        print('Unzipping file')
        with ZipFile(filename, 'r') as zipObj:
            # Extract all the contents of zip file in current directory
            zipObj.extractall(os.path.dirname(filename))<jupyter_output><empty_output><jupyter_text>***  Plotting SLSTR dataWe plot our SLSTR scenes using a function that manages data ingestion, flagging, bias correction and makes some map embelishements (e.g. adds dotted lines to the scene edges, so we can tell where the boundaries are). We call this function for each image in the boxes further down.<jupyter_code>def make_SLSTR_composite_plot(nc_files, plot_extents=None, fsz=20,\
                              land_resolution='50m', vmin=10, vmax=28,\
                              xsize=20, ysize=16, dpi=150, QMASK=4,\
                              cmap=plt.cm.RdYlBu_r):
    '''
     Plots SLSTR images from n input files, will overay first to last.
    '''
    # get cartopy land layers
    land_poly = cfeature.NaturalEarthFeature('physical', 'land', land_resolution,
                                             edgecolor='k',
                                             facecolor=cfeature.COLORS['land'])
    
    # setup figure
    fig = plt.figure(figsize=(xsize, ysize), dpi=dpi)
    plt.rc('font', size=fsz)
    
    # setup axes
    gs = gridspec.GridSpec(3, 1, height_ratios=[20,0.5,1])
    gs.update(wspace=0.01, hspace=0.01)

    # setup plot 1
    axes_m = plt.subplot(gs[0,0], projection=ccrs.PlateCarree())
    axes_m.background_patch.set_facecolor('0.5')


    if plot_extents:
        axes_m.set_extent(plot_extents, ccrs.PlateCarree())

    # read data and plot background
    for nc_file in nc_files:
        nc_fid = xr.open_dataset(nc_file)
        LON = nc_fid.lon.data
        LAT = nc_fid.lat.data
        
        x1 = 0 ; x2 = -1; y1 = 0 ; y2 = -1
        # get subset based on plot extents
        if plot_extents:
            x1, x2, y1, y2 = img.subset_image(LAT, LON, plot_extents)

        # read in and subset if required
        LON = LON[x1:x2,y1:y2]
        LAT = LAT[x1:x2,y1:y2]
        SST_raw = np.squeeze(nc_fid.sea_surface_temperature.data[0,x1:x2,y1:y2])
        SST_BIAS = np.squeeze(nc_fid.sses_bias.data[0,x1:x2,y1:y2])
        QUALITY_LEVEL = np.squeeze(nc_fid.quality_level.data[0,x1:x2,y1:y2])
        nc_fid.close()
        
        # correct SST
        SST_C = SST_raw + SST_BIAS - 273.15

        # flag data
        SST_C[QUALITY_LEVEL<=QMASK] = np.nan

        # plot the SST field
        p1 = axes_m.pcolormesh(LON, LAT, SST_C, cmap=cmap,\
                          vmin=vmin, vmax=vmax, zorder=-1)

        # add the plot edges
        xml_file = os.path.dirname(nc_file)+'/xfdumanifest.xml'
        with open(xml_file) as fd:
            doc = xmltodict.parse(fd.read())
            coords = doc['xfdu:XFDU']['metadataSection']['metadataObject'][2]\
                        ['metadataWrap']['xmlData']['sentinel-safe:frameSet']\
                        ['sentinel-safe:footPrint']['gml:posList']

        # split the coords
        lats_lons = coords.split(' ')
        lats = np.asarray(lats_lons[::2]).astype(float)
        lons = np.asarray(lats_lons[1::2]).astype(float)

        if 'S3B' in nc_file:
            plot_col = 'b'
        else:
            plot_col = 'k'

        p2 = axes_m.plot(lons, lats, c=plot_col, linewidth=1,\
                    linestyle='--', zorder=5, alpha=0.5,\
                    transform=ccrs.PlateCarree())

    # add some map embelishments
    axes_m.coastlines(resolution=land_resolution, color='black', linewidth=1)
    axes_m.add_feature(land_poly)
    g1 = axes_m.gridlines(draw_labels = True)
    g1.xlabels_top = False
    g1.ylabels_right = False

    g1.xlabel_style = {'size': fsz, 'color': 'gray'}
    g1.ylabel_style = {'size': fsz, 'color': 'gray'}

    # setup plot2: colorbar
    axes_c = plt.subplot(gs[2,0])
    cbar = plt.colorbar(p1, cax=axes_c, orientation='horizontal')
    cbar.ax.tick_params(labelsize=fsz) 
    cbar.set_label('SLSTR SST composite [$^{o}$C]',fontsize=fsz)
    
    return gs, axes_m, axes_c<jupyter_output><empty_output><jupyter_text>We start by finding all the necessary files (which glob.glob takes care of). Th files are added to a list which is then sent to our large function above.<jupyter_code># verbosity
verbose = False

# figure options
figure_font_size = 20
plot_extents = [-160, -116, 10, 45]
vmin = 10
vmax = 28

# get the files
SLSTR_files = glob.glob(os.path.join(download_dir_path,'S3*WST*201909*','*.nc'))<jupyter_output><empty_output><jupyter_text>And now we pass our list of files to the plotting routine. The plotting routine returns the handles of our axes, so that we can still make some changes once the main plot is done (e.g. add the annotations). Finally, it will save fie figure.<jupyter_code># make the plot: we will call this as a function as in contains a 'for' loop to make the plot
fig, axis, colbar = make_SLSTR_composite_plot(SLSTR_files, plot_extents=plot_extents,\
                                              fsz=figure_font_size, vmin=vmin, vmax=vmax)

# add some embellishments
plt.sca(axis)

label='Sentinel-3A SLSTR SST\n07/09/2019 (22:27 local time)'    
txt = plt.annotate(label, xy=(0.66, 1.01), xycoords='axes fraction',\
               size=figure_font_size/1.25,
               color='k', zorder=100, annotation_clip=False)

label='Sentinel-3A SLSTR SST\n08/09/2019 (00:08 local time)'    
txt = plt.annotate(label, xy=(0.005, 1.01), xycoords='axes fraction',\
               size=figure_font_size/1.25,
               color='k', zorder=100, annotation_clip=False)

label='Sentinel-3B SLSTR SST\n08/09/2019 (23:28 local time)'    
txt = plt.annotate(label, xy=(0.33, 1.01), xycoords='axes fraction',\
               size=figure_font_size/1.25,
               color='b', zorder=100, annotation_clip=False)

plt.savefig('SLSTR_All_SST_California_20190909.png',bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>You can now find the image in the folder where your code is stored in your instance of the JupyterLab. The image is a composite of three images from the SLSTR sensors aboard the Sentinel-3A and B satellites which captured warm temperatures around the eastern Pacific in September 2019. This highlights the increased coverage that can be achieved in both time and space, using multiple sensors, but in order to understand if these temperatures could be an indicator of the beginning of a marine heatwave we need some longer term context... ***  Looking at time series data of sea surface temperature To place the images above in context we'll need to look at a time series of data. For this we can access data from the Copernicus services, where multiple data sources (different satellites etc) are combined to produce data sets that cover longer time periods. In this case we are going to look at the OSTIA reprocessed SST, as well as the near real time (NRT)processing of the same product, so we can look at data right up to the time period of the SLSTR images we looked at previously. Reprocessed and NRT data streams are produced separately and we must be careful when we interpret these data sets together. As time progresses, we understand better how satellite instruments perform, algorithms improve, and data is reprocessed to ensure good quality and consistency between sources. With NRT we do not have all the information we have in hindsight, particularly about instrument characterisation, but this data is vitally important for understanding events that are happening right now. So, for example, you would not use combined NRT and reprocessed data sets for long-term trend analysis, and you would want to consider NRT measurements with a higher degree of uncertainty that you might consider with reprocessed data.  As before, we will construct and submit a query to the WEkEO Harmonised Data Access API to get this data. Here we have supplied the JSON files for the data sets of interest, but you could edit these files to look at your own time frames/regions of interest etc.<jupyter_code># find query file for OSTIA NRT and initialise dictionary
dataset_id = "EO:MO:DAT:SST_GLO_SST_L4_NRT_OBSERVATIONS_010_001"

# makes a date array...NfH
start_dates = []
end_dates = []
for ii in range(2008,2019+1):
    start_dates.append(str(ii)+"-01-01T00:00:00.000Z")
    end_dates.append(str(ii)+"-12-31T00:00:00.000Z")

# find query file
JSON_query_file = os.path.join(JSON_query_dir,dataset_id.replace(':','_')+".json")
if not os.path.exists(JSON_query_file):
    print('Query file ' + JSON_query_file + ' does not exist')
else:
    print('Found JSON query file for '+dataset_id)

# init HAPI
HAPI_dict = hapi.init(dataset_id, api_key, download_dir_path)
HAPI_dict = hapi.get_access_token(HAPI_dict)
HAPI_dict = hapi.acceptTandC(HAPI_dict)

# start loop over dates
for start_date, end_date in zip(start_dates, end_dates):
    print('Running for: ' + start_date + ' to ' + end_date)
    
    date_tag = start_date.split('T')[0] + '_' + end_date.split('T')[0]
    date_string = start_date.split('T')[0].replace('-','') \
                  + '_' + end_date.split('T')[0].replace('-','')

    # load the query
    with open(JSON_query_file, 'r') as f:
        query = f.read()
        query = query.replace("%DATE_START%",start_date)
        query = query.replace("%DATE_END%",end_date)
        query = json.loads(query)

    # launch job
    print('Launching job...')
    HAPI_dict = hapi.get_job_id(HAPI_dict, query)

    # check results
    print('Getting results...')
    HAPI_dict = hapi.get_results_list(HAPI_dict)
    HAPI_dict = hapi.get_order_ids(HAPI_dict)

    # download data
    print('Downloading data...')
    HAPI_dict = hapi.download_data(HAPI_dict, \
                  user_filename='METOFFICE-GLO-SST-L4-NRT_' + date_tag + '.nc')
# find query file for OSTIA reprocessed and initialise dictionary

dataset_id = "EO_MO_DAT_SST_GLO_SST_L4_REP_OBSERVATIONS_010_011"

# makes a date array...NfH
start_dates = []
end_dates = []
for ii in range(1988,2007+1):
    start_dates.append(str(ii)+"-01-01T00:00:00.000Z")
    end_dates.append(str(ii)+"-12-31T00:00:00.000Z")

# find query file
JSON_query_file = os.path.join(JSON_query_dir,dataset_id.replace(':','_')+".json")
if not os.path.exists(JSON_query_file):
    print('Query file ' + JSON_query_file + ' does not exist')
else:
    print('Found JSON query file for '+dataset_id)

# init HAPI
HAPI_dict = hapi.init(dataset_id, api_key, download_dir_path)
HAPI_dict = hapi.get_access_token(HAPI_dict)
HAPI_dict = hapi.acceptTandC(HAPI_dict)

# start loop over dates
for start_date, end_date in zip(start_dates, end_dates):
    print('Running for: ' + start_date + ' to ' + end_date)

    date_tag = start_date.split('T')[0] + '_' + end_date.split('T')[0]
    date_string = start_date.split('T')[0].replace('-','') \
                  + '_' + end_date.split('T')[0].replace('-','')
    
    # load the query
    with open(JSON_query_file, 'r') as f:
        query = f.read()
        query = query.replace("%DATE_START%",start_date)
        query = query.replace("%DATE_END%",end_date)
        query = json.loads(query)

    # launch job
    print('Launching job...')
    HAPI_dict = hapi.get_job_id(HAPI_dict, query)

    # check results
    print('Getting results...')
    HAPI_dict = hapi.get_results_list(HAPI_dict)
    HAPI_dict = hapi.get_order_ids(HAPI_dict)

    # download data
    print('Downloading data...')
    HAPI_dict = hapi.download_data(HAPI_dict, \
                  user_filename='METOFFICE-GLO-SST-L4-RAN_' + date_tag + '.nc')<jupyter_output><empty_output><jupyter_text>***  Plotting time series data to investigate occurrence and potential for marine heatwaves  First we will set up some parameters including times and space we want to work over when plotting.<jupyter_code># subset image: cut a relevant section out of an image for area averaging. If false, whole area will be used.
# THERE IS NO AREA WEIGHTING IN THE AVERAGING HERE! 
# Subset_extents [lon1,lon2,lat1,lat2] describes the section we cut
subset_image = False
subset_extents = [-165.0, -145.0, 10.0, 30.0]

# The data product time is measured in seconds since the following refernce date
Date_ref = datetime.datetime(1981,1,1)

# Select the times we want to use for our spatial anomaly plots [month, day]
month_day_start = [8,1] 
month_day_end = [9,24]

# Plotting font size
fsz = 30<jupyter_output><empty_output><jupyter_text>Then we will find our datasets and concatenate them.<jupyter_code>SST_files = []
my_files = glob.glob(os.path.join(download_dir_path,'METOFFICE-GLO-SST-L4-RAN*'))
for SST_file in sorted(my_files):
    SST_files.append(SST_file)
    
my_files = glob.glob(os.path.join(download_dir_path,'METOFFICE-GLO-SST-L4-NRT*'))
for SST_file in sorted(my_files):
    SST_files.append(SST_file)<jupyter_output><empty_output><jupyter_text>The next cell performs the bulk of the processing work for our plot. We start by loading in coordinates that we can use to subset the data (if required) and make our plots. Subsequently, we loop through the SST products and read in the data, correct Kelvin to Celsius and perform our averaging in either time and space.<jupyter_code>#load the co-ordinate variables we need for subsetting/plotting
ds = xr.open_dataset(SST_files[-1])
lat = ds.lat.data
lon = ds.lon.data
ds.close()

# subset coords if required, getting indices to subset out output SST products
if subset_image:
    ii = np.where((lon >= subset_extents[0]) & (lon <= subset_extents[1]))[0]
    jj = np.where((lat >= subset_extents[2]) & (lat <= subset_extents[3]))[0]
    lon = lon[ii]
    lat = lat[jj]

# initialise lists for output times series variables for MWH    
all_times = []
all_SST = []

# initialise arrays for output SST fields
iter_SST = np.ones([len(SST_files), len(lat), len(lon)])*np.nan

# now we get the area-averaged data
count = -1
for SST_file in SST_files:
    print(SST_file)
    count = count + 1

    # xarray does not read times consistency between RAN and NRT, so we load as integer
    ds = xr.open_dataset(SST_file, decode_times=False)
    this_SST = ds.analysed_sst.data
    times = ds.time.data
    ds.close()

    my_times = []
    for time in times:
        my_times.append(Date_ref + datetime.timedelta(seconds=int(time)))
    times = np.asarray(my_times)
    
    t0 = datetime.datetime(times[0].year, month_day_start[0], month_day_start[1])
    t1 = datetime.datetime(times[0].year, month_day_end[0], month_day_end[1])
    tt = np.where((times >= t0) & (times <= t1))[0]

    if np.nanmean(this_SST) > 100:
        this_SST = this_SST - 273.15

    if subset_image:
        this_SST = this_SST[:,jj[0]:jj[-1],ii[0]:ii[-1]]
    
    time_subset_SST = this_SST[tt,:,:]

    iter_SST[count,:,:] = np.squeeze(np.nanmean(time_subset_SST, axis=0))

    all_times.append(my_times)
    all_SST.append(np.nanmean(np.nanmean(this_SST, axis=1), axis=1))<jupyter_output><empty_output><jupyter_text>A little bit of formatting is necessary for the plots we want to make...<jupyter_code># flatten the SST list
SST_time_series = [item for sublist in all_SST for item in sublist]
SST_time_series = np.asarray(SST_time_series)

# format the dates for the MWH toolkit
Dates_time_series = [item for sublist in all_times for item in sublist]
Dates_time_series_formatted = [datetime.date.toordinal(tt) for tt in Dates_time_series]

# make climatology of time subset region
clim_SST = np.nanmean(iter_SST, axis=0)

# calculate heat waves
mhws, clim = mhw.detect(np.asarray(Dates_time_series_formatted), SST_time_series)

# make matrix
stripe_array = np.ones([20,len(SST_files)])*np.nan

# now make the anomaly
for ii in range(len(SST_files)):
    stripe_array[:, ii] = np.nanmean(np.squeeze(iter_SST[ii,:,:]) - clim_SST)<jupyter_output><empty_output><jupyter_text>The first plot we'll make, shows ‘stripes’ of the average sea surface temperature anomaly for the region around Hawaii. Recently, ‘climate stripes’  have been used by 'citizen scientists' all over the world to show long-term trends in regional temperatures. The plot below shows a stripes-style graphic for the eastern Pacific region around Hawaii. High anomalies are apparent during 2014, in 2015 during the previous marine heatwave, and were associated with reports in 2019.<jupyter_code>fig = plt.figure(figsize=(30,10), dpi =300)
vmax = np.nanmax(abs(stripe_array))

date_ticks = []
for ii in range(len(SST_files)):
    date_ticks.append(str(1988+ii))
    
plt.pcolormesh(stripe_array,vmin=vmax*-1,vmax=vmax,cmap=plt.cm.RdBu_r)
plt.xticks(np.arange(len(SST_files))+0.5,date_ticks, rotation='90')
plt.yticks([],[])
cbar = plt.colorbar()
cbar.set_label('SST anomaly [$^{o}$C] (1$^{st}$ Aug - 24$^{th}$ Sep)',fontsize=fsz/1.25)
plt.savefig('Stripes.png')<jupyter_output><empty_output><jupyter_text>Finally, we will look a little more quantitatively at historical marine heatwaves in this region, and see how the situation in 2019 compared (with the previously mentioned assumptions about the accuracy of NRT data for longer term analysis). The plot below is generated using a a toolkit very kindly provided as open source code by Hobday et al (you can find out more information on the toolkit and the science behind it here and here)<jupyter_code># plot MWHs
ev = np.argmax(mhws['intensity_max']) # Find largest event

fig = plt.figure(figsize=(35,15), dpi = 300)
plt.rc('font', size=fsz)

# Find indices for all n MHWs before and after event of interest and shade accordingly
n=10
for ev0 in np.arange(ev-n, ev+n, 1):
    t1 = np.where(Dates_time_series_formatted==mhws['time_start'][ev0])[0][0]
    t2 = np.where(Dates_time_series_formatted==mhws['time_end'][ev0])[0][0]
    p1 = plt.fill_between(Dates_time_series[t1:t2+1], SST_time_series[t1:t2+1],\
                          clim['thresh'][t1:t2+1], color=(1,0.85,0.85))
    
# Find indices for MHW of interest and shade accordingly
t1 = np.where(Dates_time_series_formatted==mhws['time_start'][-1])[0][0]
t2 = np.where(Dates_time_series_formatted==mhws['time_end'][-1])[0][0]
p2 = plt.fill_between(Dates_time_series[t1:t2+1], SST_time_series[t1:t2+1],\
                      clim['thresh'][t1:t2+1], color='r')

# Plot SST, seasonal cycle, threshold, shade MHWs with main event in red
p3, = plt.plot(Dates_time_series, SST_time_series, 'k-', linewidth=2)
p4, = plt.plot(Dates_time_series, clim['thresh'], 'b--', linewidth=2)
p5, = plt.plot(Dates_time_series, clim['seas'], 'b-', linewidth=2)

xmin = datetime.datetime(2014,1,1).toordinal() - datetime.datetime(1,1,1).toordinal()
xmax = datetime.datetime(2019,12,31).toordinal() - datetime.datetime(1,1,1).toordinal()
plt.xlim(xmin,xmax)

plt.ylim(clim['seas'].min() - 0.3, clim['seas'].max() + mhws['intensity_max'][ev] + 0.2)
plt.ylabel('SST [$^\circ$C]')
plt.annotate('Plotting script credit:\nhttps://github.com/ecjoliver/marineHeatWaves',\
             (0.005, 0.935), xycoords='axes fraction', color='0.5', fontsize=fsz/1.25)

leg = plt.legend([p3, p5, p4, p2, p1],\
                 ['OSTIA NRT SST','SST Seas. clim.','SST Seas. thresh.','Current heatwave','Past heatwave'],\
                 bbox_to_anchor=(1.0, -0.05), ncol=5)

leg.get_frame().set_linewidth(0.0)
plt.savefig('MHW.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10312418737694565, 'fraction_numerical': 0.0203573684015008, 'mean_word_length': 4.027077497665733, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 74, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '166399', 'license': 'MIT', 'n_tokens_mistral': 9147, 'n_tokens_neox': 8720, 'n_words': 2775, 'path': 'ocean/OceanCaseStudies/.ipynb_checkpoints/Marine_heatwaves_S3_CMEMS-WITH-REP-checkpoint.ipynb', 'repo_name': 'stewartchrisecmwf/wekeo-jupyter-lab', 'size': '26934'}