A Python function for plotting seasonal variation of an output with a depth profile

[Benjamin-Maglio]

Comparing simulation outputs that contain layers of varying depth to each other and against empirical data is required. This function is for doing just that.

We assume you have outputs read as pandas.Dataframe using load_trsc_dataframe() from output_utils.py, and while it will be possible to run this in the terminal, this has been designed for use in Jupyter Notebook.

def seasonal_profile(VAR, depth, thickness, time_range, months, z):
    '''
    VAR : variable dataframe (i.e. TLAYER - temperature by layer)
    depth : associated LAYERDEPTH dataframe
    thickness : associated LAYERDZ dataframe
    time_range : time period to be calculated over -  e.g. ['2011-01-01','2021-01-01']
    months : list of strings - months included in calculation - ['Jan', 'Feb', 'Dec'] (e.g winter season)
    z : array for depth required for interpolation to the same depth profile 
         (typically, np.linspace(min(depth), max(depth), resolution)
    '''
    month_range = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']

    startyr, endyr = time_range[0], time_range[1]

    #Limiting by time range    
    range_series = VAR[startyr:endyr]
    LD = depth[startyr:endyr]
    LZ = thickness[startyr:endyr]

    #Averaging by layer for months in month_range
    range_series = range_series[range_series.index.month.isin([i+1 for i, e in enumerate(month_range) if e in months])]
    LD = (LD[LD.index.month.isin([i+1 for i, e in enumerate(month_range) if e in months])]).mean()
    LZ = (LZ[LZ.index.month.isin([i+1 for i, e in enumerate(month_range) if e in months])]).mean()

    #Creating lists for interpolation
    mean = range_series.mean().values.tolist()    
    maxi = range_series.max().values.tolist()
    mini = range_series.min().values.tolist()
    std =  range_series.std().values.tolist()

    #Creating list with z position at the center of each layer
    zp = (LD + LZ/2).values.tolist()
    
    #Indexes of depths with fill values
    indexes = [i for i,v in enumerate(zp) if v < 0]

    #Removing indexes with fill values
    for index in sorted(indexes, reverse=True):
        del zp[index]
        del mean[index]
        del mini[index]
        del maxi[index]
        del std[index]

    #Interpolating seasonal mean, min, max, and standard deviation
    interp_mean=np.interp(z,zp,mean)
    interp_mini=np.interp(z,zp,mini)
    interp_maxi=np.interp(z,zp,maxi)
    interp_std=np.interp(z,zp,std)
    
    return z, interp_mean, interp_std, interp_mini, interp_maxi

You may generate a trumpet diagram by plotting mean annual (months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']) values:

import matplotlib.pyplot as plt
import numpy as np

fig, ax = plt.subplots()

z, mean, std, mn, mx = seasonal_profile(
                                        TLAYER, LAYERDEPTH, LAYERDZ, ['2011-01-01', '2012-12-01'], 
                                        ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'],
                                        np.linspace(min(LAYERDEPTH), max(LAYERDEPTH), 100)
                                        )


ax.plot(mean, z, label="annual mean")
ax.fill_betweenx(z, mn, mx, alpha=0.2, label="annual range")
ax.plot(np.zeros(len(z)), z, 'k--', alpha=0.5) #Zero degree line

plt.legend(loc="lower right", fontsize=12)

ax.xaxis.tick_top()
ax.set_xlabel(f"Temperature [$^\circ$C]", fontsize=14)
ax.xaxis.set_label_position('top') 

ax.set_ylabel(f"Depth [m]", fontsize=14)
ax.invert_yaxis()

plt.show()

[tobeycarman]

Wow, that's an awesome example!