freshwater_discharge_depth.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Freshwater discharge at depth in southeast Greenland fjords\n",
    "\n",
    "## Purpose\n",
    "This notebook generates figures and spreadsheets of freshwater discharge as a function of depth for each fjord in the study, and cumulatively across all fjords. Throughout this notebook, 'discharge' refers to freshwater discharge only (no solid discharge).\n",
    "\n",
    "The figures show, for each fjord, the discharge above/below a threshold depth, and binned by depth across the whole depth range, for discharge originating from either land or ice basins. There are also depth-binned figures that are cumulative across all fjords. \n",
    "\n",
    "The spreadsheets detail the discharge as a function of depth for each fjord, and the mean annual discharge for each fjord. \n",
    "\n",
    "## Requirements\n",
    "- `../freshwater/ice/{MAR, RACMO}.nc`\n",
    "- `../freshwater/land/{MAR, RACMO}.nc`\n",
    "- `./fjord_outlets`\n",
    "\n",
    "If you are generating your own data rather than reproducing our analysis, see instructions in README.\n",
    "\n",
    "## User parameters\n",
    "- data model source (MAR or RACMO)\n",
    "- list of fjord ID numbers\n",
    "- threshold depth for calculating discharge above/below threshold\n",
    "- altitude intervals for depth binning\n",
    "- start and end dates\n",
    "- colormap details\n",
    "\n",
    "## Output\n",
    "- Figures of freshwater discharge above/below a threshold depth in each fjord \\\n",
    "`../figures/fjord_FWdischarge_depththreshold_{MAR, RACMO}/fjord*_FWdischarge_threshold_{MAR, RACMO}.png`\n",
    "- Figures of freshwater discharge binned by depth across the whole depth range \\\n",
    "`../figures/fjord_FWdischarge_depthbinned_{MAR, RACMO}/fjord*_FWdischarge_depthbinned_{MAR, RACMO}.png`\n",
    "- Figure of freshwater discharge time series binned by depth and combined across all fjords, separated by land-sourced and ice-sourced discharge \\\n",
    "`../figures/allfjords_FWdischarge_depthbinned_{MAR, RACMO}_landvsice.png`\n",
    "- Figure of freshwater discharge time series binned by depth and combined across all fjords, land and ice sources summed together \\\n",
    "`../figures/allfjords_FWdischarge_depthbinned_{MAR, RACMO}_combined.png`\n",
    "- Spreadsheets of freshwater discharge time series for individual fjords, binned by depth and summed across land and ice sources \\\n",
    "`../databases/fjord_FWdischarge_depthbinned_{MAR, RACMO}_combined/fjord*_FWdischarge_depthbinned_{MAR, RACMO}.csv`\n",
    "- Spreadsheet of mean annual deep ice+land freshwater discharge for each fjord \\\n",
    "`../databases/fjord_FWdischarge_mean_annual_deep_{MAR, RACMO}.csv`\n",
    "\n",
    "## Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib.cm import ScalarMappable\n",
    "import numpy as np\n",
    "import copy\n",
    "import pathlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# USER PARAMETERS\n",
    "\n",
    "# -- Choose to load ice and land discharge from either RACMO or MAR\n",
    "source = 'MAR' # choose 'RACMO' or 'MAR'\n",
    "\n",
    "# -- Get list of fjord numbers\n",
    "fjord_number_list = list(range(1,53))\n",
    "\n",
    "# -- Set threshold depth for above/below discharge calculations (units=m)\n",
    "# -- Note: raw data convention uses negative for sub-zero \"altitude\"; we maintain this but reformat plot labels\n",
    "threshold_depth = -200\n",
    "\n",
    "# -- Define altitude intervals for depth binning (units=m)\n",
    "interval = 20\n",
    "top = 20\n",
    "shallow_bins = np.arange(threshold_depth, top+interval, interval)\n",
    "deep_bins = np.arange(-1000, threshold_depth+interval, interval)\n",
    "all_bins = np.arange(-1000, top+interval, interval)\n",
    "alt_bins = pd.interval_range(start=-1000, end=top, freq=interval)\n",
    "\n",
    "# -- Get list of dates for timespan of data (format='YYYY-MM-DD')\n",
    "time_start = pd.to_datetime('2015-01-01')\n",
    "time_end = pd.to_datetime('2019-12-31')\n",
    "time = pd.date_range(start=time_start, end=time_end)\n",
    "\n",
    "# -- Define Blues cmap with nodata color as dimgray, under_value color as silver\n",
    "cmap = copy.copy(plt.cm.get_cmap('Blues'))\n",
    "cmap.set_bad(color='dimgray')\n",
    "cmap.set_under(color='silver')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "discharge_ice = xr.open_dataset(f'../freshwater/ice/{source}.nc', engine='h5netcdf')\n",
    "discharge_land = xr.open_dataset(f'../freshwater/land/{source}.nc', engine='h5netcdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roundup(val, order=None):\n",
    "    '''Round value up to next given order of magnitude, used for setting colorbar limits'''\n",
    "    if order is None:\n",
    "        order = np.floor(np.log10(val))\n",
    "    val_roundup = np.ceil(val / 10**(order)) * 10**(order)\n",
    "    return val_roundup\n",
    "\n",
    "def major_formatter(x, pos):\n",
    "    '''Plot tick label formatting - switch negative altitudes to positive depths'''\n",
    "    label = str(-int(x)) if x < 0 else str(int(x))\n",
    "    return label\n",
    "\n",
    "def add_Gt_ticks_to_cbar(cbar):\n",
    "    '''For existing discharge colorbar in figure, add second set of ticks in Gt/yr'''\n",
    "    # move original cbar label to left side\n",
    "    cbar.ax.yaxis.set_label_position('left')\n",
    "    # copy cbar to create second cbar axes\n",
    "    cax2 = cbar.ax.twinx()\n",
    "    # set cax2 to log scale to match the m3/s scale\n",
    "    cax2.set_yscale('log')\n",
    "    # set cax2 lims to be equivalent to original cbar, converted from m3/s to Gt/yr\n",
    "    cax2.set_ylim([cbar.vmin * (60*60*24*365 / 1E9), cbar.vmax * (60*60*24*365 / 1E9)])\n",
    "    cax2.set_ylabel('Discharge ($\\mathrm{Gt\\ a^{-1}}$)')\n",
    "    cax2.yaxis.set_label_position('right')\n",
    "    return cax2\n",
    "\n",
    "def discharge_from_outlets(oi, ol):\n",
    "    # MODIFIED FROM KEN MANKOFF\n",
    "    ## oi & ol are outlets for ice and land respectively\n",
    "    \n",
    "    ## simple method, but fails if ol contains stations not in discharge_land\n",
    "    ## NOTE: not sure how that occurs\n",
    "    ## See https://github.com/pydata/xarray/issues/1990\n",
    "    # ri = discharge_ice.sel(station = oi)\n",
    "    # rl = discharge_land.sel(station = ol, method='nearest')\n",
    "\n",
    "    # update 2023-06-07: breaks if no stations present. Added if/else.\n",
    "    # If ice/land outlets are present, get discharge at them and sum\n",
    "    if not oi.empty:\n",
    "        ri = discharge_ice.reindex({'station':oi.id.values})\n",
    "        ri = ri.discharge.sum(dim='station', skipna=True).to_dataframe()\n",
    "    else:\n",
    "        ri = xr.full_like(discharge_ice.sel(station=slice(1,2)), fill_value=np.nan)\n",
    "        ri = ri.discharge.sum(dim='station', skipna=False).to_dataframe()\n",
    "    if not ol.empty:\n",
    "        rl = discharge_land.reindex({'station':ol.id.values})\n",
    "        rl = rl.discharge.sum(dim='station', skipna=True).to_dataframe()\n",
    "    else:\n",
    "        rl = xr.full_like(discharge_land.sel(station=slice(1,2)), fill_value=np.nan)\n",
    "        rl = rl.discharge.sum(dim='station', skipna=False).to_dataframe()\n",
    "\n",
    "    return ri, rl\n",
    "\n",
    "\n",
    "def discharge_from_outlets_binned(oi, ol, bins):\n",
    "    # bins can be an integer or a list (e.g. np.arange(-200, 250, 50))\n",
    "\n",
    "    # update 2023-06-07: breaks if no stations present. Added if/else.\n",
    "    # If ice/land outlets are present, get discharge at them\n",
    "    if not oi.empty:\n",
    "        ri = discharge_ice.reindex({'station':oi.id.values})\n",
    "    else:\n",
    "        ri = xr.full_like(discharge_ice.sel(station=slice(1,2)), fill_value=np.nan)\n",
    "        ri['alt'] = np.nan # force nan altitude for fake dataframe so binning works\n",
    "    if not ol.empty:\n",
    "        rl = discharge_land.reindex({'station':ol.id.values})\n",
    "    else:\n",
    "        rl = xr.full_like(discharge_land.sel(station=slice(1,2)), fill_value=np.nan)\n",
    "        rl['alt'] = np.nan # force nan altitude for fake dataframe so binning works\n",
    "\n",
    "    # bin discharge by altitude (depth)\n",
    "    if ri.alt.values.size == 0 or np.isnan(ri.alt.values).all() or (type(bins) is not int and (ri.alt.values > bins.max()).all()):\n",
    "        ri_bin = None\n",
    "    else: \n",
    "        ri.alt.values = np.clip(ri.alt.values, a_min=None, a_max=1) # clip all values above 0 so that they're in top bin, regardless of bin edge\n",
    "        ri_bin = ri.groupby_bins(group=ri.alt, bins=bins).sum()\n",
    "    if rl.alt.values.size == 0 or np.isnan(rl.alt.values).all() or (type(bins) is not int and (ri.alt.values > bins.max()).all()):\n",
    "        rl_bin = None\n",
    "    else: \n",
    "        rl.alt.values = np.clip(rl.alt.values, a_min=None, a_max=1) # clip all values above 0 so that they're in top bin, regardless of bin edge\n",
    "        rl_bin = rl.groupby_bins(group=rl.alt, bins=bins).sum()\n",
    "\n",
    "    return ri_bin, rl_bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example discharge-at-depth plot for a single fjord\n",
    "fjord_number = 33\n",
    "outlets = pd.read_csv(f'./fjord_outlets/fjord{fjord_number:02}_outlets.txt')\n",
    "# Separate land and ice outlets\n",
    "o_land = outlets.where(outlets.domain == 'land').dropna(subset=['domain'])\n",
    "o_ice = outlets.where(outlets.domain == 'ice').dropna(subset=['domain'])\n",
    "# Separate deep and shallow outlets (relative to a threshold depth)\n",
    "o_land_deep = o_land.where(o_land.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "o_land_shal = o_land.where(o_land.elev > threshold_depth).dropna(subset=['elev'])\n",
    "o_ice_deep = o_ice.where(o_ice.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "o_ice_shal = o_ice.where(o_ice.elev > threshold_depth).dropna(subset=['elev'])\n",
    "\n",
    "ris_bin, rls_bin = discharge_from_outlets_binned(o_ice_shal, o_land_shal, bins=shallow_bins)\n",
    "\n",
    "fig, ax = plt.subplots(layout='constrained')\n",
    "ris_bin.where(ris_bin.discharge != 0, other=1E-18).discharge.plot(\n",
    "    ax=ax, \n",
    "    cmap=cmap, \n",
    "    norm=colors.LogNorm(vmin=1, vmax=roundup(ris_bin.discharge.max().values), clip=False),\n",
    "    cbar_kwargs={'label': 'Discharge ($\\mathrm{m^3\\ s^{-1}}$)'}\n",
    ")\n",
    "start, end = ax.get_ylim()\n",
    "ax.yaxis.set_ticks(np.arange(start, end, interval))\n",
    "ax.yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "ax.set_ylabel('Depth (m)')\n",
    "ax.set_xlabel('Time')\n",
    "add_Gt_ticks_to_cbar(ax.collections[-1].colorbar)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate discharge plots for each fjord\n",
    "\n",
    "# -- Loop through each fjord and save plots of discharge above/below a threshold depth\n",
    "for fjord_number in fjord_number_list:\n",
    "    print(f'Processing Fjord {fjord_number}...')\n",
    "    outlets = pd.read_csv(f'./fjord_outlets/fjord{fjord_number:02}_outlets.txt')\n",
    "\n",
    "    # Separate land and ice outlets\n",
    "    o_land = outlets.where(outlets.domain == 'land').dropna(subset=['domain'])\n",
    "    o_ice = outlets.where(outlets.domain == 'ice').dropna(subset=['domain'])\n",
    "\n",
    "    # Separate deep and shallow outlets (relative to a threshold depth)\n",
    "    o_land_deep = o_land.where(o_land.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "    o_land_shal = o_land.where(o_land.elev > threshold_depth).dropna(subset=['elev'])\n",
    "    o_ice_deep = o_ice.where(o_ice.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "    o_ice_shal = o_ice.where(o_ice.elev > threshold_depth).dropna(subset=['elev'])\n",
    "\n",
    "    # THRESHOLD\n",
    "    # Get discharge from shallow/deep land/ice outlets\n",
    "    ris, rls = discharge_from_outlets(o_ice_shal, o_land_shal)\n",
    "    rid, rld = discharge_from_outlets(o_ice_deep, o_land_deep)\n",
    "    # Plot discharge above/below threshold depth for land and ice outlets\n",
    "    fig, axs = plt.subplots(nrows=1, ncols=2, sharex=True, figsize=(8,4), layout='constrained')\n",
    "    axs = axs.flatten()\n",
    "    fig.suptitle(f'Fjord {fjord_number}')\n",
    "    # -- Plot shallow discharge\n",
    "    axs[0].plot(ris, label='ice basins')\n",
    "    axs[0].plot(rls, label='land basins')\n",
    "    axs[0].legend(loc='upper left', bbox_to_anchor=(0.01, 0.99))\n",
    "    axs[0].set_title(f'Shallow FW discharge (above {-threshold_depth}m)')\n",
    "    axs[0].set_ylabel('Discharge ($\\mathrm{m^3\\ s^{-1}}$)')\n",
    "    axs[0].set_xlabel('Time')\n",
    "    axs[0].grid(True)\n",
    "    if (ris.discharge.isnull().all() & rls.discharge.isnull().all()):\n",
    "        axs[0].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center',\n",
    "            verticalalignment='center',\n",
    "            transform=axs[0].transAxes\n",
    "        )\n",
    "        axs[0].set_yticks([])\n",
    "        axs[0].grid(False)\n",
    "    # -- Plot deep discharge\n",
    "    axs[1].plot(rid, label='ice basins')\n",
    "    axs[1].plot(rld, label='land basins')\n",
    "    axs[1].set_title(f'Deep FW discharge (below {-threshold_depth}m)')\n",
    "    axs[1].set_xlabel('Time')\n",
    "    axs[1].grid(True)\n",
    "    if (rid.discharge.isnull().all() & rld.discharge.isnull().all()):\n",
    "        axs[1].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center',\n",
    "            verticalalignment='center',\n",
    "            transform=axs[1].transAxes\n",
    "        )\n",
    "        axs[1].set_yticks([])\n",
    "        axs[1].grid(False)\n",
    "    # -- Save output\n",
    "    pathlib.Path(f'../figures/fjord_FWdischarge_depththreshold_{source}/').mkdir(parents=True, exist_ok=True)\n",
    "    plt.savefig(\n",
    "        f'../figures/fjord_FWdischarge_depththreshold_{source}/fjord{fjord_number:02}_FWdischarge_threshold_{source}.png',\n",
    "        bbox_inches='tight',\n",
    "        dpi=300,\n",
    "        facecolor='white'\n",
    "    )\n",
    "    plt.close()\n",
    "\n",
    "    # BINNING\n",
    "    # Get discharge binned at depths for ice and land shallow and deep outlets\n",
    "    ris_bin, rls_bin = discharge_from_outlets_binned(o_ice_shal, o_land_shal, bins=shallow_bins)\n",
    "    rid_bin, rld_bin = discharge_from_outlets_binned(o_ice_deep, o_land_deep, bins=deep_bins)\n",
    "    max_discharge = max([da.discharge.max() for da in [ris_bin, rls_bin, rid_bin, rld_bin] if da is not None]).values\n",
    "    \n",
    "    # Plot discharge time series binned by depth\n",
    "    fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey='row', figsize=(8,7), layout='constrained')\n",
    "    axs = axes.flatten()\n",
    "    fig.suptitle(f'Fjord {fjord_number}')\n",
    "\n",
    "    # -- Plot shallow ice binned discharge\n",
    "    if ris_bin is not None:\n",
    "        ris_bin.where(ris_bin.discharge != 0, other=1E-18).discharge.plot(\n",
    "            ax=axs[0], \n",
    "            cmap=cmap,\n",
    "            norm=colors.LogNorm(vmin=1, vmax=roundup(max_discharge), clip=False),\n",
    "            add_colorbar=False\n",
    "        )\n",
    "    else:\n",
    "        axs[0].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center',\n",
    "            verticalalignment='center',\n",
    "            transform=axs[0].transAxes\n",
    "        )\n",
    "    axs[0].set_title('Shallow ice FW discharge')\n",
    "    axs[0].set_ylabel('Depth (m)')\n",
    "    axs[0].set_ylim(threshold_depth, interval)\n",
    "    axs[0].yaxis.set_ticks(np.arange(threshold_depth, 1, interval))\n",
    "    axs[0].yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "    axs[0].xaxis.label.set_visible(False)\n",
    "    \n",
    "    # -- Plot shallow land binned discharge\n",
    "    if rls_bin is not None:\n",
    "        rls_bin.where(rls_bin.discharge != 0, other=1E-18).discharge.plot(\n",
    "            ax=axs[1],\n",
    "            cmap=cmap,\n",
    "            norm=colors.LogNorm(vmin=1, vmax=roundup(max_discharge), clip=False),\n",
    "            add_colorbar=False\n",
    "        )\n",
    "    else:\n",
    "        axs[1].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center',\n",
    "            verticalalignment='center',\n",
    "            transform=axs[1].transAxes\n",
    "        )\n",
    "    axs[1].set_title('Shallow land FW discharge')\n",
    "    axs[1].yaxis.label.set_visible(False)\n",
    "    axs[1].xaxis.label.set_visible(False)\n",
    "\n",
    "    # -- Plot deep ice binned discharge\n",
    "    if rid_bin is not None:\n",
    "        rid_bin.where(rid_bin.discharge != 0, other=1E-18).discharge.plot(\n",
    "            ax=axs[2],\n",
    "            cmap=cmap,\n",
    "            norm=colors.LogNorm(vmin=1, vmax=roundup(max_discharge), clip=False),\n",
    "            add_colorbar=False\n",
    "        )\n",
    "    else:\n",
    "        axs[2].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center',\n",
    "            verticalalignment='center',\n",
    "            transform=axs[2].transAxes\n",
    "        )\n",
    "    axs[2].set_title('Deep ice FW discharge')\n",
    "    axs[2].set_ylabel('Depth (m)')\n",
    "    axs[2].set_xlabel('Time')\n",
    "    axs[2].set_ylim(-1000, threshold_depth)\n",
    "    axs[2].yaxis.set_ticks(np.arange(-1000, threshold_depth+1, 100))\n",
    "    axs[2].yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "    axs[2].xaxis.set_tick_params(reset=True, top=False)\n",
    "    axs[2].set_xlim(right=pd.to_datetime('2020-01-01'))\n",
    "    \n",
    "    # -- Plot deep land binned discharge\n",
    "    if rld_bin is not None:\n",
    "        rld_bin.where(rld_bin.discharge != 0, other=1E-18).discharge.plot(\n",
    "            ax=axs[3],\n",
    "            cmap=cmap,\n",
    "            norm=colors.LogNorm(vmin=1, vmax=roundup(max_discharge), clip=False),\n",
    "            add_colorbar=False\n",
    "        )\n",
    "    else:\n",
    "        axs[3].text(\n",
    "            x=0.5,\n",
    "            y=0.5,\n",
    "            s='NO DATA',\n",
    "            fontsize=20,\n",
    "            horizontalalignment='center', \n",
    "            verticalalignment='center', \n",
    "            transform=axs[3].transAxes\n",
    "        )\n",
    "    axs[3].set_title('Deep land FW discharge')\n",
    "    axs[3].set_xlabel('Time')\n",
    "    axs[3].yaxis.label.set_visible(False)\n",
    "    axs[3].xaxis.set_tick_params(reset=True, top=False)\n",
    "    axs[3].set_xlim(right=pd.to_datetime('2020-01-01'))\n",
    "\n",
    "    # -- Add colorbar for all subplots\n",
    "    cb = fig.colorbar(\n",
    "        ScalarMappable(\n",
    "            cmap=cmap,\n",
    "            norm=colors.LogNorm(vmin=1, vmax=roundup(max_discharge), clip=False)\n",
    "        ),\n",
    "        ax=axes.ravel().tolist(),\n",
    "        extend='min',\n",
    "        label='Discharge ($\\mathrm{m^3\\ s^{-1}}$)'\n",
    "    )\n",
    "    add_Gt_ticks_to_cbar(cb)\n",
    "    \n",
    "    # -- Save output\n",
    "    pathlib.Path(f'../figures/fjord_FWdischarge_depthbinned_{source}/').mkdir(parents=True, exist_ok=True)\n",
    "    plt.savefig(\n",
    "        f'../figures/fjord_FWdischarge_depthbinned_{source}/fjord{fjord_number:02}_FWdischarge_depthbinned_{source}.png',\n",
    "        bbox_inches='tight', \n",
    "        dpi=300, \n",
    "        facecolor='white'\n",
    "    )\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store ice and land discharge data in depth bins for each fjord\n",
    "\n",
    "# -- Initialize dataarrays for storage\n",
    "fjorddischarge_ice = xr.DataArray(\n",
    "    dims=['alt_bins', 'time', 'fjord'],\n",
    "    coords=dict(alt_bins=alt_bins, time=time, fjord=fjord_number_list)\n",
    ")\n",
    "fjorddischarge_land = xr.DataArray(\n",
    "    dims=['alt_bins', 'time', 'fjord'],\n",
    "    coords=dict(alt_bins=alt_bins, time=time, fjord=fjord_number_list)\n",
    ")\n",
    "\n",
    "# -- Loop through each fjord and save plots of discharge vs depth\n",
    "for fjord_number in fjord_number_list:\n",
    "    print(f'Processing Fjord {fjord_number}...')\n",
    "    outlets = pd.read_csv(f'./fjord_outlets/fjord{fjord_number:02}_outlets.txt')\n",
    "\n",
    "    # Separate land and ice outlets\n",
    "    o_land = outlets.where(outlets.domain == 'land').dropna(subset=['domain'])\n",
    "    o_ice = outlets.where(outlets.domain == 'ice').dropna(subset=['domain'])\n",
    "\n",
    "    # Get discharge from shallow/deep land/ice outlets\n",
    "    ri, rl = discharge_from_outlets_binned(o_ice, o_land, bins=all_bins)\n",
    "    \n",
    "    # Add discharge data to dataarrays\n",
    "    if ri is not None:\n",
    "        ri_da = xr.DataArray(ri.discharge, coords=dict(alt_bins=alt_bins, time=time, fjord=fjord_number))\n",
    "        fjorddischarge_ice.loc[dict(fjord=fjord_number)] = ri_da\n",
    "    if rl is not None:\n",
    "        rl_da = xr.DataArray(rl.discharge, coords=dict(alt_bins=alt_bins, time=time, fjord=fjord_number))\n",
    "        fjorddischarge_land.loc[dict(fjord=fjord_number)] = rl_da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot land and ice basin discharge time series binned by depth, for all fjords combined\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4), sharex=True, sharey=True, layout='constrained')\n",
    "axs = axes.flatten()\n",
    "max_sumdischarge = max([da.sum(dim='fjord').max().values for da in [fjorddischarge_ice, fjorddischarge_land] if da is not None])\n",
    "ceil_max_sumdischarge = roundup(max_sumdischarge.max(), order=4)\n",
    "\n",
    "# -- Plot binned ice discharge\n",
    "fjorddischarge_ice.where(fjorddischarge_ice != 0, other=1E-18).sum(dim='fjord').plot(\n",
    "    ax=axs[0], \n",
    "    cmap=cmap,\n",
    "    norm=colors.LogNorm(vmin=1, vmax=ceil_max_sumdischarge, clip=False),\n",
    "    add_colorbar=False\n",
    ")\n",
    "axs[0].set_title('FW discharge originating from ice basins')\n",
    "start, end = axs[0].get_ylim()\n",
    "axs[0].yaxis.set_ticks(np.arange(start, end, 100))\n",
    "axs[0].yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "axs[0].set_ylabel('Depth (m)')\n",
    "axs[0].set_xlabel('Time')\n",
    "axs[0].xaxis.set_tick_params(reset=True, top=False)\n",
    "axs[0].set_xlim(right=pd.to_datetime('2020-01-01'))\n",
    "    \n",
    "\n",
    "# -- Plot binned land discharge\n",
    "fjorddischarge_land.where(fjorddischarge_land != 0, other=1E-18).sum(dim='fjord').plot(\n",
    "    ax=axs[1],\n",
    "    cmap=cmap,\n",
    "    norm=colors.LogNorm(vmin=1, vmax=ceil_max_sumdischarge, clip=False),\n",
    "    add_colorbar=False\n",
    ")\n",
    "axs[1].set_title('FW discharge originating from land basins')\n",
    "start, end = axs[1].get_ylim()\n",
    "axs[1].yaxis.set_ticks(np.arange(start, end, 100))\n",
    "axs[1].yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "axs[1].yaxis.label.set_visible(False)\n",
    "axs[1].set_xlabel('Time')\n",
    "axs[1].xaxis.set_tick_params(reset=True, top=False)\n",
    "axs[1].set_xlim(right=pd.to_datetime('2020-01-01'))\n",
    "\n",
    "# -- Add one colorbar for both subplots\n",
    "cb = fig.colorbar(\n",
    "    ScalarMappable(\n",
    "        cmap=cmap,\n",
    "        norm=colors.LogNorm(vmin=1, vmax=ceil_max_sumdischarge, clip=False)\n",
    "    ),\n",
    "    ax=axes.ravel().tolist(),\n",
    "    extend='min',\n",
    "    label='Discharge ($\\mathrm{m^3\\ s^{-1}}$)'\n",
    ")\n",
    "add_Gt_ticks_to_cbar(cb)\n",
    "\n",
    "# -- Save output\n",
    "pathlib.Path('../figures/').mkdir(parents=True, exist_ok=True)\n",
    "plt.savefig(\n",
    "    f'../figures/allfjords_FWdischarge_depthbinned_{source}_landvsice.png',\n",
    "    bbox_inches='tight',\n",
    "    dpi=300,\n",
    "    facecolor='white'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot summed ice+land discharge time series binned by depth\n",
    "# -- concat along arbitrary axis helps to prevent nan issues when summing\n",
    "fjorddischarge_total = xr.concat((fjorddischarge_ice, fjorddischarge_land), 'z').sum(dim='z', min_count=1)\n",
    "ceil_totaldischarge = roundup(fjorddischarge_total.sum(dim='fjord').max().values.max(), 4)\n",
    "fig, ax = plt.subplots(figsize=(6, 4.5), layout='constrained')\n",
    "fjorddischarge_total.where(fjorddischarge_total != 0, other=1E-18).sum(dim='fjord', min_count=1).plot(\n",
    "    ax=ax, \n",
    "    cmap=cmap, \n",
    "    norm=colors.LogNorm(vmin=1, vmax=ceil_totaldischarge, clip=False),\n",
    "    cbar_kwargs={'label': 'Discharge ($\\mathrm{m^3\\ s^{-1}}$)'}\n",
    ")\n",
    "ax.set_title('Total FW discharge')\n",
    "start, end = ax.get_ylim()\n",
    "ax.yaxis.set_ticks(np.arange(start, end, 100))\n",
    "ax.yaxis.set_major_formatter(ticker.FuncFormatter(major_formatter))\n",
    "ax.set_ylabel('Depth (m)')\n",
    "ax.set_xlabel('Time')\n",
    "ax.xaxis.set_tick_params(reset=True, top=False)\n",
    "ax.set_xlim(right=pd.to_datetime('2020-01-01'))\n",
    "add_Gt_ticks_to_cbar(ax.collections[-1].colorbar)\n",
    "\n",
    "pathlib.Path('../figures/').mkdir(parents=True, exist_ok=True)\n",
    "plt.savefig(\n",
    "    f'../figures/allfjords_FWdischarge_depthbinned_{source}_combined.png',\n",
    "    bbox_inches='tight',\n",
    "    dpi=300,\n",
    "    facecolor='white'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save depth-binned freshwater flux time series for each fjord to CSV (round to 3 decimal places)\n",
    "pathlib.Path(f'../databases/fjord_FWdischarge_depthbinned_{source}_combined').mkdir(parents=True, exist_ok=True)\n",
    "fjorddischarge_total = fjorddischarge_total.round(decimals=3)\n",
    "for f in fjord_number_list:\n",
    "    fjorddischarge_total.sel(fjord=f).to_pandas().transpose().to_csv(\n",
    "        f'../databases/fjord_FWdischarge_depthbinned_{source}_combined/fjord{f:02}_FWdischarge_depthbinned_{source}.csv', \n",
    "        na_rep=np.nan,\n",
    "        header=True,\n",
    "        index=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get mean annual deep (ice+land) discharge for each fjord and save to CSV (round to 3 decimal places)\n",
    "\n",
    "# -- Initialize empty series to store mean annual deep discharge\n",
    "deep_discharge = pd.Series(index=fjord_number_list, name='discharge', dtype='float64')\n",
    "\n",
    "# -- Loop through each fjord and save plots of discharge vs depth\n",
    "for fjord_number in fjord_number_list:\n",
    "    print(f'Processing Fjord {fjord_number}...')\n",
    "    outlets = pd.read_csv(f'./fjord_outlets/fjord{fjord_number:02}_outlets.txt')\n",
    "\n",
    "    # Separate land and ice outlets\n",
    "    o_land = outlets.where(outlets.domain == 'land').dropna(subset=['domain'])\n",
    "    o_ice = outlets.where(outlets.domain == 'ice').dropna(subset=['domain'])\n",
    "\n",
    "    # Separate deep and shallow outlets (relative to a threshold depth)\n",
    "    o_land_deep = o_land.where(o_land.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "    o_ice_deep = o_ice.where(o_ice.elev <= threshold_depth).dropna(subset=['elev'])\n",
    "\n",
    "    # Get discharge from shallow/deep land/ice outlets\n",
    "    rid, rld = discharge_from_outlets(o_ice_deep, o_land_deep)\n",
    "\n",
    "    # Get mean annual combined (ice+land) discharge at depth\n",
    "    rd = rid.add(rld, fill_value=0) # convert nan to 0 in either df before adding\n",
    "    rd = rd * 86400 # convert from m3/s to m3/day\n",
    "    mean_annual_deep_discharge = rd.groupby(rd.index.year).sum().mean()\n",
    "    deep_discharge[fjord_number] = mean_annual_deep_discharge.round(decimals=3)\n",
    "\n",
    "pathlib.Path('../databases/').mkdir(parents=True, exist_ok=True)\n",
    "deep_discharge.to_csv(\n",
    "    f'../databases/fjord_FWdischarge_mean_annual_deep_{source}.csv',\n",
    "    index_label='fjord',\n",
    "    header=['mean annual discharge (m3)']\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}