{ "cells": [ { "cell_type": "markdown", "id": "acb38dff-3de5-4be7-91cd-dfbeedd909ab", "metadata": {}, "source": [ "# Analyse the influence of using any GCM vs W5E5 on the specific MB during calibration time period:\n", "\n", "- idea: use the GCM ISIMIP3b data and check how different the specific MB from 2000-2019 is to the W5E5 calibrated specific MB" ] }, { "cell_type": "code", "execution_count": 127, "id": "e2f0a291-a7ff-488d-ad12-64734ad327b0", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-10-10 10:21:19: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.\n", "2022-10-10 10:21:19: oggm.cfg: Multiprocessing switched OFF according to the parameter file.\n", "2022-10-10 10:21:19: oggm.cfg: Multiprocessing: using all available processors (N=32)\n", "2022-10-10 10:21:19: oggm.cfg: Multiprocessing switched ON after user settings.\n", "2022-10-10 10:21:19: oggm.cfg: PARAMS['border'] changed from `40` to `10`.\n", "2022-10-10 10:21:19: oggm.cfg: Multiprocessing: using the requested number of processors (N=28)\n", "2022-10-10 10:21:20: oggm.workflow: Execute entity tasks [GlacierDirectory] on 3927 glaciers\n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import numpy as np\n", "import os\n", "\n", "import oggm\n", "from oggm import cfg, utils, workflow, tasks, graphics\n", "from oggm.core import massbalance, climate\n", "\n", "cfg.initialize(logging_level='WARNING')\n", "cfg.PARAMS['use_multiprocessing'] = True #True\n", "\n", "cfg.PATHS['working_dir'] = '/home/users/lschuster/oggm_w5e5_vs_gcm_gdir'\n", "\n", "cfg.PARAMS['border'] = 10\n", "cfg.PARAMS['mp_processes'] = 28\n", " \n", "base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/elev_bands/w5e5/qc0/pcpwin/match_geod_pergla'\n", "pd_geod = utils.get_geodetic_mb_dataframe()\n", "pd_geod = pd_geod.loc[pd_geod.period == '2000-01-01_2020-01-01']\n", "\n", "ssps = ['ssp126', 'ssp585']\n", "gcms = ['ukesm1-0-ll_r1i1p1f2', 'gfdl-esm4_r1i1p1f1',\n", " 'ipsl-cm6a-lr_r1i1p1f1', 'mpi-esm1-2-hr_r1i1p1f1',\n", " 'mri-esm2-0_r1i1p1f1']\n", "climate_type = 'W5E5'\n", "y0 = 2000\n", "y1 = 2019\n", "reload = False\n", "if reload:\n", " gdirs = workflow.init_glacier_directories(pd_geod[pd_geod.reg == 11].index, from_prepro_level=3, prepro_base_url=base_url)\n", " for ensemble in gcms:\n", " for ssp in ssps:\n", " for correction in [True, False]:\n", " if correction is True:\n", " correct=''\n", " else:\n", " correct='_no_correction'\n", " workflow.execute_entity_task(oggm.shop.gcm_climate.process_monthly_isimip_data, gdirs, ensemble = ensemble,\n", " ssp = ssp, output_filesuffix=f'{ensemble}_{ssp}{correct}',\n", " apply_bias_correction = correction)\n", " \n", " # also do the bias correction on the calibration time period only --> does this help to reduce the differences???\n", " for ensemble in gcms:\n", " for ssp in ssps:\n", " workflow.execute_entity_task(oggm.shop.gcm_climate.process_monthly_isimip_data, gdirs, year_range=(str(y0), str(y1)),\n", " ssp = ssp, ensemble = ensemble,\n", " output_filesuffix=f'{ensemble}_{ssp}_bc_{y0}_{y1}',\n", " apply_bias_correction=True)\n", "else:\n", " gdirs = workflow.init_glacier_directories(pd_geod[pd_geod.reg == 11].index) #, from_prepro_level=3, prepro_base_url=base_url)\n", "\n", " " ] }, { "cell_type": "code", "execution_count": 169, "id": "e9774075-1463-4abd-b724-82b69ad24187", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 169, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correct" ] }, { "cell_type": "code", "execution_count": 176, "id": "1cf29e1e-6e7c-4353-a906-c51701f70bf6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "t\n" ] } ], "source": [ " if correction is True:\n", " correct=''\n", " elif correction == f'True: yr {y0}--{y1}':\n", " print('t')\n", " correct = f'_bc_{y0}_{y1}'\n", " else:\n", " correct='_no_correction'" ] }, { "cell_type": "code", "execution_count": 175, "id": "55831226-949d-4d2d-be1f-6c493785c3b6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'True: yr 2000--2019'" ] }, "execution_count": 175, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correction" ] }, { "cell_type": "code", "execution_count": null, "id": "77525f8c-bcb5-4bb3-bfb2-790a08f579df", "metadata": {}, "outputs": [], "source": [ " hgts, widths = gdir.get_inversion_flowline_hw()\n", " mb_mod = massbalance.PastMassBalance(gdir)\n", " spec_mb_2000_2019_w5e5 = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", " spec_mb_2000_2019_obs = pd_geod.loc[gdir.rgi_id].dmdtda*1000\n", " for ensemble in gcms:\n", " for ssp in ssps:\n", " for correction in [f'True: yr {y0}--{y1}']: #[True, False]\n", " if correction is True:\n", " correct=''\n", " elif correction == f'True: yr {y0}--{y1}':\n", " correct = f'_bc_{y0}_{y1}'\n", " else:\n", " correct='_no_correction'\n", " mb_mod = massbalance.PastMassBalance(gdir, filename='gcm_data',\n", " input_filesuffix = f'{ensemble}_{ssp}{correct}')\n", " spec_mb_2000_2019_gcm = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", "\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'rgi_id'] = gdir.rgi_id\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'gcm'] = ensemble\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'ssp'] = ssp\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'apply_bias_correction'] = correction\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_gcm'] = spec_mb_2000_2019_gcm\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_w5e5'] = spec_mb_2000_2019_w5e5\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_obs'] = spec_mb_2000_2019_obs\n", " j += 1" ] }, { "cell_type": "code", "execution_count": 177, "id": "4ea87699-d4fb-43a0-bcde-bbe54caddd16", "metadata": {}, "outputs": [], "source": [ "y0 = 2000\n", "y1 = 2019\n", "load = False\n", "if load:\n", " pd_spec_mb_gcm_vs_w5e5 = pd.DataFrame(np.NaN, index=np.arange(0, 1000000), #*len(gdirs))),\n", " columns=['rgi_id', 'gcm', 'ssp',\n", " 'spec_mb_2000_2019_gcm', 'spec_mb_2000_2019_w5e5',\n", " 'spec_mb_2000_2019_obs', 'apply_bias_correction'])\n", "\n", " j = 0\n", " for gdir in gdirs:\n", " try:\n", " hgts, widths = gdir.get_inversion_flowline_hw()\n", " mb_mod = massbalance.PastMassBalance(gdir)\n", " spec_mb_2000_2019_w5e5 = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", " spec_mb_2000_2019_obs = pd_geod.loc[gdir.rgi_id].dmdtda*1000\n", " for ensemble in gcms:\n", " for ssp in ssps:\n", " for correction in [f'True: yr {y0}--{y1}']: #[True, False]\n", " if correction is True:\n", " correct=''\n", " elif correction == f'True: yr {y0}--{y1}':\n", " correct = f'_bc_{y0}_{y1}'\n", " else:\n", " correct='_no_correction'\n", " mb_mod = massbalance.PastMassBalance(gdir, filename='gcm_data',\n", " input_filesuffix = f'{ensemble}_{ssp}{correct}')\n", " spec_mb_2000_2019_gcm = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", "\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'rgi_id'] = gdir.rgi_id\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'gcm'] = ensemble\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'ssp'] = ssp\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'apply_bias_correction'] = correction\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_gcm'] = spec_mb_2000_2019_gcm\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_w5e5'] = spec_mb_2000_2019_w5e5\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_obs'] = spec_mb_2000_2019_obs\n", " j += 1\n", "\n", " except:\n", " pass\n", " pd_spec_mb_gcm_vs_w5e5 = pd_spec_mb_gcm_vs_w5e5.dropna()\n", " pd_spec_mb_gcm_vs_w5e5['difference_gcm_vs_w5e5'] = pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_gcm'] - pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", " pd_spec_mb_gcm_vs_w5e5['rel_perc_difference_gcm_vs_w5e5'] = 100 * pd_spec_mb_gcm_vs_w5e5['difference_gcm_vs_w5e5'] / pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", " pd_spec_mb_gcm_vs_w5e5['ratio_gcm_vs_w5e5'] = pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_gcm'] / pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", " pd_spec_mb_gcm_vs_w5e5_bc_2000_2019 = pd_spec_mb_gcm_vs_w5e5.copy()\n", " pd_spec_mb_gcm_vs_w5e5_bc_2000_2019.to_csv(f'alpine_glacier_spec_mb_gcm_vs_w5e5_bc_{y0}_{y1}.csv')\n", "else:\n", " pd_spec_mb_gcm_vs_w5e5_bc_2000_2019 = pd.read_csv(f'alpine_glacier_spec_mb_gcm_vs_w5e5_bc_{y0}_{y1}.csv', index_col=0)\n" ] }, { "cell_type": "code", "execution_count": 73, "id": "c5fdd793-a608-4d62-bb74-ba9d413333dd", "metadata": {}, "outputs": [], "source": [ "load = False\n", "if load:\n", " pd_spec_mb_gcm_vs_w5e5 = pd.DataFrame(np.NaN, index=np.arange(0, 1000000), #*len(gdirs))),\n", " columns=['rgi_id', 'gcm', 'ssp',\n", " 'spec_mb_2000_2019_gcm', 'spec_mb_2000_2019_w5e5',\n", " 'spec_mb_2000_2019_obs', 'apply_bias_correction'])\n", "\n", " j = 0\n", " for gdir in gdirs:\n", " try:\n", " hgts, widths = gdir.get_inversion_flowline_hw()\n", " mb_mod = massbalance.PastMassBalance(gdir)\n", " spec_mb_2000_2019_w5e5 = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", " spec_mb_2000_2019_obs = pd_geod.loc[gdir.rgi_id].dmdtda*1000\n", " for ensemble in gcms:\n", " for ssp in ssps:\n", " for correction in [True, False]:\n", " if correction:\n", " correct=''\n", " else:\n", " correct='_no_correction'\n", " mb_mod = massbalance.PastMassBalance(gdir, filename='gcm_data', input_filesuffix = f'{ensemble}_{ssp}{correct}')\n", " spec_mb_2000_2019_gcm = mb_mod.get_specific_mb(heights=hgts, widths=widths,\n", " year=np.arange(2000, # yrs_seasonal_mbs[0]\n", " 2019 + 1, 1)).mean()\n", "\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'rgi_id'] = gdir.rgi_id\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'gcm'] = ensemble\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'ssp'] = ssp\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'apply_bias_correction'] = correction\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_gcm'] = spec_mb_2000_2019_gcm\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_w5e5'] = spec_mb_2000_2019_w5e5\n", " pd_spec_mb_gcm_vs_w5e5.loc[j, 'spec_mb_2000_2019_obs'] = spec_mb_2000_2019_obs\n", " j += 1\n", "\n", " except:\n", " pass\n", " pd_spec_mb_gcm_vs_w5e5 = pd_spec_mb_gcm_vs_w5e5.dropna()\n", " pd_spec_mb_gcm_vs_w5e5['difference_gcm_vs_w5e5'] = pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_gcm'] - pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", " pd_spec_mb_gcm_vs_w5e5['rel_perc_difference_gcm_vs_w5e5'] = 100 * pd_spec_mb_gcm_vs_w5e5['difference_gcm_vs_w5e5'] / pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", " pd_spec_mb_gcm_vs_w5e5['ratio_gcm_vs_w5e5'] = pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_gcm'] / pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs']\n", "\n", " pd_spec_mb_gcm_vs_w5e5.to_csv('alpine_glacier_spec_mb_gcm_vs_w5e5.csv')\n", "else:\n", " pd_spec_mb_gcm_vs_w5e5 = pd.read_csv('alpine_glacier_spec_mb_gcm_vs_w5e5.csv', index_col=0)\n" ] }, { "cell_type": "code", "execution_count": 179, "id": "90598449-38fb-4cd1-9115-0f48603f83eb", "metadata": {}, "outputs": [], "source": [ "_pd_spec_mb_gcm_vs_w5e5 = pd.read_csv('alpine_glacier_spec_mb_gcm_vs_w5e5.csv', index_col=0)\n", "_pd_spec_mb_gcm_vs_w5e5_bc_2000_2019 = pd.read_csv(f'alpine_glacier_spec_mb_gcm_vs_w5e5_bc_{y0}_{y1}.csv', index_col=0)\n", "pd_spec_mb_gcm_vs_w5e5 = pd.concat([_pd_spec_mb_gcm_vs_w5e5, _pd_spec_mb_gcm_vs_w5e5_bc_2000_2019])\n" ] }, { "cell_type": "code", "execution_count": 180, "id": "e333a649-3a0c-4765-b316-5d38b442e91b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['RGI60-11.03778' 'RGI60-11.03818']\n" ] } ], "source": [ "# remove those glaciers with \"0\" observed specific MB out of the statistic\n", "rgi_w_0_mb = pd_spec_mb_gcm_vs_w5e5.iloc[np.where(pd_spec_mb_gcm_vs_w5e5['spec_mb_2000_2019_obs'] == 0)].rgi_id.unique()\n", "print(rgi_w_0_mb)\n", "for r in rgi_w_0_mb:\n", " pd_spec_mb_gcm_vs_w5e5 = pd_spec_mb_gcm_vs_w5e5.iloc[np.where(pd_spec_mb_gcm_vs_w5e5.index!=r)]" ] }, { "cell_type": "code", "execution_count": 181, "id": "d58d98fb-c046-403d-a610-ae858b9d60b8", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1875130/1366677533.py:4: FutureWarning: Dropping invalid columns in DataFrameGroupBy.quantile is deprecated. In a future version, a TypeError will be raised. Before calling .quantile, select only columns which should be valid for the function.\n", " pd_spec_mb_gcm_vs_w5e5.groupby(by='apply_bias_correction').quantile([0.05, 0.25,0.50,0.75, 0.95])\n" ] }, { "data": { "text/html": [ "
| \n", " | \n", " | spec_mb_2000_2019_gcm | \n", "spec_mb_2000_2019_w5e5 | \n", "spec_mb_2000_2019_obs | \n", "difference_gcm_vs_w5e5 | \n", "rel_perc_difference_gcm_vs_w5e5 | \n", "ratio_gcm_vs_w5e5 | \n", "
|---|---|---|---|---|---|---|---|
| apply_bias_correction | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
| False | \n", "0.05 | \n", "-1395.838741 | \n", "-1.156900e+03 | \n", "-1156.900000 | \n", "-501.315651 | \n", "-87.099273 | \n", "0.129007 | \n", "
| 0.25 | \n", "-893.899883 | \n", "-7.943000e+02 | \n", "-794.300000 | \n", "-164.533389 | \n", "-21.227245 | \n", "0.787728 | \n", "|
| 0.50 | \n", "-630.625525 | \n", "-5.963719e+02 | \n", "-596.371893 | \n", "-15.245944 | \n", "2.332115 | \n", "1.023321 | \n", "|
| 0.75 | \n", "-371.430970 | \n", "-4.184000e+02 | \n", "-418.400000 | \n", "122.811305 | \n", "29.414121 | \n", "1.294141 | \n", "|
| 0.95 | \n", "68.291317 | \n", "-1.264766e-13 | \n", "0.000000 | \n", "320.563072 | \n", "124.210839 | \n", "2.242108 | \n", "|
| True | \n", "0.05 | \n", "-1351.567789 | \n", "-1.156900e+03 | \n", "-1156.900000 | \n", "-494.344949 | \n", "-94.705103 | \n", "0.052949 | \n", "
| 0.25 | \n", "-852.582931 | \n", "-7.943000e+02 | \n", "-794.300000 | \n", "-129.348306 | \n", "-26.050633 | \n", "0.739494 | \n", "|
| 0.50 | \n", "-597.043709 | \n", "-5.963719e+02 | \n", "-596.371893 | \n", "17.549061 | \n", "-2.379314 | \n", "0.976207 | \n", "|
| 0.75 | \n", "-343.222261 | \n", "-4.184000e+02 | \n", "-418.400000 | \n", "152.655583 | \n", "22.497070 | \n", "1.224971 | \n", "|
| 0.95 | \n", "97.625642 | \n", "-1.264766e-13 | \n", "0.000000 | \n", "341.586502 | \n", "113.835016 | \n", "2.138350 | \n", "|
| True: yr 2000--2019 | \n", "0.05 | \n", "-1166.958733 | \n", "-1.156900e+03 | \n", "-1156.900000 | \n", "-73.725384 | \n", "-12.439902 | \n", "0.875601 | \n", "
| 0.25 | \n", "-811.165777 | \n", "-7.943000e+02 | \n", "-794.300000 | \n", "-37.753931 | \n", "-1.780337 | \n", "0.982197 | \n", "|
| 0.50 | \n", "-617.705575 | \n", "-5.963719e+02 | \n", "-596.371893 | \n", "-15.063724 | \n", "2.169556 | \n", "1.021696 | \n", "|
| 0.75 | \n", "-431.358401 | \n", "-4.184000e+02 | \n", "-418.400000 | \n", "8.896489 | \n", "6.350496 | \n", "1.063505 | \n", "|
| 0.95 | \n", "-8.802774 | \n", "-1.264766e-13 | \n", "0.000000 | \n", "44.701444 | \n", "20.759904 | \n", "1.207599 | \n", "