wgms_mbprofiles_constant_dh.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "<img src=\"https://raw.githubusercontent.com/OGGM/oggm/master/docs/_static/logo.png\" width=\"40%\"  align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Creating constant step-size mass-balance profiles via interpolation for OGGM \n",
    "\n",
    "In this notebook, we use the created selected link file together with the raw MB profile data to create mass-balance profiles with a constant step size (i.e. altitude difference for all profile levels, here: dh = 50m) via interpolation. \n",
    "- Those MB profiles of a year for glaciers with just one measurement are removed as we can not create a profile with just a single measurement. \n",
    "- in addition, two glaciers had single measurements for some years with very large differences between Upper and Lower BOUNDARY (>150m) while they had other measurements with lower differences in that MB elevation altitudinal band. These measurements were discarded for the MB profiles!\n",
    "\n",
    "This notebook was initially created by Moritz Oberrauch and was updated and modified by Lilian Schuster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Read the WGMS files\n",
    "\n",
    "Change paths to necesary folders and data files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just download the newest data, change the path_to_download_data and the year and month accordingly. If you run the entire notebook, the new WGMS MB data should be processed for OGGM\n",
    "year = '2021'\n",
    "month = '05'\n",
    "path_to_download_data = '/home/lilianschuster/Downloads/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# working directory\n",
    "idir = f'{path_to_download_data}DOI-WGMS-FoG-{year}-{month}'\n",
    "# WMGS lookup table\n",
    "#df_links = pd.read_csv(os.path.join(idir, f'WGMS-FoG-{year}-{month}-AA-GLACIER-ID-LUT.csv'), encoding='iso8859_15')\n",
    "# WGMS FoG data file\n",
    "df_mb_all = pd.read_csv(os.path.join(idir, f'WGMS-FoG-{year}-{month}-EE-MASS-BALANCE.csv'), encoding='iso8859_15')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# output directory\n",
    "odir = '/home/lilianschuster/oggm-sample-data/wgms/'\n",
    "#odir = '/Users/oberrauch/Documents/Studium/work/WGMS_Profiles'\n",
    "df_links_sel = pd.read_csv(os.path.join(odir, 'rgi_wgms_links_20220112.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics about the height levels\n",
    "In order to choose one standardized step size for all mb profile, we'll take a closer look..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get bounds and mid levels for all glaciers and years\n",
    "Iterate over all (profiled) mass balance records and get the different bounds for each glacier and every year. Store them in multi-dimenstional dictionaries (bounds, mid_level) with wgms id and years as primary and secondary keys, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# create dictionaries as containers\n",
    "bounds = dict()\n",
    "mid_levels = dict()\n",
    "\n",
    "# iterate over all selceted glaciers\n",
    "for rid, wid in zip(df_links_sel.RGI60_ID, df_links_sel.WGMS_ID):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI60_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "\n",
    "    # create empty sets to store profile altitudes without duplicates\n",
    "    bd = set() # (lower and upper) bounds\n",
    "    ml = set() # mid levels\n",
    "    \n",
    "    # create containers\n",
    "    bd_ = dict()\n",
    "    ml_ = dict()\n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # save bounds for the current year\n",
    "        tmp = np.array(list(bd))\n",
    "        tmp.sort()\n",
    "        bd_[str(yr)] = tmp\n",
    "        \n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids += df_mb_sel_yr.UPPER_BOUND.values[:len(mids)]\n",
    "        mids *= 0.5\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids:\n",
    "            ml.add(int(m))\n",
    "        \n",
    "        # save mid levels for the current year\n",
    "        tmp = np.array(list(ml))\n",
    "        tmp.sort()\n",
    "        ml_[str(yr)] = tmp\n",
    "        \n",
    "    \n",
    "    # convert, sort and store set in dictionaries\n",
    "    bounds[str(wid)] = bd_\n",
    "    mid_levels[str(wid)] = ml_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate step size\n",
    "Now calculate the \"step size\" (i.e. altitude difference for all profile levels) and perform some basics statistics (median, min/max, ...) stored in a DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# empty DataFrame as container\n",
    "df_stat = pd.DataFrame()\n",
    "\n",
    "# iterate over every glacier and every year\n",
    "for wid in bounds.keys():\n",
    "    # subset for current glacier\n",
    "    bounds_gl = bounds[wid]\n",
    "    \n",
    "    # empty DataFrame for current glacier\n",
    "    df_stat_gl = pd.DataFrame()\n",
    "    \n",
    "    # iterate over all years\n",
    "    for yr in bounds_gl.keys():\n",
    "        \n",
    "        # subset for current year\n",
    "        bounds_gl_yr = bounds_gl[yr]\n",
    "        # compute step size\n",
    "        difs = bounds_gl_yr[1:] - bounds_gl_yr[:-1]\n",
    "        \n",
    "        # compute some statistical properties\n",
    "        df_stat_yr = pd.DataFrame({'min': difs.min(),\n",
    "                                'max': difs.max(),\n",
    "                                'median': int(round(np.median(difs))),\n",
    "                                'mean' : int(round(difs.mean())),\n",
    "                                'std' : int(round(difs.std()))},\n",
    "                                columns=list(['min', 'max', 'median', 'mean', 'std']),\n",
    "                                index=[yr])\n",
    "\n",
    "\n",
    "        # concat all years into one DataFrame\n",
    "        if (df_stat_gl.empty):\n",
    "            df_stat_gl = df_stat_yr.copy()\n",
    "        else: \n",
    "            df_stat_gl = pd.concat([df_stat_gl, df_stat_yr], axis=0)\n",
    "    \n",
    "    # sort by years\n",
    "    df_stat_gl = df_stat_gl.sort_index()\n",
    "    # create multi level index with WGMS ID and year\n",
    "    index = pd.MultiIndex.from_product([[wid],df_stat_gl.index.values], names=['wid', 'year'])\n",
    "    df_stat_gl.index = index\n",
    "    \n",
    "    # concat all glaciers into one DataFrame\n",
    "    if (df_stat.empty):\n",
    "        df_stat = df_stat_gl.copy()\n",
    "    else: \n",
    "        df_stat = pd.concat([df_stat, df_stat_gl], axis=0)\n",
    "    \n",
    "    df_stat = df_stat.sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data saved in a multi indexed DataFrame for each glacier (identified by WGMS Id) and each year, sorted in ascending order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_stat.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a first look at the overall outcome of the step sizes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stat.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the outliers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_stat[df_stat['min'] > 500]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The outliers with large stepsizes are STORBREEN, HOFSJOKULL (SW, E, N):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_links_sel[df_links_sel['WGMS_ID']==302]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all['WGMS_ID'] == 302][df_mb_all['LOWER_BOUND'] != 9999]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all['WGMS_ID'] == 3090][df_mb_all['LOWER_BOUND'] != 9999]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all['WGMS_ID'] == 3088][df_mb_all['LOWER_BOUND'] != 9999]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all['WGMS_ID'] == 3089][df_mb_all['LOWER_BOUND'] != 9999]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems like this profile covers the whole glacier. I'd suggest to check for 'single-profile' glaciers in a prior (or further) step and eliminate them!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excorsion: Single mb profiles - error in wgms data\n",
    "\n",
    "Now, since we know that these lines are errors in the wgms data I'll search for all of them and report the issue...\n",
    "\n",
    "-> Remove all measurements when there is only one measurement for that year and glacier, because then we can not use it for MB profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# create empty container to store indices of glaciers with just one profile\n",
    "single_ind = list()\n",
    "# iterate over all selceted glaciers\n",
    "for rid, wid in zip(df_links_sel.RGI60_ID, df_links_sel.WGMS_ID):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        df_mb_sel_yr = df_mb_sel[df_mb_sel.YEAR == yr]\n",
    "        if len(df_mb_sel_yr) == 1:\n",
    "            single_ind.append([wid, yr])\n",
    "            continue\n",
    "        \n",
    "# convert list into DataFrame\n",
    "single_ind = pd.DataFrame(single_ind, columns=['WGMS_ID', 'YEAR'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_stat.loc[str(302)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print DataFrame\n",
    "single_ind.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select information for single mb-profile glaciers from look up table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# create list as container\n",
    "df_mb_single = list()\n",
    "# iterate over all found wgms ids and years\n",
    "for _,si in single_ind.iterrows():\n",
    "    row_tmp = df_links_sel[df_links_sel.WGMS_ID == si.WGMS_ID]\n",
    "    df_mb_single.append(row_tmp.values[0])        \n",
    "\n",
    "# store in DataFrame\n",
    "df_mb_single = pd.DataFrame(df_mb_single, columns=row_tmp.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# print Data Frame\n",
    "df_mb_single.WGMS_ID.unique()\n",
    "#for si in df_mb_single.WGMS_ID.unique():\n",
    "#    print(df_mb_all[df_mb_all.WGMS_ID == si][df_mb_all[df_mb_all.WGMS_ID == si].LOWER_BOUND!=9999])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the mb table for some of those glaciers with only single-measurements for one year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all.WGMS_ID == 357][df_mb_all[df_mb_all.WGMS_ID == 357].LOWER_BOUND!=9999]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all.WGMS_ID == 3088][df_mb_all[df_mb_all.WGMS_ID == 3088].LOWER_BOUND!=9999]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all[df_mb_all.WGMS_ID == 302][df_mb_all[df_mb_all.WGMS_ID == 302].LOWER_BOUND!=9999]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the found 'one-profile' measurements for that specific year and glacier from the statistics data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only remove those years of the glaciers, not the entire glaciers!!!\n",
    "for _,si in single_ind[::-1].iterrows():\n",
    "    #midx = pd.MultiIndex(levels=['wid', 'year'],\n",
    "    #                     codes=[[str(si.WGMS_ID)],[yr]])\n",
    "    if (not df_stat.loc[str(si.WGMS_ID)].empty):\n",
    "            print('Single mb-profile glacier wgms-id ', si.WGMS_ID, ' of year {} deleted'.format(si.YEAR))\n",
    "            df_stat = df_stat.drop(index = (str(si.WGMS_ID), str(si.YEAR)))\n",
    "            #df_stat.drop(index=str(si.WGMS_ID), level=yr).copy()\n",
    "            #df_stat = df_stat.drop(midx).copy()\n",
    "            #print('t')\n",
    "    else:\n",
    "        print('Glacier wgms-id ', si.WGMS_ID, ' for the year {} already deleted or not in DataFrame!'.format(si.YEAR))\n",
    "        \n",
    "        \n",
    "#for _,si in single_ind[::-1].iterrows():\n",
    "#    try:\n",
    "#        if (not df_stat.loc[str(si.WGMS_ID)].empty):\n",
    "#            print('Single mb-profile glacier wgms-id ', si.WGMS_ID, ' deleted')\n",
    "#            df_stat = df_stat.drop(str(si.WGMS_ID)).copy()\n",
    "#            print('t')\n",
    "#    except:\n",
    "#        print('Glacier wgms-id ', si.WGMS_ID, ' already deleted or not in DataFrame!')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average step size\n",
    "Let's take a look at the average 'step size' over all glaciers...\n",
    "The mean is probaly less meaningfull than the median, even though the IQR seems to be quiet similar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# plot simple box plot\n",
    "df_stat[['mean','median']].plot.box()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The average step size is around 50 meter, and more than half of the step sizes are below 100 meters.\n",
    "The question is now, wheter to use 50 meters as step size and interpolate over quiet some altitude for a bunch of glaciers or use 100 meters as step size and loose some data..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calc max values of all mean and median values\n",
    "max_ = max(df_stat[['mean','median']].max())\n",
    "# set bin size and create bins\n",
    "bin_size = 5\n",
    "bins = np.arange(0,max_+bin_size, bin_size)\n",
    "# plot simple histogram\n",
    "df_stat[['mean','median']].plot.hist(stacked='True', bins=bins) #, normed=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram shows that there is a substatial number of glaciers and years with a step size around 50 meters and one around 100 meters. So let's take a look at the extrem values and foremost the maxima..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extrem values\n",
    "\n",
    "These show the smallest and biggest step size for every glacier and every year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# plot simple boxplot\n",
    "df_stat[['min','max']].plot.box()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# calc max value for mins and maxs\n",
    "max_ = max(df_stat[['min', 'max']].max())\n",
    "# set bin size and set bins\n",
    "bin_size = 10\n",
    "bins = np.arange(0,max_+bin_size, bin_size)\n",
    "# print simple histogram\n",
    "df_stat[['min', 'max']].plot.hist(stacked='True', bins=bins) #, normed=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of glaciers and years with step sizes over 100 is small and most are even around or below 50 meters. \n",
    "Given that the interpolation between bigger section is not a problem, but a big step size means more data loss we will use 50 meters. This can easily be changed in a second step..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some utility functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Linear interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lin_inter(f0,f1,h0,h1,h):\n",
    "    return f0 + (h-h0)/(h1-h0)*(f1-f0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find closest bounds or mid levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_bound(x, step_size=50):\n",
    "    return int(math.ceil(x/step_size))*step_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lower_mid(x, step_size=50):\n",
    "    if type(x) == np.ndarray:\n",
    "        return (x/step_size+0.5).astype(int)*step_size-step_size/2\n",
    "    return int(x/step_size+0.5)*step_size-step_size/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def upper_mid(x, step_size=50):\n",
    "    if type(x) == np.ndarray:\n",
    "        return ((x+step_size/2-1)/step_size).astype(int)*step_size+step_size/2\n",
    "    return int((x+step_size/2-1)/step_size)*step_size+step_size/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Actual routine to create the interpolated MB profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify the step size\n",
    "step_size = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_all.loc[df_mb_all.WGMS_ID == wid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Iterate over all glaciers\n",
    "df_mb_problems = []\n",
    "for rid, wid in zip(df_links_sel.RGI60_ID, df_links_sel.WGMS_ID):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set HAS_PROFILE flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI60_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "        \n",
    "    # save years in array\n",
    "    years = df_mb_sel.YEAR.unique()\n",
    "        \n",
    "    # create empty set to store profile mid level altitudes without duplicates\n",
    "    lb = set()\n",
    "    bd = set()\n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in years:\n",
    "        \n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        if wid == 942 or wid == 299:\n",
    "            # one probably false measurement with an upper / lower bound difference of more than 300m (2838 -> 3213)\n",
    "            #if len(df_mb_sel_yr.loc[df_mb_sel_yr.UPPER_BOUND - df_mb_sel_yr.LOWER_BOUND<200])>= 1:\n",
    "            df_mb_sel_yr = df_mb_sel_yr.loc[df_mb_sel_yr.UPPER_BOUND - df_mb_sel_yr.LOWER_BOUND<150]\n",
    "            \n",
    "            #print(wid, yr)\n",
    "        \n",
    "        \n",
    "        # check if more than one profile per year\n",
    "        if (len(df_mb_sel_yr) < 2):\n",
    "            # delete year from list of 'profiled' years\n",
    "            years = np.delete(years, np.argwhere(years==yr))\n",
    "            # print user message\n",
    "            print('Glacier wgms-id ', wid, ' has just one profile in ', yr)\n",
    "            # continue with next year\n",
    "            continue\n",
    "        \n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids_yr = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids_yr += df_mb_sel_yr.UPPER_BOUND.values[:len(mids_yr)]\n",
    "        mids_yr *= 0.5\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids_yr:\n",
    "            lb.add(int(m))\n",
    "            \n",
    "    # continue to next glacier if lb is emty\n",
    "    if not lb:\n",
    "        df_links_sel.loc[df_links_sel.RGI60_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "    \n",
    "    # create array with all possible mid levels (for a given step size)\n",
    "    mids_tot = np.arange(upper_mid(min(lb), step_size=step_size),\n",
    "                       lower_mid(max(lb), step_size=step_size)+step_size/2,\n",
    "                       step_size)\n",
    "    \n",
    "        \n",
    "    # create a DataFrame with one column for each mid level\n",
    "    # index by the year of the mb recors\n",
    "    prof = pd.DataFrame(np.nan, columns=mids_tot.astype(int), index=sorted(years))\n",
    "\n",
    "    # iterate over all years\n",
    "    for yr in prof.index:\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        if wid == 942 or wid == 299:\n",
    "            # one probably false measurement with an upper / lower bound difference of more than 300m (2838 -> 3213)\n",
    "            #if len(df_mb_sel_yr.loc[df_mb_sel_yr.UPPER_BOUND - df_mb_sel_yr.LOWER_BOUND<200])>= 1:\n",
    "            if len(df_mb_sel_yr.loc[df_mb_sel_yr.UPPER_BOUND - df_mb_sel_yr.LOWER_BOUND>=150])>0:\n",
    "                df_mb_problems.append(df_mb_sel_yr)\n",
    "            df_mb_sel_yr = df_mb_sel_yr.loc[df_mb_sel_yr.UPPER_BOUND - df_mb_sel_yr.LOWER_BOUND<150]\n",
    "        # check if more than one profile per year\n",
    "        # otherwise skip this year\n",
    "        if (len(df_mb_sel_yr) < 2):\n",
    "            continue\n",
    "            \n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids_yr = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids_yr += df_mb_sel_yr.UPPER_BOUND.values[:len(mids_yr)]\n",
    "        mids_yr *= 0.5\n",
    "        # create Series to store mb measurements at corresponding mid points for current year\n",
    "        prof_yr = pd.Series(df_mb_sel_yr.ANNUAL_BALANCE.values, index=mids_yr)\n",
    "        \n",
    "        # iterate over all mid levels\n",
    "        for mid in prof.columns:\n",
    "            \n",
    "            # mid level in DataFrame is mid level in Series\n",
    "            if (mid in prof_yr.index):\n",
    "                # mb value can be directly copied\n",
    "                # Value ERROR\n",
    "                prof.loc[yr,int(mid)] = prof_yr[mid]\n",
    "                #prof.loc[yr,int(mid)] = prof_yr[int(mid)]\n",
    "                # continue loop with next mid level\n",
    "                continue\n",
    "                \n",
    "            # check if there is a lower and upper mid level\n",
    "            if (any(prof_yr.index[prof_yr.index < mid]) & any(prof_yr.index[prof_yr.index > mid])):\n",
    "                low = prof_yr.index[prof_yr.index < mid][-1]\n",
    "                upp = prof_yr.index[prof_yr.index > mid][0]\n",
    "                # linear interpolation between closest bounds\n",
    "                prof.loc[yr,int(mid)] = lin_inter(prof_yr[low], prof_yr[upp], low, upp, mid)\n",
    "                continue\n",
    "                \n",
    "            # if mid level from DataFrame is bigger or smaller than any of the mid levels\n",
    "            # of the current year do nothing, since  extrapolation is not physical!\n",
    "            # eventually add lower and upper bounds to interpolation...\n",
    "                \n",
    "    # write DataFrame to *.csv file\n",
    "    prof.to_csv(os.path.join(odir, 'mb_profiles_constant_dh', 'profile_constant_dh_WGMS-{:05d}.csv'.format(wid)))\n",
    "    # set the HAS_PROFILE flag to TRUE for the current glacier\n",
    "    df_links_sel.loc[df_links_sel.RGI60_ID == rid, 'HAS_PROFILE'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those are the measurements that are removed because there are at the same time measurements at finer Upper-Lower bound scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_problems[0] # measurement 23514 is removed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_problems[1] # measurement 23530 is removed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mb_problems[2] # measurement 33812 is removed!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# select all glaciers with profile\n",
    "df_links_prof = df_links_sel[df_links_sel['HAS_PROFILE']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_links_prof[df_links_prof['RGI60_ID'] == 'RGI60-11.00897']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# randomly select one glacier\n",
    "#wid = df_links_prof.sample(1)['WGMS_ID'].values[0]\n",
    "wid = 491\n",
    "df_links_prof[df_links_prof['WGMS_ID'] == wid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read profiles\n",
    "prof_unsort = pd.read_csv(os.path.join(odir, 'mb_profiles', 'profile_WGMS-{:05d}.csv'.format(wid)), index_col=0)\n",
    "prof_sort = pd.read_csv(os.path.join(odir, 'mb_profiles_constant_dh', 'profile_constant_dh_WGMS-{:05d}.csv'.format(wid)), index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "prof_unsort.head() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "prof_sort.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = prof_sort.index\n",
    "n_yr = len(years)\n",
    "cols = 4 if n_yr > 8 else 2\n",
    "rows = math.ceil(n_yr/cols)\n",
    "f, axes = plt.subplots(rows, cols, sharex='col', sharey='row', figsize=[20,16])\n",
    "for i,axes_ in enumerate(axes):\n",
    "    for j,ax in enumerate(axes_):\n",
    "        loc = len(axes_)*i+j\n",
    "        if (loc < n_yr):\n",
    "            ax.plot(prof_unsort.loc[years[loc]].values, list(map(int, prof_unsort.loc[years[loc]].index)), color='k', lw=3, label='unsorted')\n",
    "            ax.plot(prof_sort.loc[years[loc]].values, list(map(int, prof_sort.loc[years[loc]].index)), color='r', lw=3, label='sorted')\n",
    "            x_min = ax.get_xlim()\n",
    "            ax.axvline(0, linestyle=':')\n",
    "            ax.text(0.5,0.9, years[loc], transform=ax.transAxes, horizontalalignment='center')\n",
    "        plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot MB - profiles per Latitude and Longitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# select all glaciers with profile\n",
    "df_links_prof = df_links_sel[df_links_sel['HAS_PROFILE']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create figure\n",
    "fig = plt.figure(figsize=[15,15])\n",
    "ax = fig.add_axes([0.1,0.1,0.8,0.8])\n",
    "\n",
    "# get all wgms ids\n",
    "wids = df_links_prof.WGMS_ID.values\n",
    "# iterate over all wgms ids\n",
    "for wid in wids:\n",
    "    if wid == 2660:\n",
    "        continue\n",
    "    # get profile from csv file\n",
    "    prof = pd.read_csv(os.path.join(odir, 'mb_profiles_constant_dh', 'profile_constant_dh_WGMS-{:05d}.csv'.format(wid)), index_col=0)\n",
    "    link = df_links_prof[df_links_prof.WGMS_ID == wid]\n",
    "    mb_mean = prof.mean(axis=0)\n",
    "    height = list(map(int, mb_mean.index.values))\n",
    "    mb = mb_mean.values\n",
    "    ax.plot(mb, height, label=wid)\n",
    "\n",
    "# add zero line\n",
    "ax.axvline(0, color='k', linestyle=':')\n",
    "# add axis labels\n",
    "ax.set_xlabel('mass balance [mm w.e.]')\n",
    "ax.set_ylabel('altitude [m a.s.l.]')\n",
    "# add label\n",
    "ax.legend(bbox_to_anchor=(1.35, 1.02), ncol=3, title='WGMS ID')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- we don't plot Claridenfirn (wgmsid 2660) because it has extremely negative values and would make the plot difficult to read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_links_sel[df_links_sel['WGMS_ID'] == 2660]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "prof = pd.read_csv(os.path.join(odir, 'mb_profiles_constant_dh', 'profile_constant_dh_WGMS-{:05d}.csv'.format(2660)), index_col=0)\n",
    "prof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## old stuff from Moritz\n",
    "### Actual routine - First try, didn't work out quite as planned... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IDEA: if there are more measurement points within one 'slice' of 100 meters, the interpolated mb value will be obtained by averaging the linear interpolation between the outer and inner slice..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# White glacier in Canada (WGMS ID: 0)\n",
    "df_links[df_links['WGMS_ID'] == 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_bound(x, step_size=100):\n",
    "    return int(math.ceil(x/step_size))*step_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_mid(x, step_size=100):\n",
    "    mid_size = step_size/2\n",
    "    return int(math.ceil(x/mid_size))*mid_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "step_size = 100\n",
    "\n",
    "# I just use the White glacier (WGMS ID 0) as test case\n",
    "for rid, wid in zip(['RGI50-03.04539'], [0]):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set HAS_PROFILE flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "        \n",
    "    # create empty set to store profile mid level altitudes without duplicates\n",
    "    lb = set()\n",
    "    bd = set()\n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids += df_mb_sel_yr.UPPER_BOUND.values[:len(mids)]\n",
    "        mids *= 0.5\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids:\n",
    "            lb.add(int(m))\n",
    "            \n",
    "    \n",
    "    print('All altitude steps:')\n",
    "    print(sorted(list(bd)))\n",
    "    \n",
    "    min_bd = min(bd)\n",
    "    max_bd = max(bd)\n",
    "    print('Min: '+str(min_bd))\n",
    "    print('Max: '+str(max_bd))\n",
    "    \n",
    "    tmp = np.arange(closest_bound(min_bd), closest_bound(max_bd)+step_size, step_size)\n",
    "    print(tmp)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works, but just for this one glacier WGMS ID 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "step_size = 100\n",
    "\n",
    "# I just use the White glacier (WGMS ID 0) as test cas\n",
    "for rid, wid in zip(['RGI50-03.04539'], [2665]):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set HAS_PROFILE flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "        \n",
    "    # create empty set to store profile mid level altitudes without duplicates\n",
    "    lb = set()\n",
    "    bd = set()\n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids += df_mb_sel_yr.UPPER_BOUND.values[:len(mids)]\n",
    "        mids *= 0.5\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids:\n",
    "            lb.add(int(m))\n",
    "            \n",
    "    # create array with all possible lower and upper bounds (for step size)\n",
    "    bounds = np.arange(closest_bound(min(bd), step_size), closest_bound(max(bd), step_size)+step_size, step_size)\n",
    "    # calculate mid levels\n",
    "    mids = (bounds[:-1] + bounds[1:])/2\n",
    "    \n",
    "    # create a DataFrame with one column for each mid level\n",
    "    # index by the year of the mb recors\n",
    "    prof = pd.DataFrame(np.nan, columns=mids, index=sorted(df_mb_sel.YEAR.unique()))\n",
    "    \n",
    "    # iterate over all years\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        \n",
    "        # iterate over all mid levels\n",
    "        for mid in prof.columns:\n",
    "            # fetch levels above and below the mid level altitude on which a linear interpolation can be performed\n",
    "            lower = df_mb_sel_yr[df_mb_sel_yr['LOWER_BOUND'] < mid][df_mb_sel_yr['UPPER_BOUND'] > (mid-step_size/2)]\n",
    "            upper = df_mb_sel_yr[df_mb_sel_yr['UPPER_BOUND'] > mid][df_mb_sel_yr['LOWER_BOUND'] < (mid+step_size/2)]\n",
    "            # number of slices found\n",
    "            l_low = len(lower)\n",
    "            l_up = len(upper)\n",
    "            \n",
    "            # if there are no lower or upper slices, skip this mid level and continue with loop\n",
    "            if(l_low == 0 | l_up == 0):\n",
    "                continue\n",
    "            \n",
    "            # check if same levels above and below mid level\n",
    "            if(l_low == l_up):\n",
    "                # containers for interpolated mb and distance as weight for the average\n",
    "                mb = np.zeros(l_low)\n",
    "                weight = np.zeros(l_low)\n",
    "                # iterate over all steps\n",
    "                for i in np.arange(0,l_low):\n",
    "                    # get value\n",
    "                    mb_low = lower.iloc[i]['ANNUAL_BALANCE']\n",
    "                    bd_low = lower.iloc[i]['LOWER_BOUND']\n",
    "                    mb_up = lower.iloc[(l_up-1-i)]['ANNUAL_BALANCE']\n",
    "                    bd_up = lower.iloc[i]['UPPER_BOUND']\n",
    "                    # interpolate linearaly\n",
    "                    mb[i] = lin_inter(mb_low, mb_up, bd_low, bd_up, mid)\n",
    "                    # calculate distance from mid point\n",
    "                    weight[i] = step_size/2 - abs((bd_up+bd_low)/2-mid)\n",
    "                    \n",
    "                # average the mb inversely weighted with distance from mid level\n",
    "                mb_avg = np.average(mb, weights=weight)\n",
    "                \n",
    "                # save in DataFrame \n",
    "                prof.loc[yr, int(mid)] = mb_avg\n",
    "                \n",
    "    # write DataFrame to *.csv file\n",
    "    prof.to_csv(os.path.join(odir, 'profiles_mo', 'profile_WGMS-{:05d}.csv'.format(wid)))\n",
    "    # set the HAS_PROFILE flag to TRUE for the current glacier\n",
    "    df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_bound(x, step_size=100):\n",
    "    return int(round(x/step_size))*step_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_mid(x, step_size=100):\n",
    "    mid_size = step_size/2\n",
    "    return int(round(x/mid_size))*mid_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "step_size = 50\n",
    "\n",
    "# I just use the White glacier (WGMS ID 0) as test cas\n",
    "for rid, wid in zip(['RGI50-03.04539'], [2665]):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set HAS_PROFILE flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "        \n",
    "    # create empty set to store profile mid level altitudes without duplicates\n",
    "    lb = set()\n",
    "    bd = set()\n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids_ = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids_ += df_mb_sel_yr.UPPER_BOUND.values[:len(mids_)]\n",
    "        mids_ *= 0.5\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids_:\n",
    "            lb.add(int(m))\n",
    "            \n",
    "    # create array with all possible lower and upper bounds (for step size)\n",
    "    bounds = np.arange(closest_bound(min(bd), step_size), closest_bound(max(bd), step_size)+step_size, step_size)\n",
    "    # calculate mid levels\n",
    "    mids = (bounds[:-1] + bounds[1:])/2\n",
    "    \n",
    "    # create a DataFrame with one column for each mid level\n",
    "    # index by the year of the mb recors\n",
    "    prof = pd.DataFrame(np.nan, columns=mids.astype(int), index=sorted(df_mb_sel.YEAR.unique()))\n",
    "    \n",
    "    # iterate over all years\n",
    "    for yr in [2000]:\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        \n",
    "        # iterate over all mid levels\n",
    "        for mid in prof.columns:\n",
    "            print(mid)\n",
    "            # fetch levels above and below the mid level altitude on which a linear interpolation can be performed\n",
    "            slices = df_mb_sel_yr[df_mb_sel_yr['LOWER_BOUND'] <= mid][df_mb_sel_yr['UPPER_BOUND'] >= mid]\n",
    "            \n",
    "            print('Slices:')\n",
    "            print(slices[['LOWER_BOUND', 'UPPER_BOUND', 'ANNUAL_BALANCE']])\n",
    "\n",
    "            # number of slices\n",
    "            l_slices = len(slices)\n",
    "            \n",
    "            # just continue if there are slices around the mid level \n",
    "            if(l_slices != 0):\n",
    "                \n",
    "                # containers for interpolated mb and distance as weight for the average\n",
    "                mb = np.zeros(l_slices)\n",
    "                weight = np.zeros(l_slices)\n",
    "                # iterate over all steps\n",
    "                for i,s in enumerate(slices):\n",
    "                    # get values\n",
    "                    mb_ = s['ANNUAL_BALANCE']\n",
    "                    bd_low = s['LOWER_BOUND']\n",
    "                    bd_up = s['UPPER_BOUND']\n",
    "                    # interpolate linearaly\n",
    "                    mb[i] = lin_inter(mb_low, mb_up, bd_low, bd_up, mid)\n",
    "                    # calculate distance from mid point\n",
    "                    weight[i] = step_size/2 - abs((bd_up+bd_low)/2-mid)\n",
    "                    \n",
    "                # average the mb inversely weighted with distance from mid level\n",
    "                if (sum(weight)==0):\n",
    "                    weight = None\n",
    "                mb_avg = np.average(mb, weights=weight)\n",
    "                \n",
    "                # save in DataFrame \n",
    "                prof.loc[yr, int(mid)] = mb_avg\n",
    "                \n",
    "    # write DataFrame to *.csv file\n",
    "    prof.to_csv(os.path.join(odir, 'profiles_mo', 'profile_WGMS-{:05d}.csv'.format(wid)))\n",
    "    # set the HAS_PROFILE flag to TRUE for the current glacier\n",
    "    df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "[df_mb_sel_yr['UPPER_BOUND'] < (mid+step_size/2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "bd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "prof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sorted(list(lb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == 2000]\n",
    "mid=150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "lower = df_mb_sel_yr[df_mb_sel_yr['LOWER_BOUND'] < mid][df_mb_sel_yr['UPPER_BOUND'] > (mid-step_size/2)]\n",
    "print(lower)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "upper = df_mb_sel_yr[df_mb_sel_yr['UPPER_BOUND'] > mid][df_mb_sel_yr['LOWER_BOUND'] < (mid+step_size/2)]\n",
    "print(upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for special case, for that write in a file..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "step_size = 100\n",
    "\n",
    "\n",
    "# open file\n",
    "f = open('workfile', 'w')\n",
    "\n",
    "\n",
    "# iterate over all glaciers\n",
    "for rid, wid in zip(df_links_sel.RGI_ID, df_links_sel.WGMS_ID):\n",
    "    # copy the mb records for each glacier\n",
    "    df_mb_sel = df_mb_all.loc[df_mb_all.WGMS_ID == wid].copy()\n",
    "    # further just use glaciers with a 'profiled' mb record\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.LOWER_BOUND != 9999]\n",
    "    df_mb_sel = df_mb_sel.loc[df_mb_sel.UPPER_BOUND != 9999]\n",
    "    # if glacier has no profiles, set HAS_PROFILE flag to FALSE\n",
    "    # and continue with the next glacier\n",
    "    if len(df_mb_sel) == 0:\n",
    "        df_links_sel.loc[df_links_sel.RGI_ID == rid, 'HAS_PROFILE'] = False\n",
    "        continue\n",
    "        \n",
    "    # create empty set to store profile mid level altitudes without duplicates\n",
    "    lb = set()\n",
    "    bd = set()\n",
    "    \n",
    "    f.write('Glacier: '+str(wid)+'\\n')\n",
    "        \n",
    "    \n",
    "    # iterate over all (unique) years with profiled mb records\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        # calculate mean between lower and upper bounds\n",
    "        mids = df_mb_sel_yr.LOWER_BOUND.values*1.\n",
    "        mids += df_mb_sel_yr.UPPER_BOUND.values[:len(mids)]\n",
    "        mids *= 0.5\n",
    "        \n",
    "        # save all unique bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.LOWER_BOUND.values] # lower bounds\n",
    "        [bd.add(int(b)) for b in df_mb_sel_yr.UPPER_BOUND.values] # upper bounds\n",
    "        \n",
    "        # keep all unique mid levels\n",
    "        for m in mids:\n",
    "            lb.add(int(m))\n",
    "            \n",
    "    # create array with all possible lower and upper bounds (for step size)\n",
    "    bounds = np.arange(closest_bound(min_bd, step_size), closest_bound(max_bd, step_size)+step_size, step_size)\n",
    "    # calculate mid levels\n",
    "    mids = (bounds[:-1] + bounds[1:])/2\n",
    "    \n",
    "    # create a DataFrame with one column for each mid level\n",
    "    # index by the year of the mb recors\n",
    "    prof = pd.DataFrame(np.nan, columns=mids, index=sorted(df_mb_sel.YEAR.unique()))\n",
    "    \n",
    "    # iterate over all years\n",
    "    for yr in df_mb_sel.YEAR.unique():\n",
    "        # subset dataframe for each year\n",
    "        df_mb_sel_yr = df_mb_sel.loc[df_mb_sel.YEAR == yr]\n",
    "        \n",
    "        f.write('Year: '+str(yr)+'\\n')\n",
    "        \n",
    "        \n",
    "        # iterate over all mid levels\n",
    "        for mid in prof.columns:\n",
    "            \n",
    "            f.write('Mid: '+str(mid)+'\\t')\n",
    "            \n",
    "            # fetch levels above and below the mid level altitude on which a linear interpolation can be performed\n",
    "            lower = df_mb_sel_yr[df_mb_sel_yr['LOWER_BOUND'] < mid][df_mb_sel_yr['UPPER_BOUND'] > (mid-step_size/2)]\n",
    "            upper = df_mb_sel_yr[df_mb_sel_yr['UPPER_BOUND'] > mid][df_mb_sel_yr['LOWER_BOUND'] < (mid+step_size/2)]\n",
    "            # number of slices found\n",
    "            l_low = len(lower)\n",
    "            l_up = len(upper)\n",
    "            \n",
    "            # if there are no lower or upper slices, skip this mid level and continue with loop\n",
    "            if(l_low == 0 | l_up == 0):\n",
    "                f.write('No lower or upper slices\\n')\n",
    "                continue\n",
    "            \n",
    "            # check if same levels above and below mid level\n",
    "            if(l_low == l_up):\n",
    "                f.write('Same number of upper and lower slices\\n')\n",
    "            elif (l_low > l_up):\n",
    "                f.write('More lower than upper slices\\n')\n",
    "            else: # l_low < l_up\n",
    "                f.write(r'More upper than lower slices\\n')\n",
    "                \n",
    "        \n",
    "        f.write('\\n')\n",
    "    \n",
    "    f.write('--------------------------\\n')\n",
    "    \n",
    "    \n",
    "# close file\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_mid(x, step_size=100):\n",
    "    mid_size = step_size/2\n",
    "    return int(math.ceil(x/mid_size))*mid_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "closest_mid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Links: add some stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Handle various RGI versions\n",
    "df_links_sel.rename(columns = {'RGI_ID':'RGI50_ID'}, inplace = True)\n",
    "df_links_sel['RGI40_ID'] = df_links_sel['RGI50_ID']\n",
    "df_links_sel['RGI40_ID'] = [rid.replace('RGI50', 'RGI40') for rid in df_links_sel['RGI40_ID']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Get the RGI\n",
    "import geopandas as gpd\n",
    "import glob, os\n",
    "frgi = '/home/mowglie/Documents/rgi50_allglaciers.csv'\n",
    "if not os.path.exists(frgi):\n",
    "    # one time action only\n",
    "    fs = list(sorted(glob.glob(\"/home/mowglie/disk/Data/GIS/SHAPES/RGI/RGI_V5/*/*_rgi50_*.shp\")))[2:]\n",
    "    out = []\n",
    "    for f in fs:\n",
    "        sh = gpd.read_file(f).set_index('RGIId')\n",
    "        del sh['geometry']\n",
    "        del sh['GLIMSId']\n",
    "        del sh['Name']\n",
    "        out.append(sh)\n",
    "    mdf = pd.concat(out)\n",
    "    mdf.to_csv(frgi)\n",
    "mdf = pd.read_csv(frgi, index_col=0, converters={'GlacType': str, 'RGIFlag':str, 'BgnDate':str, \n",
    "                                                 'O1Region': str, 'O2Region':str})\n",
    "mdf['RGI_REG'] = [rid.split('-')[1].split('.')[0] for rid in mdf.index]\n",
    "# add region names\n",
    "sr = gpd.read_file('/home/mowglie/disk/Data/GIS/SHAPES/RGI/RGI_V5/00_rgi50_regions/00_rgi50_O1Regions.shp')\n",
    "sr = sr.drop_duplicates('Secondary_').set_index('Secondary_')[['Primary_ID']]\n",
    "sr['Primary_ID'] = [i + ': ' + s for i, s in sr.Primary_ID.iteritems()]\n",
    "mdf['RGI_REG_NAME'] = sr.loc[mdf.RGI_REG].Primary_ID.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Read glacier attrs\n",
    "key1 = {'0': 'Glacier',\n",
    "        '1': 'Ice cap',\n",
    "        '2': 'Perennial snowfield',\n",
    "        '3': 'Seasonal snowfield',\n",
    "        '9': 'Not assigned',\n",
    "        }\n",
    "\n",
    "key2 = {'0': 'Land-terminating',\n",
    "        '1': 'Marine-terminating',\n",
    "        '2': 'Lake-terminating',\n",
    "        '3': 'Dry calving',\n",
    "        '4': 'Regenerated',\n",
    "        '5': 'Shelf-terminating',\n",
    "        '9': 'Not assigned',\n",
    "        }\n",
    "\n",
    "def is_tidewater(ttype):\n",
    "    return \n",
    "\n",
    "mdf['GlacierType'] = [key1[gtype[0]] for gtype in mdf.GlacType]\n",
    "mdf['TerminusType'] = [key2[gtype[1]] for gtype in mdf.GlacType]\n",
    "mdf['IsTidewater'] = [ttype in ['Marine-terminating', 'Lake-terminating'] for ttype in mdf.TerminusType]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# add lons and lats and other attrs to the WGMS ones\n",
    "smdf = mdf.loc[df_links_sel.RGI50_ID]\n",
    "df_links_sel['CenLon'] = smdf.CenLon.values\n",
    "df_links_sel['CenLat'] = smdf.CenLat.values\n",
    "df_links_sel['GlacierType'] = smdf.GlacierType.values\n",
    "df_links_sel['TerminusType'] = smdf.TerminusType.values\n",
    "df_links_sel['IsTidewater'] = smdf.IsTidewater.values\n",
    "df_links_sel['RGI_REG_NAME'] = smdf.RGI_REG_NAME.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "df_links_sel = df_links_sel[['CenLon', 'CenLat',\n",
    "                             'POLITICAL_UNIT', 'NAME', 'WGMS_ID', 'PSFG_ID', 'WGI_ID', 'GLIMS_ID',\n",
    "                             'RGI40_ID', 'RGI50_ID', 'RGI_REG', 'RGI_REG_NAME', 'GlacierType', 'TerminusType', \n",
    "                             'IsTidewater', 'N_MB_YRS', 'HAS_PROFILE', 'REMARKS']]\n",
    "df_links_sel.to_csv(os.path.join(odir, 'rgi_wgms_links_20170217_RGIV5.csv'.format(wid)), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Some plots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.set_context('poster')\n",
    "sns.set_style('whitegrid')\n",
    "pdir = '/home/mowglie/Documents/git/fmaussion.github.io/images/blog/wgms-links'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_links_sel['N_MB_YRS'].plot(kind='hist', color='C2', bins=np.arange(21)*5);\n",
    "plt.xlim(5, 100);\n",
    "plt.ylabel('Number of glaciers')\n",
    "plt.xlabel('Length of the timeseries (years)');\n",
    "plt.tight_layout();\n",
    "plt.savefig(os.path.join(pdir, 'nglacier-hist.png'), dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "f = plt.figure(figsize=(12, 7))\n",
    "ax = plt.axes(projection=ccrs.Robinson())\n",
    "# mark a known place to help us geo-locate ourselves\n",
    "ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())\n",
    "ax.stock_img()\n",
    "ax.add_feature(cartopy.feature.COASTLINE);\n",
    "s = df_links_sel.loc[df_links_sel.N_MB_YRS < 10]\n",
    "print(len(s))\n",
    "ax.scatter(s.CenLon, s.CenLat, label='< 10 MB years', s=50,\n",
    "           edgecolor='k', facecolor='C0', transform=ccrs.PlateCarree(), zorder=99)\n",
    "s = df_links_sel.loc[(df_links_sel.N_MB_YRS >= 10) & (df_links_sel.N_MB_YRS < 30)]\n",
    "print(len(s))\n",
    "ax.scatter(s.CenLon, s.CenLat, label='$\\geq$ 10 and < 30 MB years', s=50,\n",
    "           edgecolor='k', facecolor='C1', transform=ccrs.PlateCarree(), zorder=99)\n",
    "s = df_links_sel.loc[df_links_sel.N_MB_YRS >= 30]\n",
    "print(len(s))\n",
    "ax.scatter(s.CenLon, s.CenLat, label='$\\geq$ 30 MB years', s=50,\n",
    "           edgecolor='k', facecolor='C2', transform=ccrs.PlateCarree(), zorder=99)\n",
    "plt.title('WGMS glaciers with at least 5 years of mass-balance data')\n",
    "plt.legend(loc=4, frameon=True)\n",
    "plt.tight_layout();\n",
    "plt.savefig(os.path.join(pdir, 'glacier-map.png'), dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df_links_sel.TerminusType.value_counts().to_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ax = sns.countplot(x='RGI_REG', hue=\"TerminusType\", data=df_links_sel);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "md = pd.concat([mdf.GlacierType.value_counts().to_frame(name='RGI V5').T, \n",
    "          df_links_sel.GlacierType.value_counts().to_frame(name='WGMS').T]\n",
    "          ).T\n",
    "md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "md = pd.concat([mdf.TerminusType.value_counts().to_frame(name='RGI V5').T, \n",
    "          df_links_sel.TerminusType.value_counts().to_frame(name='WGMS').T]\n",
    "          ).T\n",
    "md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "area_per_reg = mdf[['Area', 'RGI_REG_NAME']].groupby('RGI_REG_NAME').sum()\n",
    "area_per_reg['N_WGMS'] = df_links_sel.RGI_REG_NAME.value_counts()\n",
    "area_per_reg = area_per_reg.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sns.barplot(x=\"Area\", y=\"RGI_REG_NAME\", data=area_per_reg);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "area_per_reg['N_WGMS_PER_UNIT'] = area_per_reg.N_WGMS / area_per_reg.Area * 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sns.barplot(x=\"N_WGMS\", y=\"RGI_REG_NAME\", data=area_per_reg);  # , palette=sns.husl_palette(19, s=.7, l=.5)\n",
    "plt.ylabel('')\n",
    "plt.xlabel('')\n",
    "plt.title('Number of WGMS glaciers per RGI region');\n",
    "plt.tight_layout();\n",
    "plt.savefig(os.path.join(pdir, 'barplot-ng.png'), dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sns.barplot(x=\"N_WGMS_PER_UNIT\", y=\"RGI_REG_NAME\", data=area_per_reg);\n",
    "plt.ylabel('')\n",
    "plt.xlabel('')\n",
    "plt.title('Number of WGMS glaciers per 1,000 km$^2$ of ice, per RGI region');\n",
    "plt.tight_layout();\n",
    "plt.savefig(os.path.join(pdir, 'barplot-perice.png'), dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "nmb_yrs = df_links_sel[[\"RGI_REG\", 'N_MB_YRS']].groupby(\"RGI_REG\").sum()\n",
    "i = []\n",
    "for k, d in nmb_yrs.iterrows():\n",
    "     i.extend([k] * d.values[0])\n",
    "df = pd.DataFrame()\n",
    "df[\"RGI_REG\"] = i\n",
    "ax = sns.countplot(x=\"RGI_REG\", data=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}