{ "cells": [ { "cell_type": "markdown", "id": "bb8851ab-cd1c-495a-bfee-3411d36bfa83", "metadata": { "tags": [] }, "source": [ "# VAS data\n", "\n", "This is just an example workflow that shows how to :\n", "1. [Analyse amount of errors and relative failing glacier area when no glacier_statistics.csv files were calculated beforehand](#For-simplifications,-let's-just-look-at-the-errors-from-a-single-GCM-and-SSP!)\n", "2. Compare the future common non-failing glacier volumes of different calibration methods (globally and regionally):\n", " - a) [for one gcm and one ssp](#Volume-projections-for-common-non-failing-glaciers-for-one-gcm-and-one-ssp)\n", " - b) [for all gcms and one ssp by looking at the mean and standard deviation](#Volume-projections-for-common-non-failing-glaciers-for-all-gcms-that-work-for-the-lowest-emission-ssp-scenario-ssp126)" ] }, { "cell_type": "markdown", "id": "3613e525-4472-439e-a008-1e442d762f0f", "metadata": {}, "source": [ "**Here we look at these three different calibration options:**\n", "\n", "cmip6_old_output:\n", " - CRU with regional geodetic MB calib and QC3\n", " - match_geod_qc3\n", " \n", "cmip6_output: \n", " - CRU with glacier-specific geodetic MB calib and QC0 \n", " - match_geod_per_glac_qc0\n", " \n", "cmip6_output_massredis:\n", " - CRU with glacier specific geodetic MB calib and mass redistribution instead of SIA and QC0\n", " - match_geod_per_glac_massredis_qc0" ] }, { "cell_type": "markdown", "id": "c101a3c5-785e-4511-95fc-079f48d25a92", "metadata": {}, "source": [ "---" ] }, { "cell_type": "code", "execution_count": 1, "id": "bdc2c077-0f2b-4eed-a9dd-64131a2fc5f0", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/users/lschuster/.local/lib/python3.9/site-packages/xarray/backends/cfgrib_.py:27: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message\n", " warnings.warn(\n" ] } ], "source": [ "from oggm import cfg, workflow, utils, shop\n", "import pandas as pd\n", "import os, glob\n", "import numpy as np\n", "import xarray as xr\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "sns.set_style(\"whitegrid\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "4770befc-5970-4e97-821e-27c20201cfb4", "metadata": {}, "outputs": [], "source": [ "gcms = os.listdir('/home/www/fmaussion/vas_paper/runs_cmip6_cru/cmip6_output_massredis/RGI01')\n", "ssps = ['ssp126','ssp245', 'ssp370', 'ssp585'] \n", "\n" ] }, { "cell_type": "markdown", "id": "d7cc2287-6d75-465a-b301-6f3621b5bb50", "metadata": {}, "source": [ "### For simplifications, let's just look at the errors from a single GCM and SSP!" ] }, { "cell_type": "code", "execution_count": 4, "id": "d969c154-b4f8-4d9c-9a36-4692e4c95219", "metadata": {}, "outputs": [], "source": [ "# as an example:\n", "#gcm = 'BCC-CSM2-MR'\n", "gcm = 'NorESM2-MM'\n", "#gcm = 'CESM2'\n", "ssp = 'ssp126'" ] }, { "cell_type": "code", "execution_count": 5, "id": "35461fff-33fc-4de4-83f2-d90dfbdd4d2a", "metadata": {}, "outputs": [], "source": [ "def error_stats_working_rgi(gcm = 'BCC-CSM2-MR', ssp = 'ssp585'):\n", " \n", " convert_dict = {'cmip6_output_massredis': 'match_geod_per_glac_massredis_qc0',\n", " 'cmip6_output':'match_geod_per_glac_qc0',\n", " 'cmip6_old_output': 'match_geod_qc3'}\n", " pd_geodetic = utils.get_geodetic_mb_dataframe()[utils.get_geodetic_mb_dataframe().period=='2000-01-01_2020-01-01']\n", " # we do not model RGI region 19 here ... \n", " pd_geodetic = pd_geodetic[pd_geodetic.reg != 19]\n", "\n", " \n", " total_area = pd_geodetic.area.sum()\n", " total_counts = len(pd_geodetic)\n", " # pd_geodetic.area\n", " pd_error = pd.DataFrame(index = pd_geodetic.index,\n", " columns=['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output', 'area', 'all_running'])\n", " pd_error['area'] = pd_geodetic.area\n", " pd_error['all_running'] = np.NaN\n", "\n", "\n", " for cmip6_output_type in ['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']:\n", " working_rgis = []\n", " # get working glaciers for each rgi and get a global list of working glaciers\n", " for n in np.arange(1,19,1):\n", " if n<10:\n", " rgi_reg = f'RGI0{n}'\n", " else:\n", " rgi_reg = f'RGI{n}'\n", " ds = xr.open_dataset(f'/home/www/fmaussion/vas_paper/runs_cmip6_cru/{cmip6_output_type}/{rgi_reg}/{gcm}/{gcm}_{ssp}.nc').volume.dropna(dim='rgi_id')\n", " working_rgi = ds.isel(time=0).dropna(dim='rgi_id').rgi_id.values\n", " working_rgis.append(working_rgi)\n", " ds.close()\n", " working_rgis = np.concatenate(working_rgis)\n", "\n", " # pd_error[cmip6_output_type] = False\n", " # set those glaciers that worked to true, the other are NaN\n", " pd_error.loc[working_rgis,cmip6_output_type] = True\n", " # check if it worked\n", " assert len(working_rgis) == len(pd_error.loc[pd_error[cmip6_output_type] == True])\n", " #print(len(working_rgis))\n", " # this is slower\n", " #for rgi in working_rgis:\n", " # if rgi in pd_test.index:\n", " # pd_test.loc[rgi][cmip6_output_type] = True\n", "\n", " # get those glaciers that work for all experiment variants ('but only for that ssp and gcm'!!!)\n", " all_running_rgis = pd_error[['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']].dropna()\n", " pd_error.loc[all_running_rgis.index, 'all_running'] = True\n", " # len(pd_error['all_running'].dropna()) <= \n", "\n", " # get relative glacier area that works\n", " pd_error_stats = {}\n", " for cmip6_output_type in ['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']:\n", " pd_error_stats[f'{convert_dict[cmip6_output_type]}_area'] = (1-(pd_error.loc[pd_error[cmip6_output_type].dropna().index]['area'].sum())/total_area)*100\n", " pd_error_stats['all_running_area'] = (1-(pd_error.loc[pd_error['all_running'].dropna().index]['area'].sum())/total_area)*100\n", " for cmip6_output_type in ['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']:\n", " pd_error_stats[f'{convert_dict[cmip6_output_type]}_counts'] = (1-len(pd_error[cmip6_output_type].dropna().index)/total_counts)*100\n", " pd_error_stats['all_running_counts'] = (1-len(pd_error['all_running'].dropna().index)/total_counts)*100\n", " \n", " pd_error_rgi_stats = {}\n", " for reg in set(pd_geodetic.reg):\n", " pd_error_rgi_stats_reg = {}\n", " rgi_reg = pd_geodetic.loc[pd_geodetic.reg==reg].index\n", " total_reg_area = pd_geodetic.loc[pd_geodetic.reg==reg].area.sum()\n", " total_reg_counts = len(pd_geodetic.loc[pd_geodetic.reg==reg])\n", "\n", " pd_error.loc[rgi_reg, 'reg'] = reg\n", " for cmip6_output_type in ['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']:\n", " pd_error_rgi_stats_reg[f'{convert_dict[cmip6_output_type]}_area'] = (1-(pd_error.loc[pd_error.loc[rgi_reg][cmip6_output_type].dropna().index]['area'].sum())/total_reg_area)*100\n", " pd_error_rgi_stats_reg['all_running_area'] = (1-(pd_error.loc[pd_error.loc[rgi_reg]['all_running'].dropna().index]['area'].sum())/total_reg_area)*100\n", " for cmip6_output_type in ['cmip6_output_massredis', 'cmip6_output', 'cmip6_old_output']:\n", " pd_error_rgi_stats_reg[f'{convert_dict[cmip6_output_type]}_counts'] = (1-len(pd_error.loc[pd_error.loc[rgi_reg][cmip6_output_type].dropna().index])/total_reg_counts)*100\n", " pd_error_rgi_stats_reg['all_running_counts'] = (1-len(pd_error.loc[pd_error.loc[rgi_reg]['all_running'].dropna().index])/total_reg_counts)*100\n", " pd_error_rgi_stats[str(reg)] = pd_error_rgi_stats_reg\n", "\n", " #return pd_error_stats, pd_error, pd_error_rgi_stats\n", " pd_error_rgi_stats = pd.DataFrame(pd_error_rgi_stats)\n", " pd_error_rgi_stats['all'] = pd.DataFrame(pd_error_stats, index=[0]).T[0] #pd_error_stats\n", " return pd_error_stats, pd_error, pd_error_rgi_stats" ] }, { "cell_type": "code", "execution_count": 6, "id": "d74ab2b9-2825-4b5a-a9e7-1832423b801d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/users/lschuster/.local/lib/python3.9/site-packages/xarray/backends/plugins.py:61: RuntimeWarning: Engine 'cfgrib' loading failed:\n", "Cannot find the ecCodes library\n", " warnings.warn(f\"Engine {name!r} loading failed:\\n{ex}\", RuntimeWarning)\n" ] }, { "data": { "text/html": [ "
<xarray.DataArray 'volume' (time: 82, rgi_id: 27080)>\n", "array([[5.755895e+06, 1.527717e+07, 4.210734e+07, ..., 1.067689e+08,\n", " 3.436057e+07, 1.805127e+06],\n", " [5.627826e+06, 1.499158e+07, 4.002556e+07, ..., 1.009200e+08,\n", " 3.218732e+07, 1.687626e+06],\n", " [5.598141e+06, 1.484801e+07, 3.837280e+07, ..., 9.553950e+07,\n", " 3.032915e+07, 1.597355e+06],\n", " ...,\n", " [1.763794e+06, 3.596061e+06, 1.764006e+05, ..., 1.712339e+07,\n", " 5.545750e+06, 0.000000e+00],\n", " [1.722843e+06, 3.471366e+06, 1.568110e+05, ..., 1.672752e+07,\n", " 5.437562e+06, 0.000000e+00],\n", " [1.691012e+06, 3.376146e+06, 1.378086e+05, ..., 1.707277e+07,\n", " 5.520620e+06, 0.000000e+00]], dtype=float32)\n", "Coordinates:\n", " * time (time) float64 2.02e+03 2.021e+03 ... 2.1e+03 2.101e+03\n", " * rgi_id (rgi_id) object 'RGI60-01.00001' ... 'RGI60-01.27112'\n", " hydro_year (time) int64 2020 2021 2022 2023 ... 2098 2099 2100 2101\n", " hydro_month (time) int64 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1\n", " calendar_year (time) int64 2020 2021 2022 2023 ... 2098 2099 2100 2101\n", " calendar_month (time) int64 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1\n", "Attributes:\n", " description: Total glacier volume\n", " unit: m 3
array([[5.755895e+06, 1.527717e+07, 4.210734e+07, ..., 1.067689e+08,\n", " 3.436057e+07, 1.805127e+06],\n", " [5.627826e+06, 1.499158e+07, 4.002556e+07, ..., 1.009200e+08,\n", " 3.218732e+07, 1.687626e+06],\n", " [5.598141e+06, 1.484801e+07, 3.837280e+07, ..., 9.553950e+07,\n", " 3.032915e+07, 1.597355e+06],\n", " ...,\n", " [1.763794e+06, 3.596061e+06, 1.764006e+05, ..., 1.712339e+07,\n", " 5.545750e+06, 0.000000e+00],\n", " [1.722843e+06, 3.471366e+06, 1.568110e+05, ..., 1.672752e+07,\n", " 5.437562e+06, 0.000000e+00],\n", " [1.691012e+06, 3.376146e+06, 1.378086e+05, ..., 1.707277e+07,\n", " 5.520620e+06, 0.000000e+00]], dtype=float32)
array([2020., 2021., 2022., 2023., 2024., 2025., 2026., 2027., 2028., 2029.,\n", " 2030., 2031., 2032., 2033., 2034., 2035., 2036., 2037., 2038., 2039.,\n", " 2040., 2041., 2042., 2043., 2044., 2045., 2046., 2047., 2048., 2049.,\n", " 2050., 2051., 2052., 2053., 2054., 2055., 2056., 2057., 2058., 2059.,\n", " 2060., 2061., 2062., 2063., 2064., 2065., 2066., 2067., 2068., 2069.,\n", " 2070., 2071., 2072., 2073., 2074., 2075., 2076., 2077., 2078., 2079.,\n", " 2080., 2081., 2082., 2083., 2084., 2085., 2086., 2087., 2088., 2089.,\n", " 2090., 2091., 2092., 2093., 2094., 2095., 2096., 2097., 2098., 2099.,\n", " 2100., 2101.])
array(['RGI60-01.00001', 'RGI60-01.00002', 'RGI60-01.00003', ...,\n", " 'RGI60-01.27110', 'RGI60-01.27111', 'RGI60-01.27112'], dtype=object)
array([2020, 2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030, 2031,\n", " 2032, 2033, 2034, 2035, 2036, 2037, 2038, 2039, 2040, 2041, 2042, 2043,\n", " 2044, 2045, 2046, 2047, 2048, 2049, 2050, 2051, 2052, 2053, 2054, 2055,\n", " 2056, 2057, 2058, 2059, 2060, 2061, 2062, 2063, 2064, 2065, 2066, 2067,\n", " 2068, 2069, 2070, 2071, 2072, 2073, 2074, 2075, 2076, 2077, 2078, 2079,\n", " 2080, 2081, 2082, 2083, 2084, 2085, 2086, 2087, 2088, 2089, 2090, 2091,\n", " 2092, 2093, 2094, 2095, 2096, 2097, 2098, 2099, 2100, 2101])
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
array([2020, 2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030, 2031,\n", " 2032, 2033, 2034, 2035, 2036, 2037, 2038, 2039, 2040, 2041, 2042, 2043,\n", " 2044, 2045, 2046, 2047, 2048, 2049, 2050, 2051, 2052, 2053, 2054, 2055,\n", " 2056, 2057, 2058, 2059, 2060, 2061, 2062, 2063, 2064, 2065, 2066, 2067,\n", " 2068, 2069, 2070, 2071, 2072, 2073, 2074, 2075, 2076, 2077, 2078, 2079,\n", " 2080, 2081, 2082, 2083, 2084, 2085, 2086, 2087, 2088, 2089, 2090, 2091,\n", " 2092, 2093, 2094, 2095, 2096, 2097, 2098, 2099, 2100, 2101])
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
\n", " | 1 | \n", "2 | \n", "3 | \n", "4 | \n", "5 | \n", "6 | \n", "7 | \n", "8 | \n", "9 | \n", "10 | \n", "11 | \n", "12 | \n", "13 | \n", "14 | \n", "15 | \n", "16 | \n", "17 | \n", "18 | \n", "all | \n", "
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
match_geod_per_glac_massredis_qc0_area | \n", "0.004060 | \n", "0.012689 | \n", "0.020457 | \n", "0.031843 | \n", "0.388694 | \n", "0.002450 | \n", "0.445790 | \n", "0.003391 | \n", "1.784461 | \n", "3.576024 | \n", "0.017876 | \n", "9.125319 | \n", "6.134486 | \n", "0.178508 | \n", "0.167026 | \n", "0.000641 | \n", "4.780309 | \n", "0.009812 | \n", "1.081192 | \n", "
match_geod_per_glac_qc0_area | \n", "0.002438 | \n", "0.007773 | \n", "0.014666 | \n", "0.029387 | \n", "0.247575 | \n", "9.655949 | \n", "1.728935 | \n", "0.003391 | \n", "0.063121 | \n", "3.597891 | \n", "0.007934 | \n", "9.124293 | \n", "5.419207 | \n", "0.072792 | \n", "0.099130 | \n", "0.000641 | \n", "0.270851 | \n", "0.000000 | \n", "0.864427 | \n", "
match_geod_qc3_area | \n", "0.004393 | \n", "0.029007 | \n", "0.143415 | \n", "0.000942 | \n", "12.742051 | \n", "0.000000 | \n", "0.000000 | \n", "0.003391 | \n", "0.000000 | \n", "3.513287 | \n", "0.007934 | \n", "9.236786 | \n", "1.744749 | \n", "0.037249 | \n", "0.758887 | \n", "7.549307 | \n", "5.023888 | \n", "0.000000 | \n", "2.520817 | \n", "
all_running_area | \n", "0.008453 | \n", "0.041999 | \n", "0.161030 | \n", "0.034607 | \n", "12.984928 | \n", "9.656826 | \n", "2.174725 | \n", "0.003391 | \n", "1.800247 | \n", "3.615193 | \n", "0.017876 | \n", "9.237812 | \n", "6.850611 | \n", "0.217342 | \n", "0.924373 | \n", "7.549307 | \n", "9.802862 | \n", "0.009812 | \n", "3.743180 | \n", "