Preamble

This notebook deals with compiling all available estimates and measurements of vertical nitrate fluxes (in the Arctic Ocean).

All figures exported from this notebook are prefixed by FIGURE_NO3-COMP_.

In [ ]:
%load_ext autoreload
%autoreload 2
%run imports.py

dims = dict(
    FN = hv.Dimension('FN', label='Nitrate flux', unit='mmol N m⁻² d⁻¹', range=(.005, 10)),
    no3_sfc_winter_uM = hv.Dimension('no3_sfc_winter_uM', label='Pre-bloom sfc. NO₃ ', unit='µM', range=(-.5, 13)),
    strat=hv.Dimension('strat', label='Brunt-Väisälä frequency', unit='s⁻²')
)

Maps

In [ ]:
df = pd.read_csv('../data/fn-compilation.csv')

options_bk = [
    opts.Points(
        backend='bokeh',
        color=dim('FN'),
        tools=['hover'], cmap=default_cmap, colorbar=True, 
        size=dim('samplesize').categorize({'single': 9, 'aggregate': 20}),
        line_color='k',
        width=500, height=400, show_legend=False),
    opts.Feature(scale='50m'),
]
             

options_mpl = translate_options(options_bk, bokeh2mpl)
options = options_mpl + options_bk

data = gv.Points(df, 
                 kdims=['Longitude', 'Latitude'], 
                 vdims=['FN','samplesize', 'Reference', 'Season'], 
                 crs=ccrs.PlateCarree()
                )

l = land.clone() * data
# l = data
l = l.redim.range(Latitude=(60,90), Longitude=(-180, 180))

ll = (
    l
    .opts(*(opts_map+options), clone=True)
    .opts(
        opts.Points(
            backend='bokeh', 
            color=np.log10(dim('FN')), hooks=[logcolor_ticks([1e-2, 1e-1, 1, 5])]
        )
    )
)
hv.renderer('bokeh').theme = theme
hv.save(ll, '../nb_fig/FIGURE_FN-COMP_map.html')
ll
In [ ]:
def finalize(fig):
    fig.set_size_inches(8,6)
    ax = fig.axes[0]
    ax.background_img()
    ax.set_title(dims['FN'].pprint_label)
    
    ax = fig.axes[1]
    ax.yaxis.set_label_text('')
    ax.yaxis.set_ticks([1e-2, 1e-1, 1, 5])
    ax.yaxis.set_ticklabels([1e-2, 1e-1, 1, 5])
    
fig = mplrender_map(
    l
    .opts(
        *opts_map_mpl, *options_mpl,
        opts.Points(
            logz=True, 
            cmap=cmocean.cm.matter_r, zlim=(1e-2, 5), clipping_colors={'min': 'grey', 'NaN': 'grey'}, 
            backend='matplotlib'
        ),
        opts.Feature(scale='50m'),
        clone=True
    ),
    '../nb_fig/FIGURE_FN-COMP_map',
    hooks=[finalize],
)
fig

Seasonal cycle

In [ ]:
df = pd.read_csv('../data/fn-compilation.csv')

df_nonwrap = df.loc[df.validity_start<df.validity_end]
df_wrap = df.loc[df.validity_start>df.validity_end]

df = pd.concat([
    df_nonwrap,
    df_wrap.assign(measured_end=13, validity_end=13),
    df_wrap.assign(measured_start=1, validity_start=1)
])

# reshape the start and end times of the different periods into a tidy tabular format 
# for easier ingestion in holoviews
def melt_df(x):
    dfr = df.melt(id_vars=['Reference', 'FN', 'samplesize', 'overturning'], value_vars=['measured_'+x, 'validity_'+x], 
                  value_name=x, var_name='type')
    dfr.loc[dfr.type.str.contains('measured'), 'type']='measured'
    dfr.loc[dfr.type.str.contains('valid'), 'type']='valid'
    return dfr

df = melt_df('start').join(melt_df('end')[['end']])

df = df.replace(dict(
    aggregate='Larger-scale aggregate', single='Few stations', overturning='Winter overturning', perennial='Perennial stratification', 
    measured='Measured', valid='Stipulated'
))

df.loc[df.start==df.end, 'end'] += 1

df.head()
In [ ]:
t = pd.date_range('2000-1-1', '2001-1-1', freq='2MS')
xticks = list(zip(t.month+.5, t.month_name().str[:3]))

od = {
    'overturning': {'color': ['orange', 'blue']},
    'type': {'line_width': [6, 3]},
    'samplesize': {'alpha': [.3, 1]}
}

options = (
    opts.NdOverlay(legend_position='right', show_legend=False, width=800),
    opts.Segments(logy=True, padding=.1, xticks=xticks, show_title=False, title_format='',
                  yformatter='%g', **faceted_legend_opts(od), xlabel=''),
)

l = hv.Dataset(df.assign(fn2=df.FN), ['start', 'FN', 'end', 'fn2'] + ['overturning', 'type', 'samplesize']).to(Segments).overlay()
l = l.redim(**dims)
l = l.opts(*options) + faceted_legend(l, od, hw={'height':200, 'width':250})

fname = '../nb_fig/FIGURE_FN-COMP_chart_seasonal_cycle'
ll = l.opts(toolbar=None, clone=True)
hv.save(ll, fname, fmt='png')
save_bokeh_svg_multipanel(ll, fname+'.svg', orientation='v')
hv.save(l, fname, fmt='html')
l

Surface nitratre concentration vs. stratification vs. bathymetry

FN- winter NO3

In [ ]:
df = pd.read_csv('../data/fn-compilation.csv')
df.overturning = df.overturning.str.capitalize()

options = [
    opts.Scatter(
        tools=['hover'],
        color='black',
        marker=hv.Cycle(['triangle', 'square']),
        show_grid=False, 
        legend_position='top_left', show_legend=True,
        frame_width=220, aspect=1, #color=dim('label'),
        size=10, 
        logx=True, xformatter='%g'),
]

l = (
    hv.Dataset(df, ['FN', 'overturning'], ['no3_sfc_winter_uM', 'Reference'])
    .to(hv.Scatter).overlay('overturning')
    .redim(**dims)
)

l = l.opts(*options).redim.range(FN=((5e-3,6))).relabel('A')

hv.renderer('bokeh').theme = theme
fname = '../nb_fig/FIGURE_FN-COMP_chart_vs_NO3_panelA'
hv.save(l, fname, fmt='html')
ll = l.opts(toolbar=None, clone=True)
hv.save(ll, fname, fmt='png')
save_bokeh_svg(ll, fname+'.svg')
panel_A = l.clone()
panel_A

Panel B: Surface nitratre vs. stratification grouped by bathymetry

In [ ]:
unicode_dict = {i:j for i, j in zip(range(10), list('⁰¹²³⁴⁵⁶⁷⁸⁹'))}

def format_exponent_unicode(n):
    if n>0:
        return '10'+unicode_dict[n]
    elif n<0:
        return '10⁻'+unicode_dict[abs(n)]
In [ ]:
#load data
df = (
    pd.read_csv('../data/no3-compilation/database-per-stn.csv', parse_dates=['date'])
    .rename(columns=dict(N2_30_60='strat'))
)
df = df.assign(logstrat=np.log10(df.strat))
df = df.loc[df.month.between(3, 5)]

# add bathymetry from IBCAO
ds = xr.load_dataarray('/Users/doppler/database/IBCAO/IBCAO_1min_bathy.nc')
df['btm'] = - ds.sel(lat=df.lat.values, lon=df.lon.values, method='nearest').values.diagonal()

# bin data by stratification and bathymetry
df = df.dropna(subset=['reg_name'])
df['btm_cat'] = pd.cut(df.btm, bins=[0, 200, 1500, np.inf], labels=['Shelf', 'Slope', 'Basin'])

strat_bins = np.arange(-6, -3, .5)
df['strat_cat'] = pd.cut(df.logstrat, bins=strat_bins)

# define data: surface nitrate conc. as function of stratification and bottom depth category
sc = df.hvplot.scatter('logstrat', 'ntr0', by=['btm_cat']).opts()
btm_cats = sc.data.keys()

def binavg(el):
    el = bin_average(el, bins=strat_bins, avg_fun=np.nanmean)
    return el.to(hv.Curve).clone(data=el.dframe().dropna())# * el

# colors = hv.Cycle(colorcet.b_glasbey_hv)
colors = hv.Cycle(colorcet.b_glasbey_category10)
# colors = hv.Cycle(['#282E3E', '#3F8196', '#41E2DD'])
# colors = hv.Cycle(['#2E6353', '#D1C461', '#E78A9C'])
options = [
    opts.Scatter(color=colors, tools=['hover']),
    opts.Curve(
        color=colors, line_width=5, frame_height=220, frame_width=300, tools=['hover'],
        xticks=[(logx, format_exponent_unicode(logx)) for logx in range(-6, -2)],
        xlim=(-6,-3.5),
        ylabel='',
    ),
    opts.RGB(),
    opts.Contours(
        show_legend=False,
        cmap=hv.Cycle([[c] for c in colors.values]),
        alpha=dim('Density').norm(), 
        line_dash='', line_width=3*np.sqrt(dim('Density').norm()),
    ),
    opts.Overlay(legend_position='bottom_left'),
]

avg = hv.NdOverlay({x: binavg(sc[x]) for x in btm_cats})

kde_args = dict(
#     levels=np.arange(0.01, 0.5, 0.02),
    bandwidth=0.3
)
kde = hv.NdOverlay({var: bivariate_kde(sc[var], **kde_args) for var in btm_cats})

l = kde * avg
l = l.opts(*options).redim(ntr0='no3_sfc_winter_uM', logstrat='strat').redim(**dims)

panel_B = l.clone().relabel('B')
panel_B

The overall statistical significance is hard to assess, but a category-wise box/whisker plot also consistently shows the ordering Basin < Shelf < Slope:

In [ ]:
hv.BoxWhisker(df, ['strat_cat', 'btm_cat'], 'ntr0').opts(xrotation=45, width=1000)

Merge panels

In [ ]:
l = panel_A + panel_B

fname = '../nb_fig/FIGURE_FN-COMP_chart_vs_NO3'

hv.save(l, fname, fmt='html')
save_bokeh_svg_multipanel(l, fname+'.svg')
os.system(f'convert -density 200 {fname}.svg {fname}.png')
In [ ]:
for panel in ['A', 'B']:
    os.system(f'rm {fname}_panel{panel}*')