%% Antarctic Ice Sheet grounding line discharge from 1996 through 2023
% Authors: Benjamin Davison, Anna Hogg, Thomas Slater, Richard Rigby
% Institution: University of Leeds
%
% Table of Contents
% Overview
%     Definitions
%         Grounding line discharge:
%         Basins:
%         What is in the data package?
% Using the data
%     Variables
% MATLAB examples
%         First look
%         Antarctic grounding line discharge
%         Example basins
%         Discharge trends
% Python examples
%         Antarctic grounding line discharge
%
%
%% Overview
% We provide grounding line discharge time-series for the Antarctic Ice Sheet and all of its drainage basins, from 1996 through 2023. We have provided several different estimates of grounding line discharge, calculated using different bed topography products, different estimates of firn air content, and calculated through different flux gates. In our manuscript, we focus on one bed topography product (FrankenBed) and one firn model (the IMAU FDM) and use the inverse-error weighted mean discharge through all of the flux gates. Henceforth, we refer to this as our 'main discharge estimate'. There isn't a strong scientific justification for preferring the IMAU FDM over the other firn models, and there is only a modest justification for preferring FrankenBed over the other bed topography products, so we have provided the others as well in case they better suit your research.
% You can find our main discharge estimate in AIS_discharge_FrankenBed.nc (the variable is discharge_mean), or if you prefer .csv, then it's in AIS_discharge_FrankenBed_IMAU.csv.
%
% Before we get going, here are a few quick definitions (see the manuscript for more details):
%
%% Definitions
%% Grounding line discharge:
% Grounding line discharge, D, is the flux of ice across the grounding line, given by:
% D=\textrm{VHw}\rho,
% where V is the ice velocity, H is the ice thickness, \rho is the density of ice and w is the width of the region over which V and H are given. In our dataset, we extract V and H over 100 m or 200 m pixels and interpolate across all data gaps, so as to provide the total discharge.
%% Basins:
% Our definition of 'drainage basins' encompasses all of the usual ways you may want to split up Antarctica, including: West Antarctica, East Antarctica and the Antarctic Peninsula, down to individual glaciers on the Antarctica Peninsula. As such, some of our basins cover the same area - for example, Pine Island Glacier falls within West Antarctica, but we provide both. Similarly, some people prefer to use the Zwally drainage basins (https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems) whilst others prefer to use the Rignot basins (https://nsidc.org/data/nsidc-0709/versions/2), so we have provided discharge time-series through both sets of basins.
%
%% What is in the data package?
% Our primary dataset is four NetCDF files of grounding line discharge:
% AIS_discharge_FrankenBed.nc
% AIS_discharge_FrankenBedAdj.nc
% AIS_discharge_BedMachine.nc
% AIS_discharge_BedMap2.nc
% These four files have a corresponding NetCDF file with an "_ns" suffix, meaning "no shift". In our manuscript, we describe an offset in velocities between lower resolution and higher resolution velocity datasets. In our primary dataset, we have increased the velocities in the lower resolution datasets so that they are better aligned with the higher resolution data (though you can still see a jump at the change in data sources, which is partly due to the differing temporal resolution of the data sources). In the "_ns" versions, we do the opposite: we reduce the velocities of the higher resolution data to better align with the lower resolution data.
% Basins
% We also provide information about the basins over which discharge has been calculated. This information is provided in NetCDF, shapefile and geopackage formats.
% Discharge in CSV format
% The discharge data are also provided in .csv format, with the naming convention:
% AIS_<VAR>_gate_<DISTANCE>m_<BEDPRODUCT>_<FIRN>.csv, where:
% <VAR> can be "discharge", "discharge_error", "discharge_standard_error", "discharge_coverage", "SMB_correction", "FDM_correction" and "dH correction".
% <DISTANCE> is an integer from 0 to 3000 giving the distance in metres upstream from the first flux gate. These correspond to the flux gate through which discharge was calculated. The csv's containing gate mean quantity do not include this name part.
% <BEDPRODUCT> can be FrankenBed, FrankenBedAdj, BedMaching, or BedMap2.
% <FIRN> can be IMAU, GSFC or GEMB.
% Balance discharge
% The discharge required to balance the long-term SMB over each basin. Provided as NetCDF and csv formats.
% Discharge statistics
% A .mat file containing discharge trends for each basin from 1996 to 2023. To help with reproducing Figure 8.
% Flux gates
% A .mat file, a NetCDF, a geopackage and a shapefile containing the coordinates of each of the 16 flux gates used in the analysis. The coordinates are in South Polar Stereographic metres, EPSG 3031: https://epsg.io/3031
% FrankenBedAdj
% FrankenBedAdj is only calculated at each flux gate pixel - it is not a grid. We provide a NetCDF containing the pixel-wise thickness in FrankenBed and the 'Adj' (or adjustment) applied to FrankenBed in order to create FrankenBedAdj. Since FrankenBedAdj is calculated based on the observed velocities and surface elevation change rates, we provide two versions: one with and without the velocity shift described above.
% Reference grids
% We include several gridded datasets used or generated in the analysis:
% FrankenBed.tif - essentially BedMachine but with the Huss and Farinotti (2014) data on the Peninsula and the Cui et al. (2020) data in Princess Elizabeth Land. See manuscript for details.
% FrankenBed_error.tif - error in FrankenBed
% reference_vx.tif - reference easting ice velocity
% refernece_vy.tif - reference northing ice velocity
% reference_vx_error.tif - error in easting ice velocity
% reference_vy_error.tif - erorr in northing ice velocity
% fac_clim_2D.csv - 1979-2008 mean firn air content from the IMAU FDM
% sec_quadratic_1992_2023.nc - grids of the coefficients to the polynomial fitted to the surface elevation change data
% For the SEC grids, the polynomial of the form: y = ax^2 + bx + c , is fitted to the altimetry timeseries in 5x5 km regions. When reconstructing the time series use e.g. x = [0:1/12:32] for t = 1991-2023. One can estimate the uncertainty by taking the max/min dh from the coefficients +/- their uncertainty. i.e. (courtesy of Thomas Slater):
% t_sampling_original = 140/365
% t_sampling_hires = 1/12
% t_monthly = np.arange(1,32, t_sampling_hires) # 1 = year 1991
% model_dh = (t_monthly**2 * sec_quadratic) + (t_monthly * sec_linear) + sec_icept
% model_dh_upper = ((t_monthly**2 * (sec_quadratic + sec_quadratic_uncertainty)) + (t_monthly * (sec_linear + sec_linear_uncertainty)) + (sec_icept + sec_icept_uncertainty))
% model_dh_lower = ((t_monthly**2 * (sec_quadratic - sec_quadratic_uncertainty)) + (t_monthly * (sec_linear - sec_linear_uncertainty)) + (sec_icept - sec_icept_uncertainty))
% model_dh_uncertainty = (t_sampling_hires  / t_sampling_original) * (((model_dh_upper - model_dh) + (model_dh - model_dh_lower)) / 2)
%
% Where (t_sampling_hires  / t_sampling_original)  is a scaling factor based on the ratio between the temporal sampling of the higher resolution x (monthly here) and the original input x (140 days) which the model is based on, otherwise the uncertainty scales based on the input x
%
% Summary figs
% .pngs summarising the discharge and related data for each basin. Lots of interesting things going on.
%
%% Using the data
% This section focuses on using our main discharge estimate presented in the manuscript - calculated using our mishmash of bed topography products ('FrankenBed') along with the IMAU FDM and averaged through all flux gates.
% The main discharge dataset is the variable discharge_mean in the NetCDF file AIS_discharge_FrankenBed.nc. discharge_mean is an MxN array corresponding to M basins and N measurement times. That is, the first row in discharge_mean is a time-series of grounding line discharge corresponding to the first basin. The basins are given by the variable 'basin_name', which have a corresponding 'basin_type', basin_type_ID' and 'basin_ID'. You can find more information about these basins in the basins output (shapefiles, geopackages and netcdfs of each basin and basin type).
% For example, the basin 'Antarctica' is the first basin in the 'basin_name' variable. Therefore, the first row in 'discharge_mean' is a timeseries of grounding line discharge for Antarctica. Similarly, the second basin name is called 'West_Antarctica_Rignot', so the second row in 'discharge_mean' is a discharge time-series for West Antarctica as defined by the Rignot (MEaSUREs) basin dataset.
% In general, to get grounding line discharge for any basin or region (including Antarctica), you just need to work out which row your basin of interest is on. The easiest way to do this is the load the basin geopackage or shapefile into QGIS or your chosen software, zoom to your basin of interest (assuming you know where it is!) and get the basin ID, which will always correspond to the row you need. Alternatively, you can get the basin name and search for that name in the list of basin names in the NetCDF file.
% Note that the basin names stored in the NetCDF files contain upper and lowercase letters, some underscores and some basins are numbered where names were duplicated between basin datasets. Therefore, even if you know a basin by one name, you might need to take a bit of care to find the matching basin in this dataset.
%
%% Variables
% Each discharge NetCDF file contains the main discharge time-series described above - discharge_mean - along with several other discharge time-series and errors. These are:
% Variable
% Size
% Dimensions
% Datatype
% Units
% Description
% basin_name
% 966x1
% basin
% string
% N/A
% The name of each basin
% basin_type
% 966x1
% basin
% string
% N/A
% The type of basin, which
% can be Full ice sheet,
% MEaSUREs regional,
% MEaSUREs glaciers,
% Zwally, Cook, ice shelf,
% Ice Sheet Rignot,
% Ice Sheet Zwally.
% basin_type_ID
% 966x1
% basin
% single
% N/A
% A number corresponding
% to each of the basin types
% in basin_type
% basin_ID
% 966x1
% basin
% single
% N/A
% A number from 1 to 966
% to identify each basin
% time
% 170x1
% time
% double
% 'days since
% 1950-01-01
% 00:00:00'
% Measurement time.
% Proleptic Gregorian
% distance
% 16x1
% gate
% single
% m
% Relating to flux gates.
% Distance upstream from
% the first flux gate.
% datasource
% 170x1
% time
% single
% N/A
% A number corresponding
% to each velocity datasource.
% ice velocity dataset source.
% 1 = MEaSUREs annual;
% 2 = Sentinel-1;
% 3 = ERS/BISICLES;
% 4 = ITS_LIVE;
% 5 = MEaSUREs multiyear;
% 6 = ITS_LIVE ASE sim;
% 7 = Joughin Pine Island TSX
% discharge_all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% Timeseries of grounding line
% discharge for every basin
% through each individual flux
% gate.
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_error
% _all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% The corresponding error for
% discharge_all_gates
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_coverage
% _all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% The percentage of
% discharge_all_gates that was
% actually observed, rather
% than based on interpolated
% velocities
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_mean
% 966x170
% basin x time
% double
% Gt per year
% The average discharge
% through all 16 flux gates.
% Gate discharge is weighted
% by 1/discharge_error_all_gates
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_standard
% _error
% 966x170
% basin x time
% double
% Gt per year
% Standard error of the gate
% mean discharge. Defined as
% one sixteenth of the standard
% deviation of discharge across
% each of the 16 flux gates.
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_error
% 966x170
% basin x time
% double
% Gt per year
% The average discharge error
% through all 16 flux gates.
% Gate errors are weighted by
% 1./discharge_error_all_gates.
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_coverage
% 966x170
% basin x time
% double
% Gt per year
% The percentage of discharge
% observed, on average through
% all gates.
% Calculated using the IMAU
% FDM and with all surface
% height corrections applied.
% discharge_all_gates
% _GSFC
% 966x170x16
% basin x time x gate
% double
% Gt per year
% Identical to
% discharge_all_gates, but using
% the GSFC-FDMv1.2
% discharge_mean_GSFC
% 966x170
% basin x time
% double
% Gt per year
% Identical to
% discharge_mean, but using
% the GSFC-FDMv1.2
% discharge_all_gates
% _GEMB
% 966x170x16
% basin x time x gate
% double
% Gt per year
% Identical to
% discharge_all_gates, but using
% the FEMB
% discharge_mean_GEMB
% 966x170
% basin x time
% double
% Gt per year
% Identical to
% discharge_mean, but using
% the GEMB
% SMB_correction
% _all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% The change in surface height
% due to SMB between each
% flux gate and the grounding line
% FDM_correction
% _all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% The firn air content correction
% from the IMAU FDM, for each
% flux gate
% dH_correction_all_gates
% 966x170x16
% basin x time x gate
% double
% Gt per year
% The change in surface height
% applied from the quadratic
% fit to surface elevation change
% observations
%
% Note that the time dimension is likely to grow as new data are added.
%
%
%% MATLAB examples
%% First look
% Before we dive into loading and plotting the discharge data, let's orientate ourselves:

% First, make sure the data are on the MATLAB path
addpath(genpath('.'));
addpath(genpath('D:/work/Fellowship_Leeds/Data')); % I keep some handy MATLAB tools in here, which I've linked to whenever they are used below

% Load grounding line discharge
fn = 'AIS_discharge_FrankenBed.nc';

% What variables do we have?
ncdisp(fn)

% The FrankenBed discharge NetCDF file contains some information about the basins - basin_name, basin_type, basin_type_ID and basin_ID - and various different discharge variables, along with some reference information, including time, datasource and distance. These variables are detailed in the NetCDF file and in the table above. Their definitions are identical for each bed topography-discharge product.
% In this version of the dataset, there are 966 individual basins. Each basin has it's own discharge timeseries. We can see which basins are included by looking at the list of basin names:

names = ncread(fn,'basin_name') % a list of strings containing the names of all of the basins. There are 966 different names, each corresponding to a basin and a row in the discharge array below.

% As you can see, I have used the term 'basins' rather loosely, refering to almost any way previous studies have split up the ice sheet. Note that the only island in this datset is Wilkins Island - so my definition of 'Antarctica' may differ to others.
% Nevertheless, the first 'basin' is "Antarctica", the second is "West_Antarctica_Rignot" (meaning the West Antarctic Ice Sheet as defined by the Rignot or MEaSUREs basins), the third is "East_Antarctica_Rignot" and so on and so forth...
% Each of these basins has it's own discharge time-series. Let's load them to take a look:

time = ncread(fn,'time') + datenum(1950,1,1,0,0,0); % time in MATLABs datenum format (i.e. the number of days since the stoke of midnight on Jan 1st 0000). 
discharge = ncread(fn,'discharge_mean'); % An array of grounding line discharge, in Gt/yr, where each row corresponds to a basin or region and each columns corresponds to a measurement time.

disp(['Discharge is an ' num2str(size(discharge,1)) 'x' num2str(size(discharge,2)) ' array']);
disp('where the first dimension corresponds to each basin and the second to each measurement time');

% So the first row of discharge measurments corresponds to the discharge time-series for "Antarctica", the second row corresponds to the time-series for "West_Antarctica_Rignot" and so on and so forth. The same structure is used for each of the various discharge variables.
% To make this even clearer, we can tabulate the data:

TT = array2table([round(doy(time,'decimalyear'),1),round(discharge',1)],'VariableNames',["Time"; names]')

% You will notice that there are some repeated dates, sometimes with repeated measurements. This is because I have used multiple velocity products, some of which have the same mid-date. They produce the same measurement only where I have interpolated the velocity through that measurement time in that basin.
%
%% Antarctic grounding line discharge
% If we want to get grounding line discharge for all of Antarctica, like in Figure 6 of the manuscript, then in MATLAB we just need to do:

% Load grounding line discharge
fn = 'AIS_discharge_FrankenBed.nc';
names = ncread(fn,'basin_name'); % a list of strings containing the names of all of the basins. There are 966 different names, each corresponding to a basin and a row in the discharge array below.
time = ncread(fn,'time') + datenum(1950,1,1,0,0,0); % time in MATLABs datenum format (i.e. the number of days since the stoke of midnight on Jan 1st 0000). 
discharge = ncread(fn,'discharge_mean'); % An array of grounding line discharge, in Gt/yr, where each row corresponds to a basin or region and each columns corresponds to a measurement time.

% So now we have an array of grounding line discharge, in Gt/yr. 'discharge' is an MxN array corresponding to M basins, given by 'names' and N measurement times, given by 'time'.
% If we just want the discharge for Antarctica, then we need to find Antarctica in the list of basin names.

idx = find(strcmp('Antarctica',names)==1); % which row is Antarctica? It should be the first one
AIS_discharge = discharge(idx,:); % Get the discharge time-series for Antarctica

% This gives us a vector of grounding line discharge, in Gt/yr, for all of Antarctica. which we can plot with

figure('color','w','units','inches','position',[1 1 7 3]);
scatter(time,AIS_discharge,15,'k','filled'); 
datetick('x','keeplimits'); 

% That's fine, but it doesn't look much like Figure 6a in the manuscript, which also includes the error in the discharge estimate. Let's add that too.

se = ncread(fn,'discharge_standard_error'); % An array of grounding line discharge standard error, in Gt/yr, where each row corresponds to a basin or region and each columns corresponds to a measurement time.
AIS_confbounds = 1.96.*se(idx,:); % multiply the standard error by 1.96 to estimate the 95% confidence limits

% Add those to the plot
hold on;
shadedErrorBar(time,AIS_discharge,AIS_confbounds); % https://uk.mathworks.com/matlabcentral/fileexchange/26311-raacampbell-shadederrorbar

% That's a little bit better, but so far we have only plotted the average discharge through all of the flux gates. What if we want to also see the discharge through each individual gate?

discharge_all_gates = ncread(fn,'discharge_all_gates'); % An MxNxP array of grounding line discharge, in Gt/yr, where each row corresponds to a basin or region, each columns corresponds to a measurement time and each page corresponds to a different flux gate
distance = ncread(fn,'distance'); % distance in metres of each flux gate upstream from the first gate
for i = 1:length(distance)
    plot(time,squeeze(discharge_all_gates(idx,:,i)),'Color',[.7 .7 .7]);
end

% make it pretty
box on
ylabel('Discharge (Gt yr^{-1})')
set(gca,'XLim',[datenum(1996,1,1) datenum(2023,10,1)],...
    'XTick',datenum(1996:2:2024,1,1),...
    'XTickLabel',{'','1997','','2001','','2005','','2009','','2013','','2017','','2021',''});

% With that time-series, it is straightforward to reproduce the discharge changes quoted in the manuscript:

disp(['Antarctic discharge in 1996 was: ' num2str(AIS_discharge(1)) ' +/- ' num2str(AIS_confbounds(1)) ' Gt/yr ']);
disp(['Antarctic discharge in ' datestr(time(end),'mmm yyyy') ' was: ' num2str(AIS_discharge(end)) ' +/- ' num2str(AIS_confbounds(end)) ' Gt/yr ']);

%%
%%
%% Example basins
% In Figure 7 of the manuscript, we show lots of example discharge time-series for basins of different sizes and behaviours. Let's recreate one of those panels.

% load the data
fn = 'AIS_discharge_FrankenBed.nc';
names = ncread(fn,'basin_name'); % a list of strings containing the names of all of the basins. There are 966 different names, each corresponding to a basin and a row in the discharge array below.
time = ncread(fn,'time') + datenum(1950,1,1,0,0,0); % time in MATLABs datenum format (i.e. the number of days since the stoke of midnight on Jan 1st 0000). 
discharge = ncread(fn,'discharge_mean'); % An array of grounding line discharge, in Gt/yr, where each row corresponds to a basin or region and each columns corresponds to a measurement time.
se = ncread(fn,'discharge_standard_error'); % An array of grounding line discharge, in Gt/yr, where each row corresponds to a basin or region and each columns corresponds to a measurement time.

% Set your target basins
target_basins = {'Bleriot Glacier','Niepce Glacier','Luke Glacier','Pyke Glacier'};

% make some pretty colours
cmap = [0.8500 0.3250 0.0980;
    0.4940 0.1840 0.5560;
    0.9290 0.6940 0.1250;
    0 0.4470 0.7410];

% plot the target basins
figure('color','w','units','inches','position',[1 1 7 3]);
for i = 1:numel(target_basins)
    idx = find(strcmp(target_basins{i},names)==1); % get the index for this basin
    scatter(time,discharge(idx,:),8,cmap(i,:),'filled'); hold on
    ptch = shadedErrorBar(time,discharge(idx,:),1.96.*se(idx,:)); % 95% confidence bounds
    ptch.mainLine.Color = 'none';
    ptch.edge(1).Color = 'none';
    ptch.edge(2).Color = 'none';
    ptch.patch.FaceColor = cmap(i,:);
    ptch.patch.FaceAlpha = 0.3;
end

% make it pretty
box on
ylabel('Discharge (Gt yr^{-1})')
set(gca,'XLim',[datenum(1996,1,1) datenum(2023,10,1)],...
    'XTick',datenum(1996:2:2024,1,1),...
    'XTickLabel',{'','1997','','2001','','2005','','2009','','2013','','2017','','2021',''},...
    'YTick',0.5:0.25:1,'YTickLabel',{'','0.75','1'},'YLim',[0.56 1.1]);

%
% If you want to plot some different basins, you will need to browse the basins (shapefiles, geopackage or netcdf) in your chosen software and extract it's name or basin ID. With that, you can pull out the discharge time-series for that basin using simple indexing.
%
%
%% Discharge trends
% In Figure 8 of the manuscript, I show some discharge trends in some basins around Antarctica. Here's how.
% First, we need the information for all basins. I prefer to work with the basins in a MATLAB structure, which the tool readBasinInfo does for me:

% load the data
basins = readBasinInfo; % basin structure
bQ = ncread('AIS_discharge_FrankenBed.nc','discharge_mean'); % gete-mean basin discharge
time = ncread('AIS_discharge_FrankenBed.nc','time')+datenum(1950,1,1,0,0,0); % time in MATLAB datenum format

% Now, in every basin let's calculate the trend, and the trend with respect to the time-average discharge, through the full time period.

% pre-allocate
Qtr_lin = zeros([numel(basins.name),1]);
Qtr_rlin = zeros([numel(basins.name),1]);
R_SLF = zeros([numel(basins.name),1]);
P_SLF = zeros([numel(basins.name),1]);


for i = 1:numel(basins.name)
    ts = bQ(i,:);
    if any(ts~=0)
        % Get linear fit
        [p, ~, mu] = polyfit(time(:), ts(:), 1);
        fit_ts = polyval(p, time(:), [], mu);
        Qtr_lin(i) = diff([fit_ts(end), fit_ts(1)])./(diff([max(time(:)), min(time(:))])/365.25); % the linear trend in Gt/yr/yr
        Qtr_rlin(i) = Qtr_lin(i)/nanmedian(fit_ts).*100; % the relative trend, in %/yr
        [R_SLF(i), P_SLF(i)] = corr(time(:), ts(:), 'type', 'Pearson'); % Pearsons R and p values
    end
end

% Now that we have the trends we can plot the result.

% prepare colormap
cpts = -2.5:.025:2.5;
cmap = cmocean('balance',length(cpts));

% Plot for MEaSUREs glacier basins
ind = find(basins.type_ID==2); % gert the index of those basins
figure('color','w','units','inches','position',[1 1 6 6]);
antbounds('coast','color','k'); hold on; % https://uk.mathworks.com/matlabcentral/fileexchange/60246-antarctic-boundaries-grounding-line-and-masks-from-insar
antbounds('coast','polyshape','facecolor','w');
%   - colour each basin based on its relative trend
for i = 1:numel(ind)
    [~,color_ind] = findNearest(cpts,Qtr_rlin(ind(i))); % https://uk.mathworks.com/matlabcentral/fileexchange/30029-findnearest-algorithm
    fill(basins.x{ind(i)}(1:end-1),basins.y{ind(i)}(1:end-1),cmap(color_ind,:),'edgecolor','none');
end
%   - plot the basin outline
for i = 1:length(ind)
    plot(basins.x{ind(i)},basins.y{ind(i)},'k','linestyle','-','marker','none');
end
antbounds('coast','color','k'); hold on; % replot the coastline
axis off
colormap(cmocean('balance'));

% add a colorbar
cb = colorbar('location','northoutside');
colormap(cmocean('balance'));
cb.Ticks = 0:.5:1;
cb.TickLabels = linspace(min(cpts),max(cpts),length(cb.Ticks));
cb.FontWeight = 'bold';
pos = cb.Position;
cb.Position = [pos(1)+0.05 pos(2)-0.06 pos(3)*0.25 pos(4)*0.5];
xlabel(cb,'Discharge trend (% yr^{-2})','fontsize',10,'fontname','Helvetica','fontweight','bold');

% And you can repeat that, but set the target basin_ID to 4 in order to plot the Peninsula basins.
%
%
%% Python examples
% MATLAB really doesn't like me doing this, so all of the code below is commented out, but it should work in Python if you copy it to your own script.
%% Antarctic grounding line discharge
% If we want to get grounding line discharge for all of Antarctica, like in Figure 6 of the manuscript, then in Python we just need to do:

%{
% import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from datetime import datetime, timedelta
from matplotlib.dates import num2date, date2num

# Load NetCDF file
fn = 'AIS_discharge_FrankenBed.nc'
dataset = Dataset(fn, 'r')

# Extract data
names = dataset.variables['basin_name'][:]
time = dataset.variables['time'][:] + date2num(datetime(1950, 1, 1, 0, 0, 0))  # Convert to MATLAB datenum format
discharge = dataset.variables['discharge_mean'][:]
se = dataset.variables['discharge_standard_error'][:]
discharge_all_gates = dataset.variables['discharge_all_gates'][:]
distance = dataset.variables['distance'][:]

# Find row index for Antarctica
idx = np.where(names == 'Antarctica')[0][0]

# Get discharge time-series for Antarctica
AIS_discharge = discharge[idx, :]

# Nested function for shaded error bar plot
def shaded_error_bar(ax, x, y, yerr, color, alpha=0.5, label=None):
    ax.plot(x, y, color=color, label=label)
    ax.fill_between(x, y - yerr, y + yerr, color=color, alpha=alpha)

# Plotting
fig, ax = plt.subplots(figsize=(7, 3))
ax.scatter(num2date(time), AIS_discharge, s=15, color='k', marker='o', label='Antarctica')
ax.xaxis_date()
plt.xticks(rotation=45)
ax.set_xlabel('Time')
ax.set_ylabel('Discharge (Gt/yr)')

# Estimate 95% confidence limits
AIS_confbounds = 1.96 * se[idx, :]

# Use the shaded error bar function
shaded_error_bar(ax, num2date(time), AIS_discharge, AIS_confbounds, color='k', label='_nolegend_')

# Plot discharge for all gates
for i in range(len(distance)):
    ax.plot(num2date(time), discharge_all_gates[idx, :, i], color=[0.7, 0.7, 0.7])

# Show the plot
plt.legend()
plt.show()

# Close the NetCDF file
dataset.close()
%}

%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%