{
 "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",
    "\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",
    "\n",
    "import pyproj\n",
    "\n",
    "os.environ[\"PROJ_LIB\"] = pyproj.datadir.get_data_dir()\n",
    "\n",
    "import math\n",
    "from datetime import datetime\n",
    "\n",
    "import contextily as ctx\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import shapely\n",
    "from dotenv import load_dotenv\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(\n",
    "        \"Please generate your own key using https://specklia.earthwave.co.uk/ApiKeys and paste it here:\"\n",
    "    )\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.0, 63.0], [-26.0, 68.0], [-11.0, 68.0], [-11.0, 63.0], [-26.0, 63.0]))\n",
    "\n",
    "dataset_name = \"CryoTEMPO-EOLIS Gridded Product\"\n",
    "available_datasets = client.list_datasets()\n",
    "gridded_product_dataset = available_datasets[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",
    "\n",
    "print(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",
    "\n",
    "print(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",
    ") -> None:\n",
    "    gridded_product_data_grouped_by_source_id = gridded_product_data.groupby([\"source_id\"], sort=False)\n",
    "\n",
    "    n_fig_rows = math.ceil(gridded_product_data_grouped_by_source_id.ngroups / n_fig_cols)\n",
    "    _, axs = plt.subplots(n_fig_rows, n_fig_cols, 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(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,\n",
    "            column=data_for_source_id[\"elevation\"],\n",
    "            s=2,\n",
    "            cmap=\"viridis\",\n",
    "            legend=True,\n",
    "            legend_kwds={\"label\": \"Elevation [m]\", \"orientation\": \"horizontal\", \"shrink\": 0.9, \"pad\": 0.1},\n",
    "        )\n",
    "        ax.set_title(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(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\"] == jan_2012_source_id]\n",
    "jan_2020_data = gridded_product_data_2020[gridded_product_data_2020[\"source_id\"] == 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\"] - merged_gdf[\"elevation_left\"]\n",
    "\n",
    "ax = merged_gdf.to_crs(epsg=3413).plot(\n",
    "    column=merged_gdf[\"elevation_difference\"],\n",
    "    figsize=(10, 10),\n",
    "    legend=True,\n",
    "    cmap=\"magma_r\",\n",
    "    s=3,\n",
    "    legend_kwds={\"label\": \"Elevation change [m]\", \"orientation\": \"horizontal\"},\n",
    ")\n",
    "ax.set_xlabel(\"x [m]\")\n",
    "ax.set_ylabel(\"y [m]\")\n",
    "ax.set_title(\"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
}