{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CryoTEMPO-EOLIS Gridded product - elevation difference across 8 years"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example retrieves two sets of [CryoTEMPO-EOLIS Gridded Product](https://cryotempo-eolis.org/gridded-product/) data 8 years apart\n",
    "and plots the elevation difference over Iceland.\n",
    "\n",
    "You can view the code below or [click here](https://colab.research.google.com/github/earthwave/specklia_demo_notebooks/blob/main/cryotempo_gridded_product.ipynb/) to run it yourself in Google Colab!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Environment Setup\n",
    "To run this notebook, you will need to make sure that the folllowing packages are installed in your python environment (all can be installed via pip/conda):\n",
    "- matplotlib\n",
    "- geopandas\n",
    "- contextily\n",
    "- ipykernel\n",
    "- shapely\n",
    "- specklia\n",
    "\n",
    "If you are using the Google Colab environment, these packages will be installed in the next cell. Please note this step may take a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    %pip install rasterio --no-binary rasterio\n",
    "    %pip install specklia\n",
    "    %pip install matplotlib\n",
    "    %pip install geopandas\n",
    "    %pip install contextily\n",
    "    %pip install shapely\n",
    "    %pip install python-dotenv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "# fix an issue that can sometimes occur with rasterio using the wrong version of proj\n",
    "import os\n",
    "import pyproj\n",
    "os.environ['PROJ_LIB'] = pyproj.datadir.get_data_dir()\n",
    "\n",
    "from datetime import datetime\n",
    "import math\n",
    "\n",
    "import contextily as ctx\n",
    "from dotenv import load_dotenv\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import shapely\n",
    "\n",
    "from specklia import Specklia\n",
    "\n",
    "# load a demonstration API key from a .env file.\n",
    "load_dotenv()\n",
    "\n",
    "# To run this code yourself, first generate your own key using https://specklia.earthwave.co.uk.\n",
    "if 'DEMO_API_KEY' in os.environ:\n",
    "    client = Specklia(os.environ['DEMO_API_KEY'])\n",
    "else:\n",
    "    user_api_key = input('Please generate your own key using https://specklia.earthwave.co.uk/ApiKeys and paste it here:')\n",
    "    client = Specklia(user_api_key)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Query dataset"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's use Specklia to retrieve Gridded Product data between January and March of 2012."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iceland_polygon = shapely.Polygon(([-26., 63.], [-26., 68.], [-11., 68.], [-11., 63.], [-26., 63.]))\n",
    "\n",
    "dataset_name = 'CryoTEMPO-EOLIS Gridded Product'\n",
    "available_datasets = client.list_datasets()\n",
    "gridded_product_dataset = available_datasets[\n",
    "    available_datasets['dataset_name'] == dataset_name].iloc[0]\n",
    "\n",
    "gridded_product_data_2012, sources = client.query_dataset(\n",
    "    dataset_id=gridded_product_dataset['dataset_id'],\n",
    "    epsg4326_polygon=iceland_polygon,\n",
    "    min_datetime=datetime(2012, 1, 1),\n",
    "    max_datetime=datetime(2012, 3, 31, 23, 59, 59),\n",
    "    columns_to_return=['timestamp', 'elevation'])\n",
    "\n",
    "print(\n",
    "    f'Query complete, {len(gridded_product_data_2012)} points returned, drawn from {len(sources)} original sources.')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's retrieve and plot the same months for 2020."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridded_product_data_2020, sources = client.query_dataset(\n",
    "    dataset_id=gridded_product_dataset['dataset_id'],\n",
    "    epsg4326_polygon=iceland_polygon,\n",
    "    min_datetime=datetime(2020, 1, 1),\n",
    "    max_datetime=datetime(2020, 3, 31, 23, 59, 59),\n",
    "    columns_to_return=['timestamp', 'elevation'])\n",
    "\n",
    "print(\n",
    "    f'Query complete, {len(gridded_product_data_2020)} points returned, drawn from {len(sources)} original sources.')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the Gridded Product  [documentation](https://cryotempo-eolis.org/gridded-product/), we know that it is generated at\n",
    "a monthly temporal resolution, one product file per region.\n",
    "\n",
    "As we have queried a single region over the span of three\n",
    "months, we received data from three source files - one for each month."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot DEMs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's quickly define a plotting function to investigate our results. The function plots the DEM of each Gridded Product\n",
    "source file in our Specklia query.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_gridded_product_per_source_id(\n",
    "        gridded_product_data: gpd.GeoDataFrame, n_fig_cols: int = 3, fig_epsg: int = 3413):\n",
    "    gridded_product_data_grouped_by_source_id = gridded_product_data.groupby(\n",
    "        ['source_id'], sort=False)\n",
    "\n",
    "    n_fig_rows = math.ceil(\n",
    "        gridded_product_data_grouped_by_source_id.ngroups / n_fig_cols)\n",
    "    _, axs = plt.subplots(n_fig_rows, n_fig_cols,\n",
    "                          figsize=(n_fig_cols*10, n_fig_rows*10))\n",
    "\n",
    "    for i, grp in enumerate(gridded_product_data_grouped_by_source_id):\n",
    "        _, data_for_source_id = grp\n",
    "        ax = axs.flatten()[i]\n",
    "        centre_date_of_product = datetime.fromtimestamp(\n",
    "            data_for_source_id['timestamp'].aggregate('mean'))\n",
    "        centre_date_of_product.strftime('%m-%Y')\n",
    "        data_for_source_id.to_crs(epsg=fig_epsg).plot(\n",
    "            ax=ax, column=data_for_source_id['elevation'], s=2, cmap='viridis', legend=True, legend_kwds={\n",
    "                \"label\": \"Elevation [m]\", 'orientation': 'horizontal', 'shrink': .9, 'pad': .1\n",
    "            })\n",
    "        ax.set_title(\n",
    "            f'DEM for CryoTEMPO-EOLIS Gridded Product, {centre_date_of_product.date()}')\n",
    "\n",
    "        ax.set_xlabel('x [m]')\n",
    "        ax.set_ylabel('y [m]')\n",
    "        ctx.add_basemap(\n",
    "            ax, source=ctx.providers.Esri.WorldImagery, crs=fig_epsg, zoom=6)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's use our function to plot the results of our Specklia queries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridded_product_data_2012 = gridded_product_data_2012.sort_values(['timestamp'])\n",
    "\n",
    "plot_gridded_product_per_source_id(gridded_product_data_2012)\n",
    "\n",
    "gridded_product_data_2020 = gridded_product_data_2020.sort_values(['timestamp'])\n",
    "\n",
    "plot_gridded_product_per_source_id(gridded_product_data_2020)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, in the Iceland region, the Gridded Product has particularly good coverage over Vatnajökull ice cap."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot elevation difference"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the elevation difference between January 2012 and January 2020 data to investigate further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jan_2012_source_id = gridded_product_data_2012.iloc[0]['source_id']\n",
    "jan_2020_source_id = gridded_product_data_2020.iloc[0]['source_id']\n",
    "\n",
    "jan_2012_data = gridded_product_data_2012[gridded_product_data_2012['source_id']\n",
    "                                          == jan_2012_source_id]\n",
    "jan_2020_data = gridded_product_data_2020[gridded_product_data_2020['source_id']\n",
    "                                          == jan_2020_source_id]\n",
    "\n",
    "merged_gdf = jan_2012_data.sjoin(df=jan_2020_data, how='inner')\n",
    "merged_gdf['elevation_difference'] = merged_gdf['elevation_right'] - \\\n",
    "    merged_gdf['elevation_left']\n",
    "\n",
    "ax = merged_gdf.to_crs(epsg=3413).plot(column=merged_gdf['elevation_difference'], figsize=(\n",
    "    10, 10), legend=True, cmap='magma_r', s=3, legend_kwds={\n",
    "        'label': 'Elevation change [m]', 'orientation': 'horizontal'\n",
    "    })\n",
    "ax.set_xlabel('x [m]')\n",
    "ax.set_ylabel('y [m]')\n",
    "ax.set_title(f\"Elevation difference between Jan 2012 and Jan 2020 in Iceland\")\n",
    "ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=3413, zoom=8)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above shows a significant decrease of ice elevation in the glacier margins at low elevations, this is particularly visible in Vatnajökull ice cap."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "specklia",
   "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.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}