{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10 minutes to... OGGM as an accelerator for machine learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we want to showcase what OGGM does best: **preparing data for your modelling workflow**.\n",
    "\n",
    "We use preprocessed directories which contain most data available in [the OGGM shop](https://docs.oggm.org/en/stable/input-data.html) to illustrate how these could be used to inform data-based workflows. The data that is available in the shop and is show cased here, is more than is required for the regular OGGM workflow, which you will see in a bit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tags:** beginner, shop, workflow "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessed directories with additional products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to use the South Glacier example taken from the [ITMIX experiment](https://www.the-cryosphere.net/11/949/2017/). It is a small (5.6 km2) glacier in Alaska."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Libs\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import salem\n",
    "\n",
    "# OGGM\n",
    "import oggm.cfg as cfg\n",
    "from oggm import utils, workflow, tasks, graphics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Initialize OGGM and set up the default run parameters\n",
    "cfg.initialize(logging_level='WARNING')\n",
    "cfg.PARAMS['use_multiprocessing'] = False\n",
    "# Local working directory (where OGGM will write its output)\n",
    "cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_Toy_Thickness_Model')\n",
    "# We use the preprocessed directories with additional data in it: \"W5E5_w_data\"  \n",
    "base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5_w_data/'\n",
    "gdirs = workflow.init_glacier_directories(['RGI60-01.16195'], from_prepro_level=3, prepro_base_url=base_url, prepro_border=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Pick our glacier\n",
    "gdir = gdirs[0]\n",
    "gdir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OGGM-Shop datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using here the preprocessed glacier directories that contain more data than the default ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:\n",
    "    ds = ds.load()\n",
    "# List all variables\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's already quite a lot! We have access to a bunch of data for this glacier, lets have a look. We prepare the map first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "smap = ds.salem.get_map(countries=False)\n",
    "smap.set_shapefile(gdir.read_shapefile('outlines'))\n",
    "smap.set_topography(ds.topo.data);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ITSLive velocity data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets start with the [ITSLIVE](https://its-live.jpl.nasa.gov/#data) velocity data: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the velocity data\n",
    "u = ds.itslive_vx.where(ds.glacier_mask)\n",
    "v = ds.itslive_vy.where(ds.glacier_mask)\n",
    "ws = ds.itslive_v.where(ds.glacier_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `.where(ds.glacier_mask)` command will remove the data outside of the glacier outline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the axes ready\n",
    "f, ax = plt.subplots(figsize=(9, 9))\n",
    "\n",
    "# Quiver only every N grid point\n",
    "us = u[1::3, 1::3]\n",
    "vs = v[1::3, 1::3]\n",
    "\n",
    "smap.set_data(ws)\n",
    "smap.set_cmap('Blues')\n",
    "smap.plot(ax=ax)\n",
    "smap.append_colorbar(ax=ax, label='ice velocity (m yr$^{-1}$)')\n",
    "\n",
    "# transform their coordinates to the map reference system and plot the arrows\n",
    "xx, yy = smap.grid.transform(us.x.values, us.y.values, crs=gdir.grid.proj)\n",
    "xx, yy = np.meshgrid(xx, yy)\n",
    "qu = ax.quiver(xx, yy, us.values, vs.values)\n",
    "qk = ax.quiverkey(qu, 0.82, 0.97, 10, '10 m yr$^{-1}$',\n",
    "                  labelpos='E', coordinates='axes')\n",
    "ax.set_title('ITS-LIVE velocity');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Millan 2022 velocity data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is more velocity data. Here follows the [Millan et al. (2022)](https://doi.org/10.1038/s41561-021-00885-z) dataset: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the velocity data\n",
    "u = ds.millan_vx.where(ds.glacier_mask)\n",
    "v = ds.millan_vy.where(ds.glacier_mask)\n",
    "ws = ds.millan_v.where(ds.glacier_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the axes ready\n",
    "f, ax = plt.subplots(figsize=(9, 9))\n",
    "\n",
    "# Quiver only every N grid point\n",
    "us = u[1::3, 1::3]\n",
    "vs = v[1::3, 1::3]\n",
    "\n",
    "smap.set_data(ws)\n",
    "smap.set_cmap('Blues')\n",
    "smap.plot(ax=ax)\n",
    "smap.append_colorbar(ax=ax, label='ice velocity (m yr$^{-1}$)')\n",
    "\n",
    "# transform their coordinates to the map reference system and plot the arrows\n",
    "xx, yy = smap.grid.transform(us.x.values, us.y.values, crs=gdir.grid.proj)\n",
    "xx, yy = np.meshgrid(xx, yy)\n",
    "qu = ax.quiver(xx, yy, us.values, vs.values)\n",
    "qk = ax.quiverkey(qu, 0.82, 0.97, 10, '10 m yr$^{-1}$',\n",
    "                  labelpos='E', coordinates='axes')\n",
    "ax.set_title('Millan 2022 velocity');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thickness data from Farinotti 2019 and Millan 2022 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also provide ice thickness data from [Farinotti et al. (2019)](https://doi.org/10.1038/s41561-019-0300-3) and [Millan et al. (2022)](https://doi.org/10.1038/s41561-021-00885-z) in the same gridded format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the axes ready\n",
    "f, ax = plt.subplots(figsize=(9, 9))\n",
    "smap.set_cmap('viridis')\n",
    "smap.set_data(ds.consensus_ice_thickness)\n",
    "smap.plot(ax=ax)\n",
    "smap.append_colorbar(ax=ax, label='ice thickness (m)')\n",
    "ax.set_title('Farinotti 2019 thickness');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get the axes ready\n",
    "f, ax = plt.subplots(figsize=(9, 9))\n",
    "smap.set_cmap('viridis')\n",
    "smap.set_data(ds.millan_ice_thickness.where(ds.glacier_mask))\n",
    "smap.plot(ax=ax)\n",
    "smap.append_colorbar(ax=ax, label='ice thickness (m)')\n",
    "ax.set_title('Millan 2022 thickness');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some additional gridded attributes "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also add some attributes that OGGM can compute for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Tested tasks\n",
    "task_list = [\n",
    "    tasks.gridded_attributes,\n",
    "    tasks.gridded_mb_attributes,\n",
    "]\n",
    "for task in task_list:\n",
    "    workflow.execute_entity_task(task, gdirs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's open the gridded data file again with xarray:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:\n",
    "    ds = ds.load()\n",
    "# List all variables\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file contains several new variables with their description. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds.oggm_mb_above_z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot a few of them (we show how to plot them with xarray and with oggm):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds.slope.plot();\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))\n",
    "graphics.plot_raster(gdir, var_name='aspect', cmap='twilight', ax=ax1)\n",
    "graphics.plot_raster(gdir, var_name='oggm_mb_above_z', ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>\n",
    "        In a few lines of code, we have used OGGM to generate or deliver a bunch of data for this glaciers. A similar workflow can be used on ALL of them! With this, we hope to facilitate access to data for many other models or machine learning workflows.\n",
    "    </b>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Not convinced yet? Lets spend 10 more minutes to apply a (very simple) machine learning workflow "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve these attributes at point locations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this glacier, we have ice thickness observations (all downloaded from the supplementary material of the [ITMIX paper](https://www.the-cryosphere.net/11/949/2017/)). Let's make a table out of it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = salem.read_shapefile(utils.get_demo_file('IceThick_SouthGlacier.shp'))\n",
    "coords = np.array([p.xy for p in df.geometry]).squeeze()\n",
    "df['lon'] = coords[:, 0]\n",
    "df['lat'] = coords[:, 1]\n",
    "df = df[['lon', 'lat', 'thick']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the longitudes and latitudes to the glacier map projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xx, yy = salem.transform_proj(salem.wgs84, gdir.grid.proj, df['lon'].values, df['lat'].values)\n",
    "df['x'] = xx\n",
    "df['y'] = yy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot these data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "geom = gdir.read_shapefile('outlines')\n",
    "f, ax = plt.subplots()\n",
    "df.plot.scatter(x='x', y='y', c='thick', cmap='viridis', s=10, ax=ax);\n",
    "geom.plot(ax=ax, facecolor='none', edgecolor='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 1: interpolation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The measurement points of this dataset are very frequent and close to each other. There are plenty of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will keep them all and interpolate the variables of interest at the point's location. We use [xarray](http://xarray.pydata.org/en/stable/interpolation.html#advanced-interpolation) for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vns = ['topo',\n",
    "       'slope',\n",
    "       'slope_factor',\n",
    "       'aspect',\n",
    "       'dis_from_border',\n",
    "       'catchment_area',\n",
    "       'lin_mb_above_z',\n",
    "       'oggm_mb_above_z',\n",
    "       'millan_v',\n",
    "       'itslive_v',\n",
    "       ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Interpolate (bilinear)\n",
    "for vn in vns:\n",
    "    df[vn] = ds[vn].interp(x=('z', df.x), y=('z', df.y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how these variables can explain the measured ice thickness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))\n",
    "df.plot.scatter(x='dis_from_border', y='thick', ax=ax1); ax1.set_title('dis_from_border');\n",
    "df.plot.scatter(x='slope', y='thick', ax=ax2); ax2.set_title('slope');\n",
    "df.plot.scatter(x='oggm_mb_above_z', y='thick', ax=ax3); ax3.set_title('oggm_mb_above_z');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a negative correlation with slope (as expected), a positive correlation with the mass-flux (oggm_mb_above_z), and a slight positive correlation with the distance from the glacier boundaries. There is also some correlaction with ice velocity, but not a strong one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n",
    "df.plot.scatter(x='millan_v', y='thick', ax=ax1); ax1.set_title('millan_v');\n",
    "df.plot.scatter(x='itslive_v', y='thick', ax=ax2); ax2.set_title('itslive_v');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 2: aggregated per grid point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are so many points that much of the information obtained by OGGM is interpolated and therefore not biring much new information to a statistical model. A way to deal with this is to aggregate all the measurement points per grid point and to average them. Let's do this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_agg = df[['lon', 'lat', 'thick']].copy()\n",
    "ii, jj = gdir.grid.transform(df['lon'], df['lat'], crs=salem.wgs84, nearest=True)\n",
    "df_agg['i'] = ii\n",
    "df_agg['j'] = jj\n",
    "# We trick by creating an index of similar i's and j's\n",
    "df_agg['ij'] = ['{:04d}_{:04d}'.format(i, j) for i, j in zip(ii, jj)]\n",
    "df_agg = df_agg.groupby('ij').mean()\n",
    "# Conversion does not preserve ints\n",
    "df_agg['i'] = df_agg['i'].astype(int)\n",
    "df_agg['j'] = df_agg['j'].astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(df_agg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Select\n",
    "for vn in vns:\n",
    "    df_agg[vn] = ds[vn].isel(x=('z', df_agg['i']), y=('z', df_agg['j']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have 9 times less points, but the main features of the data remain unchanged:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))\n",
    "df_agg.plot.scatter(x='dis_from_border', y='thick', ax=ax1);\n",
    "df_agg.plot.scatter(x='slope', y='thick', ax=ax2);\n",
    "df_agg.plot.scatter(x='oggm_mb_above_z', y='thick', ax=ax3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a statistical model "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use scikit-learn to build a very simple linear model of ice-thickness. First, we have to acknowledge that there is a correlation between many of the explanatory variables we will use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "plt.figure(figsize=(10, 8));\n",
    "sns.heatmap(df.corr(), cmap='RdBu');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a problem for linear regression models, which cannot deal well with correlated explanatory variables. We have to do a so-called \"feature selection\", i.e. keep only the variables which are independently important to explain the outcome.\n",
    "\n",
    "For the sake of simplicity, let's use the [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) method to help us out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")  # sklearn sends a lot of warnings\n",
    "# Scikit learn (can be installed e.g. via conda install -c anaconda scikit-learn)\n",
    "from sklearn.linear_model import LassoCV\n",
    "from sklearn.model_selection import KFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Prepare our data\n",
    "df = df.dropna()\n",
    "# Variable to model\n",
    "target = df['thick']\n",
    "# Predictors - remove x and y (redundant with lon lat)\n",
    "# Normalize it in order to be able to compare the factors\n",
    "data = df.drop(['thick', 'x', 'y'], axis=1).copy()\n",
    "data_mean = data.mean()\n",
    "data_std = data.std()\n",
    "data = (data - data_mean) / data_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use cross-validation to select the regularization parameter\n",
    "lasso_cv = LassoCV(cv=5, random_state=0)\n",
    "lasso_cv.fit(data.values, target.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a statistical model trained on the full dataset. Let's see how it compares to the calibration data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "odf = df.copy()\n",
    "odf['thick_predicted'] = lasso_cv.predict(data.values)\n",
    "f, ax = plt.subplots(figsize=(6, 6));\n",
    "odf.plot.scatter(x='thick', y='thick_predicted', ax=ax);\n",
    "plt.xlim([-25, 220]);\n",
    "plt.ylim([-25, 220]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not too bad. It is doing OK for low thicknesses, but high thickness are strongly underestimated. Which explanatory variables did the model choose as being the most important?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "predictors = pd.Series(lasso_cv.coef_, index=data.columns)\n",
    "predictors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is interesting. Lons and lats have a predictive power over thickness? Unfortunately, this is more a coincidence than a reality. If we look at the correlation of the error with the variables of importance, one sees that there is a large correlation between the error and the spatial variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "odf['error'] = np.abs(odf['thick_predicted'] - odf['thick'])\n",
    "odf.corr()['error'].loc[predictors.loc[predictors != 0].index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, the model is not very robust. Let's use cross-validation to check our model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = KFold(4, random_state=0, shuffle=True)\n",
    "for k, (train, test) in enumerate(k_fold.split(data.values, target.values)):\n",
    "    lasso_cv.fit(data.values[train], target.values[train])\n",
    "    print(\"[fold {0}] alpha: {1:.5f}, score (r^2): {2:.5f}\".\n",
    "          format(k, lasso_cv.alpha_, lasso_cv.score(data.values[test], target.values[test])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fact that the hyper-parameter alpha and the score change that much between iterations is a sign that the model isn't very robust."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Our model is not great, but... let's apply it "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to apply the model to our entre glacier, we have to get the explanatory variables from the gridded dataset again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add variables we are missing\n",
    "lon, lat = gdir.grid.ll_coordinates\n",
    "ds['lon'] = (('y', 'x'), lon)\n",
    "ds['lat'] = (('y', 'x'), lat)\n",
    "\n",
    "# Generate our dataset\n",
    "pred_data = pd.DataFrame()\n",
    "for vn in data.columns:\n",
    "    pred_data[vn] = ds[vn].data[ds.glacier_mask == 1]\n",
    "\n",
    "# Normalize using the same normalization constants\n",
    "pred_data = (pred_data - data_mean) / data_std\n",
    "\n",
    "# Apply the model\n",
    "pred_data['thick'] = lasso_cv.predict(pred_data.values)\n",
    "\n",
    "# For the sake of physics, clip negative thickness values...\n",
    "pred_data['thick'] = np.clip(pred_data['thick'], 0, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Back to 2d and in xarray\n",
    "var = ds[vn].data * np.NaN\n",
    "var[ds.glacier_mask == 1] = pred_data['thick']\n",
    "ds['linear_model_thick'] = (('y', 'x'), var)\n",
    "ds['linear_model_thick'].attrs['description'] = 'Predicted thickness'\n",
    "ds['linear_model_thick'].attrs['units'] = 'm'\n",
    "ds['linear_model_thick'].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take home points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- OGGM provides preprocessed data from different sources on the same grid\n",
    "- we have shown how to compute gridded attributes from OGGM glaciers such as slope, aspect, or catchments\n",
    "- we used two methods to extract these data at point locations: with interpolation or with aggregated averages on each grid point\n",
    "- as an application example, we trained a linear regression model to predict the ice thickness of this glacier at unseen locations\n",
    "\n",
    "The model we developed was quite bad and we used quite lousy statistics. If I had more time to make it better, I would:\n",
    "- make a pre-selection of meaningful predictors to avoid discontinuities\n",
    "- use a non-linear model\n",
    "- use cross-validation to better asses the true skill of the model\n",
    "- ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## What's next?\n",
    "\n",
    "- return to the [OGGM documentation](https://docs.oggm.org)\n",
    "- back to the [table of contents](../welcome.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus: how does the statistical model compare to OGGM \"out-of the box\" and other thickness products?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write our thinckness estimates back to disk\n",
    "ds.to_netcdf(gdir.get_filepath('gridded_data'))\n",
    "# Distribute OGGM thickness using default values only\n",
    "workflow.execute_entity_task(tasks.distribute_thickness_per_altitude, gdirs);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:\n",
    "    ds = ds.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_agg['oggm_thick'] = ds.distributed_thickness.isel(x=('z', df_agg['i']), y=('z', df_agg['j']))\n",
    "df_agg['linear_model_thick'] = ds.linear_model_thick.isel(x=('z', df_agg['i']), y=('z', df_agg['j']))\n",
    "df_agg['millan_thick'] = ds.millan_ice_thickness.isel(x=('z', df_agg['i']), y=('z', df_agg['j']))\n",
    "df_agg['consensus_thick'] = ds.consensus_ice_thickness.isel(x=('z', df_agg['i']), y=('z', df_agg['j']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10));\n",
    "ds['linear_model_thick'].plot(ax=ax1); ax1.set_title('Statistical model');\n",
    "ds['distributed_thickness'].plot(ax=ax2); ax2.set_title('OGGM');\n",
    "ds['millan_ice_thickness'].where(ds.glacier_mask).plot(ax=ax3); ax3.set_title('Millan 2022');\n",
    "ds['consensus_ice_thickness'].plot(ax=ax4); ax4.set_title('Farinotti 2019');\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10));\n",
    "df_agg.plot.scatter(x='thick', y='linear_model_thick', ax=ax1);\n",
    "ax1.set_xlim([-25, 220]); ax1.set_ylim([-25, 220]); ax1.set_title('Statistical model');\n",
    "df_agg.plot.scatter(x='thick', y='oggm_thick', ax=ax2);\n",
    "ax2.set_xlim([-25, 220]); ax2.set_ylim([-25, 220]); ax2.set_title('OGGM');\n",
    "df_agg.plot.scatter(x='thick', y='millan_thick', ax=ax3);\n",
    "ax3.set_xlim([-25, 220]); ax3.set_ylim([-25, 220]); ax3.set_title('Millan 2022');\n",
    "df_agg.plot.scatter(x='thick', y='consensus_thick', ax=ax4);\n",
    "ax4.set_xlim([-25, 220]); ax4.set_ylim([-25, 220]); ax4.set_title('Farinotti 2019');"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}