{ "cells": [ { "cell_type": "markdown", "id": "b2d05601", "metadata": {}, "source": [ "**Sep2026: runs without errors**\n", "\n", "# Notebook to show the steady-state Glacier model differences " ] }, { "cell_type": "markdown", "id": "084b733f-0709-4aef-ad71-5d0abae71dab", "metadata": {}, "source": [ "- same approach but with OGGM_v153 is in the \"neighbouring\" jupyter notebook..." ] }, { "cell_type": "code", "execution_count": 1, "id": "f3310431-a48e-426f-9137-2e25bda4badc", "metadata": {}, "outputs": [], "source": [ "# # these code lines above are equal to those at the beginning of 2_glacier_vs_climate_change\n", "\n", "# download it here https://cluster.klima.uni-bremen.de/~lschuster/glacierMIP3_analysis/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled.nc\n", "# and change the path to your local path\n", "DATE = 'Feb12_2024' #'Jun29' #'apr04' (has OGGM runaway effect of glaciers inside)\n", "fill_option = 'repeat_last_101yrs'\n", "\n", "import xarray as xr\n", "import numpy as np\n", "import pandas as pd\n", "import scipy\n", "import os\n", "import glob\n", "import matplotlib.pyplot as plt\n", "from datetime import date\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from matplotlib.colors import LinearSegmentedColormap\n", "from sklearn.preprocessing import MinMaxScaler\n", "import matplotlib\n", "matplotlib.__version__\n", "from matplotlib.lines import Line2D\n", "import matplotlib.patches as mpatches\n", "from matplotlib.ticker import MultipleLocator\n", "import matplotlib.colors as mcolors\n", "\n", "# for plots within plots \n", "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n", "\n", "\n", "from help_functions import (pal_models, model_order, d_reg_num_name, model_order_anonymous, compute_steady_state_yr,\n", " get_glob_temp_exp)\n", "\n", "\n", "path_merged_runs_scaled_extend = f'../data/GMIP3_reg_glacier_model_data/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}.nc'\n", "_ds_reg_models_extend = xr.open_dataset(path_merged_runs_scaled_extend)" ] }, { "cell_type": "code", "execution_count": 2, "id": "7be43c67-027f-4335-a13c-669407e78fee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['PyGEM-OGGM_v13', 'GloGEMflow', 'GloGEMflow3D', 'OGGM_v16', 'GLIMB', 'Kraaijenbrink', 'GO', 'CISM2'] ['model 1', 'model 2', 'model 3', 'model 4', 'model 5', 'model 6', 'model 7', 'model 8']\n" ] } ], "source": [ "\n", "# --> comes from isimip3b_postprocessing_analysis/isimip3b_postprocess_to_monthly.ipynb\n", "pd_global_temp_exp = pd.read_csv('../data/climate_input_data/temp_ch_ipcc_ar6_isimip3b.csv', index_col = 0)\n", "\n", "hue_order_anonymous = []\n", "\n", "pal_models_l = []\n", "hue_order = []\n", "for m, p in zip(model_order, pal_models):\n", " if (m!='OGGM-VAS') and (m!='OGGM_v153'):\n", " hue_order.append(m)\n", " pal_models_l.append(p)\n", "for m in hue_order:\n", " hue_order_anonymous.append(model_order_anonymous[m])\n", "pal_models = pal_models_l\n", "\n", "print(hue_order, hue_order_anonymous)\n", "# select the right models:\n", "pal_models = sns.color_palette(pal_models)\n", "\n", "dict_model_col = {}\n", "for c,m in zip(pal_models, hue_order):\n", " dict_model_col[m] = c\n", " \n", "# select the right models\n", "\n", "glac_models = hue_order\n", "\n", "num_dict = {0:'(a)', 1:'(b)', 2:'(c)', 3:'(d)', 4: '(e)', 5:'(f)', 6:'(g)', 7:'(h)', 8:'(i)', 9:'(j)', 10:'(k)', 11:'(l)', 12:'(m)'} \n", "\n", "ind_t = pd_global_temp_exp.groupby(['gcm','period_scenario']).mean()['temp_ch_ipcc'].index\n", "np.testing.assert_allclose(get_glob_temp_exp(region='global')['temp_ch_ipcc'].loc[ind_t].values,\n", " pd_global_temp_exp.groupby(['gcm','period_scenario']).mean()['temp_ch_ipcc'].values, rtol=1e-5)\n", "\n", "pd_global_temp_exp_m = get_glob_temp_exp(region='global')" ] }, { "cell_type": "code", "execution_count": 3, "id": "692a5921-20f3-4631-9503-77556a76c403", "metadata": {}, "outputs": [], "source": [ "\n", "glob_temp_ch = get_glob_temp_exp(region='global')['temp_ch_ipcc'].values\n", "\n", "only_12_range = False\n", "if only_12_range:\n", " colors_icci_l3 = ['#70B8FF', '#FABB00', '#F25100', '#D42300', '#B3001E', '#800040', '#550066', '#2B084D']\n", " cmap_icci_3 = LinearSegmentedColormap.from_list('bins_100', colors_icci_l3, N=1000)\n", " # IDEA --> Let's scale from 1.2 onwards\n", " scaler = MinMaxScaler()\n", " temp_ch_sel = np.arange(1.2,glob_temp_ch.max(),0.25)\n", " scaler.fit(temp_ch_sel.reshape(-1,1))\n", " cmap = cmap_icci_3\n", "else:\n", " ## '#306BAF', '#1A488A', '#002966'\n", " ## '#4D8ED6',\n", " colors_full_temp_range = ['#002966', '#306BAF', '#70B8FF','#FABB00', '#F25100', '#D42300', '#B3001E', '#800040', '#550066', '#2B084D']\n", " cmap_full_temp_range = LinearSegmentedColormap.from_list('bins_100', colors_full_temp_range, N=1000)\n", " # IDEA --> Let's scale from the entire range onwards onwards\n", " scaler = MinMaxScaler()\n", " temp_ch_sel = np.arange(-0.36,glob_temp_ch.max(),0.25)# need to manually set the minimum to a lower value to have a light blue color at ~1.2°C\n", " scaler.fit(temp_ch_sel.reshape(-1,1))\n", " cmap = cmap_full_temp_range\n", "\n", "norm = plt.Normalize(vmin=-0.36, vmax=glob_temp_ch.max()) # we need the -0.4 for the correct norm \n", "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n", "#sm.set_clim(glob_temp_ch.min(),glob_temp_ch.max())\n", "# ok, both approaches work the same \n", "assert sm.cmap == cmap\n", "\n", "t = 1.2\n", "# just want to make sure that 1.2°C is this color here: '#70B8FF'\n", "hex_color = '#70B8FF'\n", "rgba_color = mcolors.to_rgba(hex_color)\n", "c1=sm.cmap(scaler.transform(np.array(t).reshape(-1,1))).squeeze()\n", "c2=cmap(scaler.transform(np.array(t).reshape(-1,1))).squeeze()\n", "np.testing.assert_allclose(c1,c2)\n", "np.testing.assert_allclose(c1,rgba_color)" ] }, { "cell_type": "code", "execution_count": 4, "id": "1efe39cb-8194-4e58-adbd-3b44e83ccb67", "metadata": {}, "outputs": [], "source": [ "_count_models = _ds_reg_models_extend.isel(simulation_year=0).to_dataframe().dropna().groupby(['model_author']).count()\n", "global_models = list(_count_models.loc[_count_models.volume_m3==1520].index)" ] }, { "cell_type": "code", "execution_count": 5, "id": "712b25df-776b-41ee-9a2f-1643bf03384f", "metadata": {}, "outputs": [], "source": [ "# select the right models and option\n", "approach = '_via_5yravg'\n", "### todo: update paths to include wOGGM_v153 or exclude OGGM_v153\n", "_ds_reg_yr_shift= xr.open_dataset(f'../data/GMIP3_reg_glacier_model_data/all_shifted_glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}{approach}.nc') #_wOGGM_v153.nc')\n", "shift_years = True\n", "if shift_years:\n", " add = '_shifted_to_2020'\n", " # we shift by maximum by +50 years, so like that there should always be values inside \n", " ds_reg_yr_shift = _ds_reg_yr_shift.sel(year_after_2020=slice(0,4950))\n", " ds_reg_models = ds_reg_yr_shift.sel(model_author=glac_models)\n", " ds_reg_models_vol = ds_reg_models.volume_m3 #.isel(simulation_year=slice(0,2000)) #['relative volume change (in %)']\n", "else:\n", " add = '_no_shift'\n", " ds_reg_models = _ds_reg_models_extend.sel(model_author=glac_models)\n", " ds_reg_models_vol = ds_reg_models.volume_m3 #.isel(simulation_year=slice(0,2000)) #['relative volume change (in %)']\n", "\n", "# Let's take the median estimate from all glacier models for every RGI region and then do the sum:\n", "ds_reg_models_vol = ds_reg_models_vol.stack(experiments=['gcm','period_scenario'])\n", "\n", "# median estimae for every region\n", "ds_reg_models_med_vol = ds_reg_models_vol.median(dim='model_author')\n", "\n", "## now: the same for global: IMPORTANT - first the sum over the regions and then the median over the glacier models\n", "ds_reg_models_vol_glob_med = ds_reg_models_vol.sel(model_author=global_models).sum(dim='rgi_reg').median(dim='model_author')\n" ] }, { "cell_type": "markdown", "id": "d07069ce", "metadata": {}, "source": [ "----" ] }, { "cell_type": "markdown", "id": "b5a21a4b", "metadata": {}, "source": [ "- todo: maybe use regional glacier stuff instead for the \"classes\" ???" ] }, { "cell_type": "code", "execution_count": 6, "id": "8277c45d", "metadata": {}, "outputs": [], "source": [ "ds_reg_models_med_vol = ds_reg_models_med_vol.assign_coords(temp_ch_ipcc= ('experiments',\n", " pd_global_temp_exp_m.loc[ds_reg_models_med_vol.experiments, # make sure thatwe use the same experiments\n", " 'temp_ch_ipcc']))\n", "# need that to have bot volume_m3 and temp_ch_ipcc as variables\n", "ds_reg_models_med_vol = ds_reg_models_med_vol.reset_coords()\n", "\n", "ds_reg_models_vol = ds_reg_models_vol.assign_coords(temp_ch_ipcc = ('experiments',\n", " pd_global_temp_exp_m.loc[ds_reg_models_vol.experiments, # make sure thatwe use the same experiments\n", " 'temp_ch_ipcc']))\n", "# need that to have bot volume_m3 and temp_ch_ipcc as variables\n", "ds_reg_models_vol = ds_reg_models_vol.reset_coords()\n", "############\n", "ds_reg_models_vol_glob = ds_reg_models_vol.sel(model_author=global_models).sum(dim='rgi_reg')\n", "\n", "\n", "### same for global median \n", "ds_reg_models_vol_glob_med = ds_reg_models_vol_glob_med.assign_coords(temp_ch_ipcc = ('experiments',\n", " pd_global_temp_exp_m.loc[ds_reg_models_vol_glob_med.experiments, # make sure thatwe use the same experiments\n", " 'temp_ch_ipcc']))\n", "ds_reg_models_vol_glob_med = ds_reg_models_vol_glob_med.reset_coords()\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "f7558dff-5076-4b6b-9dfa-01c80a872e03", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset>\n",
"Dimensions: (model_author: 8, year_after_2020: 4951, rgi_reg: 19,\n",
" experiments: 80)\n",
"Coordinates:\n",
" * model_author (model_author) <U14 'PyGEM-OGGM_v13' ... 'CISM2'\n",
" * year_after_2020 (year_after_2020) float32 0.0 1.0 ... 4.949e+03 4.95e+03\n",
" * rgi_reg (rgi_reg) <U2 '01' '02' '03' '04' ... '16' '17' '18' '19'\n",
" * experiments (experiments) object MultiIndex\n",
" * gcm (experiments) <U13 'gfdl-esm4' ... 'ukesm1-0-ll'\n",
" * period_scenario (experiments) <U16 '1851-1870_hist' ... '2081-2100_ssp585'\n",
"Data variables:\n",
" temp_ch_ipcc (experiments) float64 0.2314 0.4783 0.3923 ... 5.84 6.884\n",
" volume_m3 (model_author, rgi_reg, year_after_2020, experiments) float64 ...| \n", " | model_author | \n", "gcm | \n", "period_scenario | \n", "temp_ch_ipcc | \n", "volume_m3 | \n", "year_after_2020 | \n", "relative volume change (in %) | \n", "delta relative volume change (in %, to median model) | \n", "delta relative volume change (in %, to median model, relative to ref period ) | \n", "
|---|---|---|---|---|---|---|---|---|---|
| 0 | \n", "GLIMB | \n", "gfdl-esm4 | \n", "1851-1870_hist | \n", "0.231409 | \n", "1.253455e+14 | \n", "0.0 | \n", "80.812475 | \n", "-5.675031 | \n", "-5.984223 | \n", "
| 1 | \n", "GLIMB | \n", "gfdl-esm4 | \n", "1901-1920_hist | \n", "0.478289 | \n", "1.146157e+14 | \n", "0.0 | \n", "74.273069 | \n", "-1.113826 | \n", "-1.423018 | \n", "
| 2 | \n", "GLIMB | \n", "gfdl-esm4 | \n", "1951-1970_hist | \n", "0.392281 | \n", "1.179973e+14 | \n", "0.0 | \n", "76.486687 | \n", "-12.982556 | \n", "-13.291747 | \n", "
| 3 | \n", "GLIMB | \n", "gfdl-esm4 | \n", "1995-2014_hist | \n", "0.901467 | \n", "1.033774e+14 | \n", "0.0 | \n", "67.609008 | \n", "0.309192 | \n", "0.000000 | \n", "
| 4 | \n", "GLIMB | \n", "gfdl-esm4 | \n", "2021-2040_ssp126 | \n", "1.493792 | \n", "8.175907e+13 | \n", "0.0 | \n", "53.517849 | \n", "2.054803 | \n", "1.745611 | \n", "
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 315 | \n", "PyGEM-OGGM_v13 | \n", "ukesm1-0-ll | \n", "2061-2080_ssp370 | \n", "4.439977 | \n", "8.614016e+12 | \n", "0.0 | \n", "5.641338 | \n", "-3.235347 | \n", "11.543704 | \n", "
| 316 | \n", "PyGEM-OGGM_v13 | \n", "ukesm1-0-ll | \n", "2061-2080_ssp585 | \n", "5.230543 | \n", "5.856102e+12 | \n", "0.0 | \n", "3.827016 | \n", "-3.290553 | \n", "11.488499 | \n", "
| 317 | \n", "PyGEM-OGGM_v13 | \n", "ukesm1-0-ll | \n", "2081-2100_ssp126 | \n", "3.038482 | \n", "2.002791e+13 | \n", "0.0 | \n", "13.095330 | \n", "-9.744234 | \n", "5.034818 | \n", "
| 318 | \n", "PyGEM-OGGM_v13 | \n", "ukesm1-0-ll | \n", "2081-2100_ssp370 | \n", "5.840495 | \n", "5.058654e+12 | \n", "0.0 | \n", "3.311065 | \n", "-1.680703 | \n", "13.098349 | \n", "
| 319 | \n", "PyGEM-OGGM_v13 | \n", "ukesm1-0-ll | \n", "2081-2100_ssp585 | \n", "6.884361 | \n", "4.742040e+12 | \n", "0.0 | \n", "3.095073 | \n", "-1.056362 | \n", "13.722690 | \n", "
320 rows × 9 columns
\n", "<xarray.DataArray 'volume_m3' (model_author: 10, rgi_reg: 19, experiments: 80)>\n",
"array([[[ 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",
" [ 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",
" [[ 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",
" [5.7725011e+12, 5.2635689e+12, 5.3572304e+12, ...,\n",
" 3.9940713e+12, 2.2937477e+12, 1.8648834e+12],\n",
" [8.1662149e+10, 7.1292092e+10, 7.7415629e+10, ...,\n",
" 5.4494708e+10, 2.7098788e+10, 2.3415669e+10],\n",
" [ nan, nan, nan, ...,\n",
" nan, nan, nan]],\n",
"\n",
" [[ nan, nan, nan, ...,\n",
" nan, nan, nan],\n",
" [2.3770835e+11, 1.5806839e+11, 2.7891958e+11, ...,\n",
" 3.2304988e+09, 0.0000000e+00, 0.0000000e+00],\n",
" [ nan, nan, nan, ...,\n",
" nan, nan, nan],\n",
" ...,\n",
" [ nan, nan, nan, ...,\n",
" nan, nan, nan],\n",
" [8.7228482e+10, 6.7772969e+10, 8.2020876e+10, ...,\n",
" 3.2495129e+10, 4.9723812e+09, 2.5286986e+09],\n",
" [ nan, nan, nan, ...,\n",
" nan, nan, nan]]], dtype=float32)\n",
"Coordinates:\n",
" * model_author (model_author) <U14 'CISM2' 'GO' ... 'GloGEMflow3D'\n",
" * rgi_reg (rgi_reg) <U2 '01' '02' '03' '04' ... '16' '17' '18' '19'\n",
" extend_option <U18 ...\n",
" * experiments (experiments) object MultiIndex\n",
" * gcm (experiments) <U13 'gfdl-esm4' ... 'ukesm1-0-ll'\n",
" * period_scenario (experiments) <U16 '1851-1870_hist' ... '2081-2100_ssp585'\n",
" temp_ch_ipcc (experiments) float64 0.2314 0.4783 0.3923 ... 5.84 6.884