{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cc496cb6-c99c-4665-b6f1-516eb4b87398",
   "metadata": {},
   "source": [
    "## Aggregate volume and area per-glacier data for the per-glacier deglaciation analysis of Lander Van Tricht"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "07cecf4e-5699-43cf-9906-130e310fcfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from oggm.utils import mkdir\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "6d6af737-bf6c-4b41-91bc-733b5c036ee5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we only aggregate here `w5e5_gcm_merged`, but you could repeat this here with `gcm_from_2000` as it is also available in the raw per-glacier datafiles\n",
    "historical_future_option = 'w5e5_gcm_merged'\n",
    "# this is the only preprocessed version (i.e., we match the geodetic observations on any individual glacier)\n",
    "exp = 'match_geod_pergla'\n",
    "oggm_version = 'oggm_v1.6.1_2023.3'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "7ee93261-c6c1-4bd4-8aed-90698b08252c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dirpath = '/home/www/oggm/oggm-standard-projections/oggm_v16/2023.3/CMIP6/2100/'\n",
    "rgi_reg_list = [f\"RGI{str(i).zfill(2)}\" for i in range(1, 20)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "2c91a8ae-dc13-420e-aac3-357a420512e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "GCMs = [\"BCC-CSM2-MR\", #\n",
    "        \"CAMS-CSM1-0\", #\n",
    "        \"CESM2\", # \n",
    "        \"CESM2-WACCM\",#\n",
    "        \"EC-Earth3\",#\n",
    "        \"EC-Earth3-Veg\", #\n",
    "        \"FGOALS-f3-L\",#\n",
    "        \"GFDL-ESM4\", #\n",
    "        \"INM-CM4-8\", #\n",
    "        \"INM-CM5-0\", #\n",
    "        \"MPI-ESM1-2-HR\",#\n",
    "        \"MRI-ESM2-0\",#\n",
    "        \"NorESM2-MM\"]\n",
    "scenarios = ['ssp126','ssp245','ssp370','ssp585']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "3aad4f85-1e5b-4765-a121-18f67d088c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "aggregate = False\n",
    "if aggregate: \n",
    "    for GCM in GCMs:\n",
    "        for rgi_reg in rgi_reg_list:\n",
    "            for scen in scenarios: \n",
    "                dirpath_reg = dirpath + rgi_reg\n",
    "                fpath = f'{dirpath_reg}/run_hydro_{historical_future_option}_{GCM}_{scen}_bc_2000_2019_*.nc'\n",
    "                with xr.open_mfdataset(fpath) as ds:\n",
    "                    ds = ds[['volume', 'area']]\n",
    "                    ds = ds.load()\n",
    "                    \n",
    "                    pd_vol = ds['volume'].to_dataframe()[['volume']].reset_index()\n",
    "                    pd_vol.time = pd_vol.time.astype(int)\n",
    "                    pd_vol['ID'] = pd_vol.rgi_id.str[-5:].astype(int).tolist()\n",
    "                    pd_vol = pd_vol.pivot(index='ID', columns = 'time', values='volume')\n",
    "                \n",
    "                    pd_area = ds['area'].to_dataframe()[['area']].reset_index()\n",
    "                    pd_area.time = pd_area.time.astype(int)\n",
    "                    pd_area['ID'] = pd_area.rgi_id.str[-5:].astype(int).tolist()\n",
    "                    pd_area = pd_area.pivot(index='ID', columns = 'time', values='area')\n",
    "                \n",
    "                    pd_vol.to_csv(f'{GCM}/{rgi_reg}_volume_{GCM}_{scen}.csv')\n",
    "                    pd_area.to_csv(f'{GCM}/{rgi_reg}_area_{GCM}_{scen}.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "fd7f9c1b-e85b-4cc6-aa25-bd6da9044a5f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "da9b0ac5-0cb3-4e77-9dc8-b7ab505f21e1",
   "metadata": {},
   "source": [
    "## Test if sum of per-glacier data gives regional data \n",
    "--> I only select the common running glaciers because the aggregated regional data used within Zekollari et al. (2024) only includes those common-running glacier data ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "76fad22b-77e3-4b3d-a04f-50cfa2da4d79",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "_base = '/home/www/oggm/oggm-standard-projections/oggm-standard-projections-csv-files/1.6.1/common_running_2100/'\n",
    "with open(f'{_base}/rgi_ids_missing.json', 'r') as f:\n",
    "    invalid_per_reg = json.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "eda4c319-1ffd-4cd5-8b19-912e0c7d94ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that correct amount of glaciers are included (some with NaNs)\n",
    "rgi_meta = pd.read_hdf('/home/www/oggm/rgi/rgi62_stats.h5')\n",
    "rgi_meta = rgi_meta.loc[rgi_meta.Connect != 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "064145b0-fad3-4623-bb47-27de17b624ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no error, so the data seems to be the same!\n"
     ]
    }
   ],
   "source": [
    "_reg_file = '/home/www/oggm/oggm-standard-projections/oggm-standard-projections-csv-files/1.6.1/common_running_2100/'\n",
    "for GCM in GCMs:\n",
    "    for rgi_reg in rgi_reg_list:\n",
    "        for scen in scenarios: \n",
    "            n_glac = len(rgi_meta.loc[rgi_meta['O1Region'] == rgi_reg[3:]])\n",
    "            pd_vol = pd.read_csv(f'{GCM}/{rgi_reg}_volume_{GCM}_{scen}.csv', index_col=0)\n",
    "            pd_area = pd.read_csv(f'{GCM}/{rgi_reg}_area_{GCM}_{scen}.csv', index_col=0)\n",
    "            reg_file_vol = _reg_file + f'volume/CMIP6/2100/{rgi_reg}/{scen}.csv'\n",
    "            reg_file_area = _reg_file + f'area/CMIP6/2100/{rgi_reg}/{scen}.csv'\n",
    "            pd_reg_vol = pd.read_csv(reg_file_vol, index_col=0)\n",
    "            pd_reg_area = pd.read_csv(reg_file_area, index_col=0)\n",
    "            # all glaciers included? \n",
    "            assert len(pd_vol) == n_glac\n",
    "            assert len(pd_area) == n_glac\n",
    "            # now select only non-failing glaciers to compare it to the regional estimate\n",
    "            invalid_per_reg_numbers = [int(item.split(f'RGI60-{rgi_reg[3:]}.')[-1]) \n",
    "                                       for item in invalid_per_reg[rgi_reg]]\n",
    "            pd_area = pd_area.loc[~pd_area.index.isin(invalid_per_reg_numbers)]\n",
    "            pd_vol = pd_vol.loc[~pd_vol.index.isin(invalid_per_reg_numbers)]\n",
    "\n",
    "            try:\n",
    "                np.testing.assert_allclose(pd_reg_vol[GCM.upper()].loc[np.arange(2000,2101,1)],\n",
    "                                       pd_vol.T.sum(axis=1), rtol=0.0001)\n",
    "            except:\n",
    "                np.testing.assert_allclose(pd_reg_vol[GCM.upper()].loc[np.arange(2000,2101,1)],\n",
    "                       pd_vol.T.sum(axis=1), rtol=0.0005)\n",
    "                print('volume', GCM, rgi_reg, scen)\n",
    "            #######\n",
    "            try:\n",
    "                np.testing.assert_allclose(pd_reg_area[GCM.upper()].loc[np.arange(2000,2101,1)],\n",
    "                                       pd_area.T.sum(axis=1), rtol=0.0001)\n",
    "            except:\n",
    "                np.testing.assert_allclose(pd_reg_area[GCM.upper()].loc[np.arange(2000,2101,1)],\n",
    "                       pd_area.T.sum(axis=1), rtol=0.0005)\n",
    "                print('area', GCM, rgi_reg, scen)\n",
    "print('no error, so the data seems to be the same!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b68ea629-91c4-4372-b8d4-8dfca7fe1aef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:oggm_gmip3_working]",
   "language": "python",
   "name": "conda-env-oggm_gmip3_working-py"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}