{ "cells": [ { "cell_type": "markdown", "id": "f43c025c-6e0b-4394-b813-8c6c3738fe12", "metadata": {}, "source": [ "## Check the weird glaciers from https://github.com/OGGM/oggm/issues/1581" ] }, { "cell_type": "code", "execution_count": 3, "id": "453ca08e-1614-4766-84ae-24a593c68802", "metadata": {}, "outputs": [], "source": [ "import logging\n", "from datetime import timedelta\n", "import os\n", "import sys\n", "import numpy as np\n", "import xarray as xr\n", "\n", "import pandas as pd\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "\n", "# Import and verify\n", "import oggm\n", "\n", "\n", "import oggm.cfg as cfg\n", "from oggm import utils, tasks, entity_task\n", "import seaborn as sns\n", "from func_add_terrafirma import SCENARIO_PARENTS" ] }, { "cell_type": "code", "execution_count": 4, "id": "9895f6c5-4935-4716-94d3-8d7e3aeaf6b2", "metadata": {}, "outputs": [], "source": [ "weird_glaciers_dict = rgi_by_region = {\n", " \"01\": [\n", " \"RGI60-01.21011\",\n", " ],\n", " \"03\": [\n", " \"RGI60-03.01053\",\n", " \"RGI60-03.01538\",\n", " \"RGI60-03.01972\",\n", " \"RGI60-03.02785\",\n", " \"RGI60-03.03419\",\n", " \"RGI60-03.03585\",\n", " \"RGI60-03.03980\",\n", " ],\n", " \"04\": [\n", " \"RGI60-04.00190\",\n", " \"RGI60-04.00318\",\n", " \"RGI60-04.02794\",\n", " \"RGI60-04.02942\",\n", " \"RGI60-04.03275\",\n", " \"RGI60-04.04920\",\n", " \"RGI60-04.05013\",\n", " \"RGI60-04.06461\",\n", " ],\n", " \"05\": [\n", " \"RGI60-05.02477\",\n", " \"RGI60-05.08792\",\n", " \"RGI60-05.09561\",\n", " \"RGI60-05.09696\",\n", " \"RGI60-05.10000\",\n", " \"RGI60-05.10033\",\n", " \"RGI60-05.10149\",\n", " \"RGI60-05.10204\",\n", " \"RGI60-05.10521\",\n", " ],\n", " \"07\": [\n", " \"RGI60-07.00079\",\n", " \"RGI60-07.00131\",\n", " \"RGI60-07.00226\",\n", " \"RGI60-07.00284\",\n", " \"RGI60-07.01136\",\n", " \"RGI60-07.01160\",\n", " \"RGI60-07.01569\",\n", " \"RGI60-07.01570\",\n", " ],\n", " \"08\": [\n", " \"RGI60-08.00005\",\n", " \"RGI60-08.00287\",\n", " \"RGI60-08.02392\",\n", " ],\n", " \"09\": [\n", " \"RGI60-09.00017\",\n", " \"RGI60-09.00244\",\n", " \"RGI60-09.00461\",\n", " \"RGI60-09.00477\",\n", " \"RGI60-09.00546\",\n", " \"RGI60-09.00564\",\n", " \"RGI60-09.00565\",\n", " \"RGI60-09.00922\",\n", " \"RGI60-09.00930\",\n", " \"RGI60-09.00933\",\n", " ],\n", " \"10\": [\n", " \"RGI60-10.04609\",\n", " \"RGI60-10.05039\",\n", " ],\n", " \"13\": [\n", " \"RGI60-13.53885\",\n", " ],\n", " \"17\": [\n", " \"RGI60-17.02945\",\n", " ],\n", "}" ] }, { "cell_type": "code", "execution_count": 39, "id": "93dd218f-1133-4c11-aff3-66ed573d5dfd", "metadata": {}, "outputs": [], "source": [ "load = False\n", "if load:\n", " ds_reg_l = []\n", " for rgi_reg in weird_glaciers_dict.keys():\n", " print(rgi_reg)\n", " rgi_ids = weird_glaciers_dict[rgi_reg]\n", " p = f'/home/www/lschuster/terrafirma_oggm_proj/output_dir/RGI{rgi_reg}/' \n", " dict_ds = []\n", " for scenario in SCENARIO_PARENTS.keys():\n", " if 'up2p0' in scenario:\n", " path = f'{p}run_terrafirma_UKESM1-2-LL_esm_{scenario}_bc_1975_2014_to_0_40_sim_year_starting_at_1979_glacier_state_from_0_sim_year_climate_Batch*.nc'\n", " with xr.open_mfdataset(path) as _d:\n", " _d = _d.volume\n", " _d = _d.sel(rgi_id=rgi_ids)\n", " \n", " _d['scenario'] = scenario\n", " dict_ds.append(_d)\n", " ds_reg = xr.concat(dict_ds, dim='scenario')\n", " ds_reg_l.append(ds_reg)\n", " ds_weird_glaciers = xr.concat(ds_reg_l, dim='rgi_id')\n", " ds_weird_glaciers = ds_weird_glaciers.load()\n", " ds_weird_glaciers.to_netcdf('aggregated_weird_glaciers.nc')" ] }, { "cell_type": "code", "execution_count": null, "id": "994c3f31-09d4-42f3-9438-d6d43d9c3399", "metadata": {}, "outputs": [], "source": [ "ds_weird_glaciers = xr.open_dataset('aggregated_weird_glaciers.nc')" ] }, { "cell_type": "code", "execution_count": 41, "id": "a567d4f3-5b39-4444-847f-d253cb7059c8", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'volume' (scenario: 31, time: 881, rgi_id: 50)> Size: 5MB\n",
       "array([[[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n",
       "         4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n",
       "        [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n",
       "         4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n",
       "        [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n",
       "         4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n",
       "        ...,\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan]],\n",
       "\n",
       "       [[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n",
       "         4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n",
       "        [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n",
       "         4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n",
       "        [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n",
       "         4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n",
       "...\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan]],\n",
       "\n",
       "       [[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n",
       "         4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n",
       "        [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n",
       "         4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n",
       "        [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n",
       "         4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n",
       "        ...,\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan],\n",
       "        [          nan,           nan,           nan, ...,\n",
       "                   nan,           nan,           nan]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * scenario        (scenario) <U23 3kB 'up2p0' ... 'up2p0-gwl6p0-200y-dn2p0'\n",
       "  * time            (time) float64 7kB 1.975e+03 1.976e+03 ... 2.855e+03\n",
       "  * rgi_id          (rgi_id) <U14 3kB 'RGI60-01.21011' ... 'RGI60-17.02945'\n",
       "    hydro_year      (scenario, time) float64 218kB 1.975e+03 1.976e+03 ... nan\n",
       "    hydro_month     (rgi_id, scenario, time) float64 11MB 4.0 4.0 ... nan nan\n",
       "    calendar_year   (scenario, time) float64 218kB 1.975e+03 1.976e+03 ... nan\n",
       "    calendar_month  (scenario, time) float64 218kB 1.0 1.0 1.0 ... nan nan nan\n",
       "Attributes:\n",
       "    description:  Total glacier volume\n",
       "    unit:         m 3
" ], "text/plain": [ " Size: 5MB\n", "array([[[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n", " 4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n", " [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n", " 4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n", " [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n", " 4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n", " ...,\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan]],\n", "\n", " [[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n", " 4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n", " [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n", " 4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n", " [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n", " 4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n", "...\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan]],\n", "\n", " [[6.2354342e+10, 3.3177219e+10, 1.7273092e+10, ...,\n", " 4.0509935e+06, 7.2912374e+09, 5.2739789e+09],\n", " [6.2425416e+10, 3.3163870e+10, 1.7288436e+10, ...,\n", " 4.0645650e+06, 7.3000556e+09, 5.2936443e+09],\n", " [6.2531195e+10, 3.3140271e+10, 1.7325959e+10, ...,\n", " 4.0834302e+06, 7.3317811e+09, 5.2856315e+09],\n", " ...,\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan],\n", " [ nan, nan, nan, ...,\n", " nan, nan, nan]]], dtype=float32)\n", "Coordinates:\n", " * scenario (scenario) " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\n", "\n", "# ------------------------------------------------------------\n", "# 1. Check that initial volume is equal across all scenarios\n", "# ------------------------------------------------------------\n", "\n", "\n", "initial_by_scenario = ds_weird_glaciers.isel(time=0)\n", "initial_reference = initial_by_scenario.isel(scenario=0)\n", "\n", "difference = abs(initial_by_scenario - initial_reference)\n", "\n", "max_difference = difference.max(\n", " dim=\"scenario\",\n", " skipna=True,\n", ")\n", "\n", "assert bool((max_difference < 1e-6).all()), (\n", " \"Initial glacier volume differs between scenarios. \"\n", " f\"Maximum difference: {difference.max().item()}\"\n", ")\n", "\n", "print(\"Check passed: initial volumes are equal across scenarios.\")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 2. Sort glaciers by initial volume, largest first\n", "# ------------------------------------------------------------\n", "\n", "sort_indices = np.argsort(initial_reference.values)[::-1]\n", "\n", "ds_sorted = ds_weird_glaciers.isel(\n", " rgi_id=sort_indices\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 3. Create the multipanel figure\n", "# ------------------------------------------------------------\n", "\n", "n_glaciers = ds_sorted.sizes[\"rgi_id\"]\n", "\n", "ncols = 5\n", "nrows = int(np.ceil(n_glaciers / ncols))\n", "\n", "fig, axes = plt.subplots(\n", " nrows=nrows,\n", " ncols=ncols,\n", " figsize=(25, 3.0 * nrows),\n", " sharex=True,\n", " sharey=False,\n", " squeeze=False,\n", ")\n", "\n", "axes = axes.ravel()\n", "\n", "legend_handles = []\n", "legend_labels = []\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 4. Plot each glacier\n", "# ------------------------------------------------------------\n", "\n", "for glacier_idx in range(n_glaciers):\n", "\n", " ax = axes[glacier_idx]\n", " ds_glacier = ds_sorted.isel(rgi_id=glacier_idx)\n", "\n", " for scenario in ds_glacier.scenario.values:\n", "\n", " ds_scenario = ds_glacier.sel(scenario=scenario)\n", " style = get_scenario_style(scenario)\n", "\n", " line, = ax.plot(\n", " ds_scenario.time,\n", " ds_scenario.values / 1e9,\n", " color=style[\"color\"],\n", " linestyle=style[\"linestyle\"],\n", " linewidth=style[\"linewidth\"] + 0.5,\n", " label=str(scenario),\n", " )\n", "\n", " # Collect legend entries only once\n", " if glacier_idx == 0:\n", " legend_handles.append(line)\n", " legend_labels.append(str(scenario))\n", "\n", " rgi_id = ds_glacier.rgi_id.item()\n", "\n", " initial_volume_km3 = (\n", " ds_glacier\n", " .isel(time=0, scenario=0)\n", " .item()\n", " / 1e9\n", " )\n", "\n", " ax.set_title(\n", " f\"{rgi_id}\\n\"\n", " f\"$V_0$ = {initial_volume_km3:.3f} km³\",\n", " fontsize=12,\n", " pad=4,\n", " )\n", "\n", " ax.tick_params(\n", " axis=\"both\",\n", " labelsize=12,\n", " )\n", "\n", " ax.grid(\n", " alpha=0.2,\n", " linewidth=0.5,\n", " )\n", "\n", "\n", "# Hide unused axes\n", "for ax in axes[n_glaciers:]:\n", " ax.set_visible(False)\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 5. Shared labels\n", "# ------------------------------------------------------------\n", "\n", "fig.supxlabel(\n", " \"Year\",\n", " fontsize=17,\n", " y=0.015,\n", ")\n", "\n", "fig.supylabel(\n", " \"Volume (km³)\",\n", " fontsize=17,\n", " x=0.02,\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 6. Figure-level legend\n", "# ------------------------------------------------------------\n", "\n", "legend_ncols = min(6, len(legend_labels))\n", "legend_nrows = int(np.ceil(len(legend_labels) / legend_ncols))\n", "\n", "fig.legend(\n", " handles=legend_handles,\n", " labels=legend_labels,\n", " loc=\"upper center\",\n", " bbox_to_anchor=(0.5, 0.99),\n", " ncol=legend_ncols,\n", " fontsize=13,\n", " frameon=False,\n", " columnspacing=1.4,\n", " handlelength=2.5,\n", " handletextpad=0.6,\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 7. Explicitly reserve space for the legend\n", "# ------------------------------------------------------------\n", "\n", "# More legend rows -> more reserved space\n", "subplot_top = 0.99 - 0.015 * legend_nrows\n", "\n", "fig.subplots_adjust(\n", " left=0.055,\n", " right=0.985,\n", " bottom=0.045,\n", " top=subplot_top,\n", " hspace=0.60,\n", " wspace=0.25,\n", ")\n", "\n", "plt.show()\n", "\n", "# Optional:\n", "fig.savefig(\n", " \"weird_glaciers_sorted_by_initial_volume.png\",\n", " dpi=200,\n", " bbox_inches=\"tight\",\n", " )" ] }, { "cell_type": "markdown", "id": "bc2a35a8-9556-485e-affd-afe31bcfb49b", "metadata": {}, "source": [ "it seems like it is fine for all previously \"weird-defined\" glaciers, specifically for the larger ones.. but there are a few interesting ones... like RGI60-03-03585... for these glaciers the scenarios stop earlier as they would grow out of the boundary" ] }, { "cell_type": "code", "execution_count": null, "id": "552546f9-eb16-45c7-bb9f-981f2eeb6e56", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:oggm_env_2025]", "language": "python", "name": "conda-env-oggm_env_2025-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.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }