{ "cells": [ { "cell_type": "markdown", "id": "53611111-53e4-400d-a4cf-2e164578622b", "metadata": {}, "source": [ "# Combine monthly and annual runoff and other information into merged per-glacier data-files for every basin\n", "\n", "- The annual or monthly glacier runoff that is computed here is the sum of annual melt and liquid precipitation on and off the glacier using a fixed-gauge with a glacier minimum reference area from year 2000.\n", "\n", "- Attention, no filling is done, those glaciers that did not run have np.NaN values\n", "- two options:\n", " - `gcm_from_2000` : gcm directly from 2000 (results in differences between GCMs from 2000-2019) \n", " - `w5e5_gcm_merged` : W5E5 from 2000-2019 and GCM from 2020 runoff\n", "\n", "(the data got updated in March 2025, because some basins were accidentally missing)" ] }, { "cell_type": "code", "execution_count": 1, "id": "330447d8-5af9-4e1f-878d-31dc6565aa30", "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np\n", "import pandas as pd\n", "from time import gmtime, strftime\n", "import oggm.cfg\n", "from oggm import utils, workflow, tasks, graphics\n", "import json\n", "import geopandas as gpd\n", "\n", "# let's also do some visualisations ...\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "gcms_cmip6 = pd.read_csv('/home/www/oggm/cmip6/all_gcm_list.csv', index_col=0) \n", "all_GCM = gcms_cmip6.gcm.unique()\n", "all_scenario = gcms_cmip6.ssp.unique()\n", "\n", "# get the dataset where coordinates of glaciers are stored\n", "frgi = utils.file_downloader('https://cluster.klima.uni-bremen.de/~oggm/rgi/rgi62_stats.h5')\n", "#frgi = '/home/users/lschuster/glacierMIP/rgi62_stats.h5'\n", "odf = pd.read_hdf(frgi, index_col=0)\n", "\n", "# get the RGI ids of the glaciers of the chosen basin\n", "pd_basin_num = gpd.read_file('/home/www/fmaussion/misc/magicc/basins_shape/glacier_basins.shp')" ] }, { "cell_type": "code", "execution_count": 2, "id": "2003ba28-35f4-498e-af6b-3457962cfd51", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "73.0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### all available gcm_scenario combinations ... \n", "len(gcms_cmip6)/2" ] }, { "cell_type": "code", "execution_count": 3, "id": "5209d22b-bcd2-4854-b4b0-e52f30a08b97", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MRBIDRIVER_BASICONTINENTOCEANSEAAREA_CALCShape_LengShape_AreaRGI_AREAgeometry
02103INDIGIRKAAsiaArctic OceanEast Siberian Sea341234.085.60405170.212499171.941POLYGON ((151.725 70.975, 151.725 70.97083, 15...
12108OBAsiaArctic OceanKara Sea3040606.1168.355250448.342075763.493POLYGON ((91.75 57.70417, 91.74542 57.70299, 9...
22302BRAHMAPUTRAAsiaIndian OceanBay of Bengal540782.566.05444049.68661110528.496POLYGON ((97.76667 28.77083, 97.76335 28.77168...
32306GANGESAsiaIndian OceanBay of Bengal1006558.6122.34998390.9468807906.081MULTIPOLYGON (((88.02555 21.5715, 88.02361 21....
42309INDUSAsiaIndian OceanArabian Sea865012.681.28629182.85227527206.546MULTIPOLYGON (((67.44167 23.97006, 67.44028 23...
.................................
706241POEuropeAtlantic OceanAdriatic Sea73289.521.5760988.425867313.417POLYGON ((9.68333 46.425, 9.68418 46.42168, 9....
716242RHINEEuropeAtlantic OceanNorth Sea163122.348.24761220.200160336.922POLYGON ((4.45638 52.14527, 4.46175 52.1411, 4...
726243RHONEEuropeAtlantic OceanMediterranean Sea96637.628.00491611.197342917.302POLYGON ((6.84167 47.82083, 6.84167 47.81667, ...
736254THJORSAEuropeAtlantic OceanNorth Atlantic8277.89.2708471.545336972.419POLYGON ((-17.5 64.625, -17.50035 64.60591, -1...
746255TORNEALVEN (also TORNIONJOKI, also TORNIONVAYLA)EuropeAtlantic OceanBaltic Sea40545.228.9190258.68565234.257POLYGON ((23.69167 68.84167, 23.70214 68.84102...
\n", "

75 rows × 10 columns

\n", "
" ], "text/plain": [ " MRBID RIVER_BASI CONTINENT \\\n", "0 2103 INDIGIRKA Asia \n", "1 2108 OB Asia \n", "2 2302 BRAHMAPUTRA Asia \n", "3 2306 GANGES Asia \n", "4 2309 INDUS Asia \n", ".. ... ... ... \n", "70 6241 PO Europe \n", "71 6242 RHINE Europe \n", "72 6243 RHONE Europe \n", "73 6254 THJORSA Europe \n", "74 6255 TORNEALVEN (also TORNIONJOKI, also TORNIONVAYLA) Europe \n", "\n", " OCEAN SEA AREA_CALC Shape_Leng Shape_Area \\\n", "0 Arctic Ocean East Siberian Sea 341234.0 85.604051 70.212499 \n", "1 Arctic Ocean Kara Sea 3040606.1 168.355250 448.342075 \n", "2 Indian Ocean Bay of Bengal 540782.5 66.054440 49.686611 \n", "3 Indian Ocean Bay of Bengal 1006558.6 122.349983 90.946880 \n", "4 Indian Ocean Arabian Sea 865012.6 81.286291 82.852275 \n", ".. ... ... ... ... ... \n", "70 Atlantic Ocean Adriatic Sea 73289.5 21.576098 8.425867 \n", "71 Atlantic Ocean North Sea 163122.3 48.247612 20.200160 \n", "72 Atlantic Ocean Mediterranean Sea 96637.6 28.004916 11.197342 \n", "73 Atlantic Ocean North Atlantic 8277.8 9.270847 1.545336 \n", "74 Atlantic Ocean Baltic Sea 40545.2 28.919025 8.685652 \n", "\n", " RGI_AREA geometry \n", "0 171.941 POLYGON ((151.725 70.975, 151.725 70.97083, 15... \n", "1 763.493 POLYGON ((91.75 57.70417, 91.74542 57.70299, 9... \n", "2 10528.496 POLYGON ((97.76667 28.77083, 97.76335 28.77168... \n", "3 7906.081 MULTIPOLYGON (((88.02555 21.5715, 88.02361 21.... \n", "4 27206.546 MULTIPOLYGON (((67.44167 23.97006, 67.44028 23... \n", ".. ... ... \n", "70 313.417 POLYGON ((9.68333 46.425, 9.68418 46.42168, 9.... \n", "71 336.922 POLYGON ((4.45638 52.14527, 4.46175 52.1411, 4... \n", "72 917.302 POLYGON ((6.84167 47.82083, 6.84167 47.81667, ... \n", "73 972.419 POLYGON ((-17.5 64.625, -17.50035 64.60591, -1... \n", "74 34.257 POLYGON ((23.69167 68.84167, 23.70214 68.84102... \n", "\n", "[75 rows x 10 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd_basin_num\n", "# we will name them after their MRBID number\n", "# here you can also get the initial glacierized area by dividing RGI_AREA through AREA_CALC" ] }, { "cell_type": "markdown", "id": "7a34da1f-06ae-409e-b7c5-d2e4fb705fe5", "metadata": {}, "source": [ "check how many glaciers there are for an example basin:" ] }, { "cell_type": "code", "execution_count": 4, "id": "2a3e3b7a-3887-4c6c-bc69-ea96e7b3ba28", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RHONE , MRBID: 6243, n=1169\n" ] } ], "source": [ "basin = 'RHONE' #'RHONE' # INDUS\n", "basin_idx = pd_basin_num[pd_basin_num['RIVER_BASI'] == basin]['MRBID'].values[0]\n", "f = open('/home/www/fmaussion/misc/magicc/rgi_ids_per_basin.json')\n", "rgis_basin = json.load(f)[str(basin_idx)]\n", "odf_basin = odf.loc[rgis_basin]\n", "print(basin, f', MRBID: {basin_idx},', f'n={len(rgis_basin)}')\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "2141e79a-c24d-4ec2-899f-cc5369807b01", "metadata": {}, "outputs": [], "source": [ "# choose only relevant variables for the basin aggregation \n", "variables = ['volume','area', \n", " 'melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier',\n", " 'melt_off_glacier_monthly', 'melt_on_glacier_monthly', 'liq_prcp_off_glacier_monthly', 'liq_prcp_on_glacier_monthly']" ] }, { "cell_type": "code", "execution_count": 12, "id": "10566cbe-eac7-4758-8bbc-d3e17587109a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6101\n", "6104\n", "6110\n", "6202\n", "6209\n", "6213\n", "6219\n", "6223\n", "6227\n", "6237\n", "6241\n", "6242\n", "6243\n", "6254\n", "6255\n" ] } ], "source": [ "creation_date = strftime(\"%Y-%m-%d %H:%M:%S\", gmtime()) # here add the current time for info\n", "\n", "bc = '_bc_2000_2019' # '_bc_1979_2014'\n", "end_yr = '2100' # you could change that to '2300'\n", "cmip = 'CMIP6' # or change that to CMIP5\n", "\n", "path = f'/home/www/oggm/oggm-standard-projections/oggm_v16/2023.3/{cmip}/{end_yr}' \n", "\n", "merge = True #True\n", "if merge:\n", " utils.mkdir(f'{path}/basins')\n", " for basin in pd_basin_num.RIVER_BASI: \n", " # get the RGI ids of the glaciers of the chosen basin\n", " pd_basin_num = gpd.read_file('/home/www/fmaussion/misc/magicc/basins_shape/glacier_basins.shp')\n", " basin_idx = pd_basin_num[pd_basin_num['RIVER_BASI'] == basin]['MRBID'].values[0]\n", "\n", " print(basin_idx)\n", " \n", " area_basin_total = pd_basin_num[pd_basin_num['RIVER_BASI'] == basin]['AREA_CALC'].values[0]\n", " area_basin_glacier = pd_basin_num[pd_basin_num['RIVER_BASI'] == basin]['RGI_AREA'].values[0]\n", " f = open('/home/www/fmaussion/misc/magicc/rgi_ids_per_basin.json')\n", " rgis_basin = json.load(f)[str(basin_idx)]\n", " odf_basin = odf.loc[rgis_basin]\n", " utils.mkdir(f'{path}/basins/{basin_idx}')\n", "\n", " for hist in ['gcm_from_2000', 'w5e5_gcm_merged']:\n", " for GCM in ['CESM2-WACCM','CESM2','EC-Earth3', 'EC-Earth3-Veg']: #all_GCM: # loop through all GCMs\n", " for scen in all_scenario: # loop through all SSPs\n", " try: # check if GCM, SCENARIO combination exists\n", " rid = '_{}_{}'.format(GCM, scen) # put together the same filesuffix which was used during the projection runs\n", " ds_all= []\n", " for rgi_reg in odf_basin.O1Region.unique():\n", " odf_basin_reg = odf_basin.loc[odf_basin.O1Region==rgi_reg]\n", " rgis_basin_rgi_reg = odf_basin_reg.index.values\n", " with xr.open_mfdataset(f'{path}/RGI{rgi_reg}/run_hydro_{hist}_{GCM}_{scen}{bc}_Batch*.nc') as ds_tmp:\n", " ds_tmp = ds_tmp[variables]\n", " ds_tmp = ds_tmp.sel(rgi_id=rgis_basin_rgi_reg).load()\n", " ds_tmp['runoff_monthly'] = (ds_tmp['melt_off_glacier_monthly'] + ds_tmp['melt_on_glacier_monthly']\n", " + ds_tmp['liq_prcp_off_glacier_monthly'] + ds_tmp['liq_prcp_on_glacier_monthly'])\n", " ds_tmp['runoff'] = (ds_tmp['melt_off_glacier'] + ds_tmp['melt_on_glacier']\n", " + ds_tmp['liq_prcp_off_glacier'] + ds_tmp['liq_prcp_on_glacier'])\n", " # for the moment, only save these variables ...\n", " ds_tmp = ds_tmp[['runoff', 'runoff_monthly', 'volume', 'area']]\n", " ds_tmp.coords['gcm'] = GCM # add GCM as a coordinate\n", " ds_tmp.coords['gcm'].attrs['description'] = 'used global circulation model' # add a description for GCM\n", " ds_tmp = ds_tmp.expand_dims(\"gcm\") # add GCM as a dimension to all Data variables\n", " ds_tmp.coords['scenario'] = scen # add scenario (here ssp) as a coordinate\n", " ds_tmp.coords['scenario'].attrs['description'] = 'used scenario (here SSPs)'\n", " ds_tmp.coords['bias_correction'] = bc[1:]\n", " ds_tmp.coords['basin_idx'] = basin_idx \n", " ds_tmp.coords['basin_name'] = basin \n", " ds_tmp.coords['glaciers_in_basin'] = len(rgis_basin)\n", " ds_tmp.coords['glacierised_area_perc'] = 100*area_basin_glacier/area_basin_total\n", " ds_tmp = ds_tmp.expand_dims(\"scenario\") # add SSO as a dimension to all Data variables\n", " ds_tmp.attrs['creation_date'] = creation_date # also add todays date for info\n", " ds_tmp.runoff.attrs = {'unit':'kg yr-1', \n", " 'description':'Annual glacier runoff: sum of annual melt and liquid precipitation on and off the glacier using a fixed-gauge with a glacier minimum reference area from year 2000'}\n", " ds_tmp.runoff_monthly.attrs = {'unit':'kg month-1',\n", " 'description':'Monthly glacier runoff from sum of monthly melt and liquid precipitation on and off the glacier using a fixed-gauge with a glacier minimum reference area from year 2000'}\n", " ds_all.append(ds_tmp) # add the dataset with extra coordinates to our final ds_all array\n", " \n", " ds_merged = xr.concat(ds_all, dim='rgi_id') \n", " ds_merged.to_netcdf(f'{path}/basins/{basin_idx}/basin_{basin_idx}_run_hydro_{hist}{bc}{rid}.nc')\n", " \n", " except OSError as err: # here we land if an error occured\n", " if str(err) == 'no files to open': # This is the error message if the GCM, SCENARIO (here ssp) combination does not exist\n", " pass #print(f'No data for gcm {GCM} with scenario {scen} found!') # print a descriptive message\n", " else:\n", " raise OSError(err) # if an other error occured we just raise it\n", " # do the actual merging\n", " # this creates too large files \n", " #ds_merged = xr.combine_by_coords(ds_all, fill_value=np.nan) # define how the missing GCM, SCENARIO combinations should be filled \n", " #ds_merged.runoff.attrs = {'unit':'kg yr-1',\n", " # 'description':'Annual glacier runoff: sum of annual melt and liquid precipitation on and off the glacier, minimum reference area: 2000'}\n", " #ds_merged.runoff_monthly.attrs = {'unit':'kg month-1',\n", " # 'description':'Monthly glacier runoff: sum of monthly melt and liquid precipitation on and off the glacier, minimum reference area: 2000'}\n", " #ds_merged.to_netcdf(f'{path}/basins/basin_{basin_idx}_combined_run_hydro_{hist}_endyr2100_CMIP6{bc}.nc')" ] }, { "cell_type": "markdown", "id": "5b07bda4-0e55-4c4d-a180-fd8ef71c248c", "metadata": {}, "source": [ "--> an example how to load and plot e.g. the annual runoff for a single basin is given in https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~oggm/oggm-standard-projections/analysis_notebooks/workflow_to_analyse_per_glacier_projection_files.ipynb?flush_cache=true " ] }, { "cell_type": "markdown", "id": "eaead815-8d57-4ac0-a7dd-035cfe7ba206", "metadata": {}, "source": [ "**Check if all basins are there**" ] }, { "cell_type": "code", "execution_count": 22, "id": "756e430f-c86e-400e-9351-66ff9e47e485", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0\n", ".ipynb_checkpoints 0\n" ] } ], "source": [ "import os\n", "\n", "def count_folders_and_files(directory):\n", " num_folders = 0\n", " num_files = 0\n", "\n", " for _, dirs, files in os.walk(directory):\n", " num_folders += len(dirs)\n", " num_files += len(files)\n", "\n", " return num_folders, num_files\n", "\n", "directory = f'{path}/basins/'\n", "folders, files = count_folders_and_files(directory)\n", "\n", "#ssert folders==75\n", "gcm_ssp_combinations = len(gcms_cmip6)/2 # / 2 because pr/tas values\n", "assert files==gcm_ssp_combinations*2*75 # *2 because there are two historical options\n", "\n", "def count_files_per_folder(directory):\n", " folder_counts = {}\n", "\n", " for root, dirs, files in os.walk(directory):\n", " folder_name = os.path.basename(root)\n", " folder_counts[folder_name] = len(files)\n", "\n", " return folder_counts\n", "\n", "file_counts = count_files_per_folder(directory)\n", "\n", "for folder, count in file_counts.items():\n", " try:\n", " assert count==int(gcm_ssp_combinations*2) # *2 because there are two historical options\n", " except:\n", " print(folder,count)" ] }, { "cell_type": "code", "execution_count": null, "id": "545978d1-4ae6-4e42-b632-1e86667f1458", "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 }