plot_utils.py

Go to the documentation of this file.

from .metrics import slope, sspb, mdsa, rmsle
from .meta import get_sensor_label
from .utils import closest_wavelength, ignore_warnings
from collections import defaultdict as dd 
from pathlib import Path 
import numpy as np 
 
 
def add_identity(ax, *line_args, **line_kwargs):
    ''' 
    Add 1 to 1 diagonal line to a plot.
    https://stackoverflow.com/questions/22104256/does-matplotlib-have-a-function-for-drawing-diagonal-lines-in-axis-coordinates
    
    Usage: add_identity(plt.gca(), color='k', ls='--')
    '''
    line_kwargs['label'] = line_kwargs.get('label', '_nolegend_')
    identity, = ax.plot([], [], *line_args, **line_kwargs)
    
    def callback(axes):
        low_x, high_x = ax.get_xlim()
        low_y, high_y = ax.get_ylim()
        lo = max(low_x,  low_y)
        hi = min(high_x, high_y)
        identity.set_data([lo, hi], [lo, hi])
 
    callback(ax)
    ax.callbacks.connect('xlim_changed', callback)
    ax.callbacks.connect('ylim_changed', callback)
 
    ann_kwargs = {
        'transform'  : ax.transAxes,
        'textcoords' : 'offset points', 
        'xycoords'   : 'axes fraction', 
        'fontname'   : 'monospace', 
        'xytext'     : (0,0), 
        'zorder'     : 25,  
        'va'         : 'top', 
        'ha'         : 'left', 
    }
    ax.annotate(r'$\mathbf{1:1}$', xy=(0.87,0.99), size=11, **ann_kwargs)
 
 
 
def _create_metric(metric, y_true, y_est, longest=None, label=None):
    ''' Create a position-aligned string which shows the performance via a single metric '''
    if label == None:   label = metric.__name__.replace('SSPB', '\\beta').replace('MdSA', '\\varepsilon\\thinspace')#.replace('Slope','S\\thinspace')
    # if label == None:   label = metric.__name__.replace('SSPB', 'Bias').replace('MdSA', 'Error')
    if longest == None: longest = len(label)
 
    ispct = metric.__qualname__ in ['mape', 'sspb', 'mdsa'] # metrics which are percentages
    diff  = longest-len(label.replace('^',''))
    space = r''.join([r'\ ']*diff + [r'\thinspace']*diff)
    prec  = (1 if abs(metric(y_true, y_est)) < 100 and metric.__name__ not in ['N'] else 0) if ispct or metric.__name__ in ['N'] else 3
    # prec  = 1 if abs(metric(y_true, y_est)) < 100 else 0
    stat  = f'{metric(y_true, y_est):.{prec}f}'
    perc  = r'$\small{\mathsf{\%}}$' if ispct else ''
    return rf'$\mathtt{{{label}}}{space}:$ {stat}{perc}'
 
def _create_stats(y_true, y_est, metrics, title=None):
    ''' Create stat box strings for all metrics, assuming there is only a single target feature '''
    longest = max([len(metric.__name__.replace('SSPB', 'Bias').replace('MdSA', 'Error').replace('^','')) for metric in metrics])
    statbox = [_create_metric(m, y_true, y_est, longest=longest) for m in metrics]
    
    if title is not None:
        statbox = [rf'$\mathbf{{\underline{{{title}}}}}$'] + statbox
    return statbox 
 
def _create_multi_feature_stats(y_true, y_est, metrics, labels=None):
    ''' Create stat box strings for a single metric, assuming there are multiple target features '''
    if labels == None: 
        labels = [f'Feature {i}' for i in range(y_true.shape[1])]
    assert(len(labels) == y_true.shape[1] == y_est.shape[1]), f'Number of labels does not match number of features: {labels} - {y_true.shape}'
    
    title   = metrics[0].__name__.replace('SSPB', 'Bias').replace('MdSA', 'Error')
    longest = max([len(label.replace('^','')) for label in labels])
    statbox = [_create_metric(metrics[0], y1, y2, longest=longest, label=lbl) for y1, y2, lbl in zip(y_true.T, y_est.T, labels)]
    statbox = [rf'$\mathbf{{\underline{{{title}}}}}$'] + statbox
    return statbox 
 
def add_stats_box(ax, y_true, y_est, metrics=[mdsa, sspb, slope], bottom_right=False, bottom=False, right=False, x=0.025, y=0.97, fontsize=16, label=None):
    ''' Add a text box containing a variety of performance statistics, to the given axis '''
    import matplotlib.pyplot as plt
    plt.rc('text', usetex=True)
    plt.rcParams['mathtext.default']='regular'
 
    create_box = _create_stats if len(y_true.shape) == 1 or y_true.shape[1] == 1 else _create_multi_feature_stats
    stats_box  = '\n'.join( create_box(y_true, y_est, metrics, label) )
    ann_kwargs = {
        'transform'  : ax.transAxes,
        'textcoords' : 'offset points', 
        'xycoords'   : 'axes fraction', 
        'fontname'   : 'monospace', 
        'xytext'     : (0,0), 
        'zorder'     : 25,  
        'va'         : 'top', 
        'ha'         : 'left', 
        'bbox'       : {
            'facecolor' : 'white',
            'edgecolor' : 'black', 
            'alpha'     : 0.7,
        }
    }
 
    ann = ax.annotate(stats_box, xy=(x,y), size=fontsize, **ann_kwargs)
 
    bottom |= bottom_right
    right  |= bottom_right
 
    # Switch location to (approximately) the bottom right corner
    if bottom or right or bottom_right:
        plt.gcf().canvas.draw()
        bbox_orig = ann.get_tightbbox(plt.gcf().canvas.renderer).transformed(ax.transAxes.inverted())
 
        new_x = bbox_orig.x0
        new_y = bbox_orig.y1
        if bottom:
            new_y = bbox_orig.y1 - bbox_orig.y0 + (1 - y)
            ann.set_y(new_y)
            new_y += 0.06
        if right:
            new_x = 1 - (bbox_orig.x1 - bbox_orig.x0) + x
            ann.set_x(new_x)
            new_x -= 0.04
        ann.xy = (new_x, new_y)
    return ann 
    
 
def draw_map(*lonlats, scale=0.2, world=False, us=True, eu=False, labels=[], ax=None, gray=False, res='i', **scatter_kws):
    ''' Helper function to plot locations on a global map '''
    import matplotlib.pyplot as plt
    from matplotlib.transforms import Bbox
    from mpl_toolkits.axes_grid1.inset_locator import TransformedBbox, BboxPatch, BboxConnector 
    from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
    from mpl_toolkits.basemap import Basemap
    from itertools import chain
 
    PLOT_WIDTH  = 8
    PLOT_HEIGHT = 6
 
    WORLD_MAP = {'cyl': [-90, 85, -180, 180]}
    US_MAP    = {
        'cyl' : [24, 49, -126, -65],
        'lcc' : [23, 48, -121, -64],
    }
    EU_MAP    = {
        'cyl' : [34, 65, -12, 40],
        'lcc' : [30.5, 64, -10, 40],
    }
 
    def mark_inset(ax, ax2, m, m2, MAP, loc1=(1, 2), loc2=(3, 4), **kwargs):
        """
        https://stackoverflow.com/questions/41610834/basemap-projection-geos-controlling-mark-inset-location
        Patched mark_inset to work with Basemap.
        Reason: Basemap converts Geographic (lon/lat) to Map Projection (x/y) coordinates
 
        Additionally: set connector locations separately for both axes:
            loc1 & loc2: tuple defining start and end-locations of connector 1 & 2
        """
        axzoom_geoLims = (MAP['cyl'][2:], MAP['cyl'][:2]) 
        rect = TransformedBbox(Bbox(np.array(m(*axzoom_geoLims)).T), ax.transData)
        pp   = BboxPatch(rect, fill=False, **kwargs)
        ax.add_patch(pp)
        p1 = BboxConnector(ax2.bbox, rect, loc1=loc1[0], loc2=loc1[1], **kwargs)
        ax2.add_patch(p1)
        p1.set_clip_on(False)
        p2 = BboxConnector(ax2.bbox, rect, loc1=loc2[0], loc2=loc2[1], **kwargs)
        ax2.add_patch(p2)
        p2.set_clip_on(False)
        return pp, p1, p2
 
 
    if world:
        MAP    = WORLD_MAP
        kwargs = {'projection': 'cyl', 'resolution': res}
    elif us:
        MAP    = US_MAP
        kwargs = {'projection': 'lcc', 'lat_0':30, 'lon_0':-98, 'resolution': res}#, 'epsg':4269}
    elif eu:
        MAP    = EU_MAP
        kwargs = {'projection': 'lcc', 'lat_0':48, 'lon_0':27, 'resolution': res}
    else:
        raise Exception('Must plot world, US, or EU')
 
    kwargs.update(dict(zip(['llcrnrlat', 'urcrnrlat', 'llcrnrlon', 'urcrnrlon'], MAP['lcc' if 'lcc' in MAP else 'cyl'])))
    if ax is None: f = plt.figure(figsize=(PLOT_WIDTH, PLOT_HEIGHT), edgecolor='w')
    m  = Basemap(ax=ax, **kwargs)
    ax = m.ax if m.ax is not None else plt.gca()
 
    if not world:
        m.readshapefile(Path(__file__).parent.joinpath('map_files', 'st99_d00').as_posix(), name='states', drawbounds=True, color='k', linewidth=0.5, zorder=11)
        m.fillcontinents(color=(0,0,0,0), lake_color='#9abee0', zorder=9)
        if not gray:
            m.drawrivers(linewidth=0.2, color='blue', zorder=9)
        m.drawcountries(color='k', linewidth=0.5)
    else:
        m.drawcountries(color='w')
    # m.bluemarble()
    if not gray:
        if us or eu: m.shadedrelief(scale=0.3 if world else 1)
        else:
            # m.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 2000, verbose= True)
            m.arcgisimage(service='World_Imagery', xpixels = 2000, verbose= True)
    else:
        pass
    # lats = m.drawparallels(np.linspace(MAP[0], MAP[1], 13))
    # lons = m.drawmeridians(np.linspace(MAP[2], MAP[3], 13))
 
    # lat_lines = chain(*(tup[1][0] for tup in lats.items()))
    # lon_lines = chain(*(tup[1][0] for tup in lons.items()))
    # all_lines = chain(lat_lines, lon_lines)
    
    # for line in all_lines:
    #   line.set(linestyle='-', alpha=0.0, color='w')
 
    if labels:
        colors = ['aqua', 'orangered',  'xkcd:tangerine', 'xkcd:fresh green', 'xkcd:clay', 'magenta', 'xkcd:sky blue', 'xkcd:greyish blue', 'xkcd:goldenrod', ]
        markers = ['o', '^', 's', '*',  'v', 'X', '.', 'x',]
        mod_cr = False
        assert(len(labels) == len(lonlats)), [len(labels), len(lonlats)]
        for i, (label, lonlat) in enumerate(zip(labels, lonlats)):
            lonlat = np.atleast_2d(lonlat)
            if 'color' not in scatter_kws or mod_cr:
                scatter_kws['color'] = colors[i]
                scatter_kws['marker'] = markers[i]
                mod_cr = True
            ax.scatter(*m(lonlat[:,0], lonlat[:,1]), label=label, zorder=12, **scatter_kws) 
        ax.legend(loc='lower left', prop={'weight':'bold', 'size':8}).set_zorder(20)
 
    else:
        for lonlat in lonlats:
            if len(lonlat):
                lonlat = np.atleast_2d(lonlat)
                s = ax.scatter(*m(lonlat[:,0], lonlat[:,1]), zorder=12, **scatter_kws)
                # plt.colorbar(s, ax=ax)
    hide_kwargs = {'axis':'both', 'which':'both'}
    hide_kwargs.update(dict([(k, False) for k in ['bottom', 'top', 'left', 'right', 'labelleft', 'labelbottom']]))
    ax.tick_params(**hide_kwargs)
 
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(1.5)
        ax.spines[axis].set_zorder(50)
    # plt.axis('off')
 
    if world:
        size = 0.35
        if us:
            loc = (0.25, -0.1) if eu else (0.35, -0.01)
            ax_ins = inset_axes(ax, width=PLOT_WIDTH*size, height=PLOT_HEIGHT*size, loc='center', bbox_to_anchor=loc, bbox_transform=ax.transAxes, axes_kwargs={'zorder': 5})
            
            scatter_kws.update({'s': 6})
            m2 = draw_map(*lonlats, labels=labels, ax=ax_ins, **scatter_kws)
            
            mark_inset(ax, ax_ins, m, m2, US_MAP, loc1=(1,1), loc2=(2,2), edgecolor='grey', zorder=3)
            mark_inset(ax, ax_ins, m, m2, US_MAP, loc1=[3,3], loc2=[4,4], edgecolor='grey', zorder=0)
 
 
        if eu:
            ax_ins = inset_axes(ax, width=PLOT_WIDTH*size, height=PLOT_HEIGHT*size, loc='center', bbox_to_anchor=(0.75, -0.05), bbox_transform=ax.transAxes, axes_kwargs={'zorder': 5})
            
            scatter_kws.update({'s': 6})
            m2 = draw_map(*lonlats, us=False, eu=True, labels=labels, ax=ax_ins, **scatter_kws)
            
            mark_inset(ax, ax_ins, m, m2, EU_MAP, loc1=(1,1), loc2=(2,2), edgecolor='grey', zorder=3)
            mark_inset(ax, ax_ins, m, m2, EU_MAP, loc1=[3,3], loc2=[4,4], edgecolor='grey', zorder=0)
 
    return m
 
 
def default_dd(d={}, f=lambda k: k):
    ''' Helper function to allow defaultdicts whose default value returned is the queried key '''
 
    class key_dd(dd):
        ''' DefaultDict which allows the key as the default value '''
        def __missing__(self, key):
            if self.default_factory is None:
                raise KeyError(key)
            val = self[key] = self.default_factory(key)
            return val 
    return key_dd(f, d)
 
 
@ignore_warnings
def plot_scatter(y_test, benchmarks, bands, labels, products, sensor, title=None, methods=None, n_col=3, img_outlbl=''):
    import matplotlib.patheffects as pe 
    import matplotlib.ticker as ticker
    import matplotlib.pyplot as plt 
    import seaborn as sns 
 
    folder = Path('scatter_plots')
    folder.mkdir(exist_ok=True, parents=True)
 
    product_labels = default_dd({
        'chl' : 'Chl\\textit{a}',
        'pc'  : 'PC',
        'aph' : '\\textit{a}_{ph}',
        'tss' : 'TSS',
        'cdom': '\\textit{a}_{CDOM}',
    })
    
    product_units = default_dd({
        'chl' : '[mg m^{-3}]',
        'pc'  : '[mg m^{-3}]',
        'tss' : '[g m^{-3}]',
        'aph' : '[m^{-1}]',
        'cdom': '[m^{-1}]',
    }, lambda k: '')
 
    model_labels = default_dd({
        'MDN' : 'MDN_{A}',
    })
 
    products = [p for i,p in enumerate(np.atleast_1d(products)) if i < y_test.shape[-1]]
 
    plt.rc('text', usetex=True)
    plt.rcParams['mathtext.default']='regular'
    # plt.rcParams['mathtext.fontset'] = 'stix'
    # plt.rcParams['font.family'] = 'cm'
 
    # Only plot certain bands
    if len(labels) > 3 and 'chl' not in products:
        product_bands = {
            'default' : [443, 482, 561, 655],
            # 'aph'     : [443, 530],
        }
 
        target     = [closest_wavelength(w, bands) for w in product_bands.get(products[0], product_bands['default'])]
        plot_label = [w in target for w in bands]
        plot_order = ['MDN', 'QAA', 'GIOP']
        plot_bands = True
    else:
        plot_label = [True] * len(labels)
        plot_order = methods
        plot_bands = False
 
        if plot_order is None:
            if 'chl' in products and len(products) == 1:
                benchmarks = benchmarks['chl']
                if 'MLP' in benchmarks:
                    plot_order = ['MDN', 'MLP', 'SVM', 'XGB', 'KNN', 'OC3']
                else:
                    plot_order = ['MDN', 'Smith_Blend', 'OC6', 'Mishra_NDCI', 'Gons_2band', 'Gilerson_2band']
            elif len(products) > 1 and any(k in products for k in ['chl', 'tss', 'cdom']):
                plot_order = {k:v for k,v in {
                    'chl'  : ['MDN', 'Gilerson_2band'],#,'Smith_Blend'],
                    'tss'  : ['MDN', 'SOLID'],#, 'Novoa'],
                    'cdom' : ['MDN', 'Ficek'],#, 'QAA_CDOM'],
                }.items() if k in products}
                plot_label = [True] * len(plot_order)
                plot_bands = True
                n_col = len(plot_order)
    plot_order = {p: ['MDN'] for p in products}
    plot_label = [True] * 4
    labels = [(p,label) for label in labels for p in products if p in label]
    print('Plotting labels:', [l for i,l in enumerate(labels) if plot_label[i]])
    assert(len(labels) == y_test.shape[-1]), [len(labels), y_test.shape]
 
    # plot_order = [p for p in plot_order if p in benchmarks]
    fig_size   = 5
    n_col      = max(n_col, sum(plot_label))
    n_row      = 1#max(1,int(not plot_bands) + int(0.5 + len(plot_order) / (1 if plot_bands else n_col)) -1)
    # if isinstance(plot_order, dict): n_row = 2#len(plot_order)
    # if plot_bands:
    #   n_col, n_row = n_row, n_col
 
    fig, axes = plt.subplots(n_row, n_col, figsize=(fig_size*n_col, fig_size*n_row+1))
    axes      = [ax for axs in np.atleast_1d(axes).T for ax in np.atleast_1d(axs)]
    colors    = ['xkcd:sky blue', 'xkcd:tangerine', 'xkcd:fresh green', 'xkcd:greyish blue', 'xkcd:goldenrod',  'xkcd:clay', 'xkcd:bluish purple', 'xkcd:reddish']
 
    print('Order:', plot_order)
    print(f'Plot size: {n_row} x {n_col}')
    print(labels)
 
    curr_idx = 0
    full_ax  = fig.add_subplot(111, frameon=False)
    full_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False, pad=10)
 
    # estimate_label = 'Estimated' #'Satellite-derived'
    # x_pre  = 'Measured'
    estimate_label = 'Measured' #'Satellite-derived'
    x_pre  = 'Estimated'
    y_pre  = estimate_label.replace('-', '\\textbf{-}')
    plabel = f'{product_labels[products[0]]} {product_units[products[0]]}'
    xlabel = fr'$\mathbf{{{x_pre} {plabel}}}$'
    ylabel = fr'$\mathbf{{{y_pre}}}$'+'' +fr'$\mathbf{{ {plabel}}}$'
    if not isinstance(plot_order, dict):
        full_ax.set_xlabel(xlabel.replace(' ', '\ '), fontsize=24, labelpad=10)
        full_ax.set_ylabel(ylabel.replace(' ', '\ '), fontsize=24, labelpad=10)
    else:
        full_ax.set_ylabel(fr'$\mathbf{{{x_pre}}}$'.replace(' ', '\ '), fontsize=24, labelpad=15)
 
    s_lbl = title or get_sensor_label(sensor).replace('-',' ')
    n_pts = len(y_test)
    title = fr'$\mathbf{{\underline{{\large{{{s_lbl}}}}}}}$' + '\n' + fr'$\small{{\mathit{{N\small{{=}}}}{n_pts}}}$'
    # full_ax.set_title(title.replace(' ', '\ '), fontsize=24, y=1.06)
 
    # if isinstance(plot_order, dict):
    #   full_ax.set_title(fr'$\mathbf{{\underline{{\large{{{s_lbl}}}}}}}$'.replace(' ', '\ '), fontsize=24, y=1.03)
 
    for plt_idx, (label, y_true) in enumerate(zip(labels, y_test.T)):
        if not plot_label[plt_idx]: continue 
 
        product, title = label 
        plabel = f'{product_labels[product]} {product_units[product]}'
 
        for est_idx, est_lbl in enumerate(plot_order[product] if isinstance(plot_order, dict) else plot_order):
            # if plt_idx >= (len(plot_order[product]) if isinstance(plot_order, dict) else benchmarks[est_lbl].shape[1]): continue
            if isinstance(plot_order, dict) and est_lbl not in benchmarks[product]: 
                axes[curr_idx].tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                axes[curr_idx].axis('off')
                curr_idx += 1
                continue 
 
            y_est = benchmarks[product][est_lbl] if isinstance(plot_order, dict) else benchmarks[est_lbl][..., plt_idx]
            ax    = axes[curr_idx]
            cidx  = int(curr_idx / n_row) #(curr_idx / n_col))# if plot_bands else curr_idx
            color = colors[cidx]
 
            # first_row = curr_idx < n_col 
            first_row = (curr_idx % n_row) == 0
            # last_row  = curr_idx >= ((n_row-1)*n_col) 
            last_row  = ((curr_idx+1) % n_row) == 0
            first_col = (curr_idx % n_col) == 0
            last_col  = ((curr_idx+1) % n_col) == 0
            print(curr_idx, first_row, last_row, first_col, last_col, est_lbl, product, plabel)
            y_est_log  = np.log10(y_est).flatten()
            y_true_log = np.log10(y_true).flatten()
            curr_idx  += 1
 
            l_kws = {'color': color, 'path_effects': [pe.Stroke(linewidth=4, foreground='k'), pe.Normal()], 'zorder': 22, 'lw': 1}
            s_kws = {'alpha': 0.4, 'color': color}#, 'edgecolor': 'grey'}
 
            if est_lbl == 'MDN':
                [i.set_linewidth(5) for i in ax.spines.values()]
                # est_lbl = 'MDN_{A}'
                # est_lbl = 'MDN-I'
            else:
                est_lbl = est_lbl.replace('Mishra_','').replace('Gons_2band', 'Gons').replace('Gilerson_2band', 'GI2B').replace('Smith_','').replace('Cao_XGB','BST')#.replace('Cao_', 'Cao\ ')
                est_lbl = est_lbl.replace('QAA_CDOM', 'QAA\ CDOM')
            if product not in ['chl', 'tss', 'cdom', 'pc'] and last_col:
                ax2 = ax.twinx()
                ax2.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False, pad=0)
                ax2.grid(False)
                ax2.set_yticklabels([])
                ax2.set_ylabel(fr'$\mathbf{{{bands[plt_idx]:.0f}nm}}$', fontsize=22)
 
            minv = -2# if product == 'cdom' else -1 # int(np.nanmin(y_true_log)) - 1 if product != 'aph' else -4
            maxv = 3#1 if product == 'cdom' else 3 # if product == 'chl' else 1 #int(np.nanmax(y_true_log)) + 1 if product != 'aph' else 1
            loc  = ticker.LinearLocator(numticks=int(round(maxv-minv+1.5)))
            fmt  = ticker.FuncFormatter(lambda i, _: r'$10$\textsuperscript{%i}'%i)
 
            ax.set_ylim((minv, maxv))
            ax.set_xlim((minv, maxv))
            ax.xaxis.set_major_locator(loc)
            ax.yaxis.set_major_locator(loc)
            ax.xaxis.set_major_formatter(fmt)
            ax.yaxis.set_major_formatter(fmt)
            fs = 22
 
            if not last_row:                   ax.set_xticklabels([])
            # elif isinstance(plot_order, dict): ax.set_xlabel(fr'$\mathbf{{{x_pre}}}$'+'' +fr'$\mathbf{{ {plabel}}}$'.replace(' ', '\ '), fontsize=18)
            # if not first_col:                  ax.set_yticklabels([])
            elif isinstance(plot_order, dict): 
                ylabel = fr'$\mathbf{{{y_pre}}}$'+'' +fr'$\mathbf{{ {plabel}}}$' + '\n' + fr'$\small{{\mathit{{N\small{{=}}}}{np.isfinite(y_true_log).sum()}}}$'
                ax.set_xlabel(ylabel.replace(' ', '\ '), fontsize=fs)
 
            valid = np.logical_and(np.isfinite(y_true_log), np.isfinite(y_est_log))
            if valid.sum():
                sns.regplot(y_true_log[valid], y_est_log[valid], ax=ax, scatter_kws=s_kws, line_kws=l_kws, fit_reg=True, truncate=False, robust=True, ci=None)
                kde = sns.kdeplot(y_true_log[valid], y_est_log[valid], shade=False, ax=ax, bw='scott', n_levels=4, legend=False, gridsize=100, color='#555')# color=color)
                # kde.collections[2].set_alpha(0)
 
            invalid = np.logical_and(np.isfinite(y_true_log), ~np.isfinite(y_est_log))
            if invalid.sum():
                ax.scatter(y_true_log[invalid], [minv]*(invalid).sum(), color='r', alpha=0.4, label=r'$\mathbf{%s\ invalid}$' % (invalid).sum())
                ax.legend(loc='lower right', prop={'weight':'bold', 'size': 18})
 
            add_identity(ax, ls='--', color='k', zorder=20)
 
            if valid.sum():
                add_stats_box(ax, y_true[valid], y_est[valid], metrics=[mdsa, sspb, slope], fontsize=18)
 
            # if first_row or not plot_bands or (isinstance(plot_order, dict) and plot_order[product][est_idx] != 'MDN'):
            if first_row or not plot_bands or (isinstance(plot_order, dict)):# and plot_order[product][est_idx] != 'MDN'):
                if est_lbl == 'BST':
                    # ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$'+'\n'+r'$\small{\textit{(Cao\ et\ al.\ 2020)}}$', fontsize=18)
                    ax.set_title(r'$\small{\textit{(Cao\ et\ al.\ 2020)}}$' + '\n' + fr'$\mathbf{{\large{{{est_lbl}}}}}$', fontsize=fs, linespacing=0.95)
                
                elif est_lbl == 'Ficek':
                    # ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$'+'\n'+r'$\small{\textit{(Cao\ et\ al.\ 2020)}}$', fontsize=18)
                    ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$' + r'$\small{\textit{\ (et\ al.\ 2011)}}$', fontsize=fs, linespacing=0.95)
                
                elif est_lbl == 'Mannino':
                    # ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$'+'\n'+r'$\small{\textit{(Cao\ et\ al.\ 2020)}}$', fontsize=18)
                    ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$' + r'$\small{\textit{\ (et\ al.\ 2008)}}$', fontsize=fs, linespacing=0.95)
 
                elif est_lbl == 'Novoa':
                    # ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$'+'\n'+r'$\small{\textit{(Cao\ et\ al.\ 2020)}}$', fontsize=18)
                    ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$' + r'$\small{\textit{\ (et\ al.\ 2017)}}$', fontsize=fs, linespacing=0.95)
 
                elif est_lbl == 'GI2B':
                    ax.set_title(fr'$\mathbf{{\large{{Gilerson}}}}$' + r'$\small{\textit{\ (et\ al.\ 2010)}}$', fontsize=fs, linespacing=0.95)
 
                elif est_lbl == 'MDN': ax.set_title(fr'$\mathbf{{{est_lbl}\ {product_labels[product]}}}$', fontsize=fs)
                else: ax.set_title(fr'$\mathbf{{\large{{{est_lbl}}}}}$', fontsize=fs)
 
            ax.tick_params(labelsize=fs)
            ax.grid('on', alpha=0.3)
 
    u_label  = ",".join([o.split('_')[0] for o in plot_order]) if len(plot_order) < 10 else f'{n_row}x{n_col}'
    filename = folder.joinpath(f'{img_outlbl}{",".join(products)}_{sensor}_{n_pts}test_{u_label}.png')
    plt.tight_layout()
    # plt.subplots_adjust(wspace=0.35)
    plt.savefig(filename.as_posix(), dpi=100, bbox_inches='tight', pad_inches=0.1,)
    plt.show()