{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CryoTEMPO-EOLIS products over Austfonna Basin 3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "This example shows how to use Specklia to retrieve [CryoTEMPO-EOLIS Point](https://cryotempo-eolis.org/point-product/) and [Gridded Product](https://cryotempo-eolis.org/gridded-product/) data over [Austfonna Basin 3](https://www.swisseduc.ch/glaciers/svalbard/austfonna/index-en.html).\n",
    "\n",
    "You can view the code below or [click here](https://colab.research.google.com/github/earthwave/specklia_demo_notebooks/blob/main/austfonna_basin_3.ipynb/) to run it yourself in Google Colab!"
   ]
  },
  {
   "attachments": {},
   "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": "markdown",
   "metadata": {},
   "source": [
    "Now, let's import the packages and create a Specklia client to work with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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 pprint\n",
    "from time import perf_counter\n",
    "\n",
    "import contextily as ctx\n",
    "from dotenv import load_dotenv\n",
    "import geopandas as gpd\n",
    "from IPython.display import display\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from pandas import Timedelta\n",
    "import shapely\n",
    "from specklia import Specklia\n",
    "\n",
    "# load a demonstration API key from a .env file.\n",
    "load_dotenv()\n",
    "\n",
    "# Specify polygon colors to match another Earthwave product, https://cs2eo.org.\n",
    "data_coverage_color = np.array([52, 211, 153])/255\n",
    "geo_filter_colour = np.array([96, 165, 250])/255\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": [
    "### Listing groups and datasets\n",
    "The demonstration user has the same permissions as a freshly created Specklia user.\n",
    "Let's list the groups they are a part of and the datasets they can access:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(client.list_groups())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "available_datasets = client.list_datasets()\n",
    "display(available_datasets)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting dataset spatial coverage\n",
    "There's a lot of information here. Let's plot the geospatial coverage of the Point and Gridded product datasets, using [Contextily](https://contextily.readthedocs.io/en/latest/) to add a satellite background map. Note that as of 30th June 2023, only products covering Svalbard have been loaded. We will be loading more in the second half of 2023.\n",
    "\n",
    "All background satellite imagery tiles are provided by the [Esri World Imagery Service](https://www.esriuk.com/en-gb/content/products?esri-world-imagery-service) via Contextily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eolis_datasets = {}\n",
    "for ds_name in ['CryoTEMPO-EOLIS Point Product', 'CryoTEMPO-EOLIS Gridded Product']:\n",
    "    eolis_datasets[ds_name] = available_datasets[\n",
    "        available_datasets['dataset_name'] == ds_name].iloc[0]\n",
    "\n",
    "    print(f\"{ds_name} contains data timestamped between \\n\"\n",
    "        f\"{eolis_datasets[ds_name]['min_timestamp']} \"\n",
    "        f\"and {eolis_datasets[ds_name]['max_timestamp']}\\n\")\n",
    "\n",
    "    print(f\"{ds_name} has the following columns:\")\n",
    "    pprint.PrettyPrinter(indent=2, width=120).pprint(eolis_datasets[ds_name]['columns'])\n",
    "    print(\"\\n\\n\")\n",
    "\n",
    "    desired_dataset_spatial_coverage = gpd.GeoDataFrame(\n",
    "        geometry=[eolis_datasets[ds_name]['epsg4326_coverage']], crs=4326)\n",
    "    \n",
    "    # we create two plots, one for each hemisphere, to minimise distortion\n",
    "    # and illustrate that multiple regions may be present within the same dataset\n",
    "    for crs, hemisphere_name, hemisphere_bounding_box in [\n",
    "            (3413, 'northern', (-180, 0, 180, 90)),\n",
    "            (3031, 'southern', (-180, -90, 180, 0))\n",
    "            ]:\n",
    "\n",
    "        cropped_data = desired_dataset_spatial_coverage.clip(hemisphere_bounding_box)\n",
    "\n",
    "        if len(cropped_data) > 0:\n",
    "            transformed_cropped_data = cropped_data.to_crs(crs)\n",
    "            ax = transformed_cropped_data.plot(\n",
    "                figsize=(10, 10), alpha=0.5, color=data_coverage_color, edgecolor=data_coverage_color)\n",
    "\n",
    "            ax.set_xlabel('x (m)')\n",
    "            ax.set_ylabel('y (m)')\n",
    "            ax.set_title(f\"{ds_name} spatial coverage, {hemisphere_name} hemisphere (EPSG {crs})\")\n",
    "            # we calculate the zoom level dynamically so that we don't have to change this code\n",
    "            # when new regions are added to the CryoTEMPO EOLIS products.\n",
    "            ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=crs, attribution=False,\n",
    "                            zoom=ctx.tile._calculate_zoom(*cropped_data.total_bounds) - 1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Querying RGI polygons\n",
    "\n",
    "In this example, we are interested in studying Austfonna Basin 3. Conveniently, Specklia contains all [RGI v7.0](https://www.glims.org/RGI/) glacier boundary definitions. Let's retrieve the RGI v7.0 boundary of Basin 3 from Specklia and take a look.\n",
    "\n",
    "Firstly, we need to define a rough polygon covering the area to query the RGIv7 polygons from Specklia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "austfonna_extent = shapely.from_wkt(\n",
    "    'POLYGON ((19.27002 79.492646, 23.862305 79.075977, 28.630371 79.843346, 25.708008 80.441282, 19.27002 79.492646))')\n",
    "austfonna_extent_gdf = gpd.GeoSeries(austfonna_extent, crs=4326)\n",
    "# We plot in the WGS 84 / NSIDC Sea Ice Polar Stereographic North CRS (EPSG 3413).\n",
    "ax = austfonna_extent_gdf.to_crs(3413).plot(\n",
    "    figsize=(10, 10), alpha=0.5, color=geo_filter_colour, edgecolor=geo_filter_colour)\n",
    "\n",
    "ax.set_xlabel('x (meters)')\n",
    "ax.set_ylabel('y (meters)')\n",
    "ax.set_title(\"Query extent for Austfonna\")\n",
    "ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=3413, attribution=False, zoom=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have defined a rough extent covering Austfonna, we can query the polygons from Specklia. Specklia will return all RGI Glaciers whose centroids fall within the query polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgiv7_dataset = available_datasets[available_datasets['dataset_name'] == 'RGIv7 Glaciers'].iloc[0]\n",
    "\n",
    "query_start_time = perf_counter()\n",
    "\n",
    "rgi_data, rgi_data_sources = client.query_dataset(\n",
    "    dataset_id=rgiv7_dataset['dataset_id'],\n",
    "    epsg4326_polygon=austfonna_extent,\n",
    "    min_datetime=rgiv7_dataset['min_timestamp'] - Timedelta(seconds=1),\n",
    "    max_datetime=rgiv7_dataset['max_timestamp'] + Timedelta(seconds=1))\n",
    "\n",
    "print(f'Query took {perf_counter()-query_start_time:.2f} seconds to complete.')\n",
    "print(f'Austfonna Polygon Query complete, {len(rgi_data)} points returned, '\n",
    "        f'drawn from {len(rgi_data_sources)} original sources.')\n",
    "print(f'Columns within the data: {rgi_data.columns}\\n\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Specklia query returns all of the original attributes in the RGI shapefile, as well as varying degrees of simplified versions of the original polygon. The simplified polygons have fewer vertices which make for faster querying, if we retrieve only those.\n",
    "\n",
    "Additionally, source information is supplied for the returned data, and a direct link to the source file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(rgi_data_sources)\n",
    "display(rgi_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the RGI v7.0 boundary definition for Austonna ice cap's Basin 3, by extracting it from the returned polygons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "austfonna_basin_3 = gpd.GeoSeries(\n",
    "    rgi_data.loc[rgi_data['rgi_id'] == 'RGI2000-v7.0-G-07-01383', 'original_polygons'], crs=4326)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "austfonna_basin_3.to_crs(3413).plot(ax=ax, facecolor=geo_filter_colour, edgecolor=geo_filter_colour, alpha=0.5)\n",
    "\n",
    "ax.set_xlim(982000, 1106000)\n",
    "ax.set_ylim(-484000, -322000)\n",
    "\n",
    "ax.set_xlabel('x (meters)')\n",
    "ax.set_ylabel('y (meters)')\n",
    "ax.set_title(f\"Austfonna Basin 3 ({rgi_data['rgi_id'].values[0]}, EPSG 3413)\")\n",
    "ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=3413, attribution=False, zoom=8)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Querying data\n",
    "We can now query the two Specklia datasets for our study region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cryosat_data = {}\n",
    "sources = {}\n",
    "additional_filters = {\n",
    "    'CryoTEMPO-EOLIS Point Product': [\n",
    "        {'column': 'uncertainty', 'operator': '<=', 'threshold': 4}],\n",
    "    'CryoTEMPO-EOLIS Gridded Product': []}\n",
    "\n",
    "for ds_name, data in eolis_datasets.items():\n",
    "    query_start_time = perf_counter()\n",
    "    cryosat_data[ds_name], sources[ds_name] = client.query_dataset(\n",
    "        dataset_id=data['dataset_id'],\n",
    "        epsg4326_polygon=austfonna_basin_3.iloc[0],\n",
    "        min_datetime=datetime(2015, 1, 1),\n",
    "        max_datetime=datetime(2015, 12, 1),\n",
    "        columns_to_return=['timestamp', 'elevation', 'uncertainty'],\n",
    "        additional_filters=additional_filters[ds_name])\n",
    "\n",
    "    print(f'Query took {perf_counter()-query_start_time:.2f} seconds to complete.')\n",
    "    print(f'{ds_name} Query complete, {len(cryosat_data[ds_name])} points returned, '\n",
    "          f'drawn from {len(sources[ds_name])} original sources.')\n",
    "    print(f'Columns within the data: {cryosat_data[ds_name].columns}\\n\\n')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that within `sources` we have all of the source information from the original product files, including both the Earth Explorer and the NetCDF4 headers, and with the `source_id` and `source_row_id` columns, the ability to trace each point we have just retrieved from its original source file.\n",
    "\n",
    "If we wanted to, we could run a query that would give us back 100% of the information available in the originally ingested file.\n",
    "\n",
    "Lastly, we plot the CryoSat-2 data over the top of our study polygon, illustrating the filtering that has been performed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ds_name, cs_data in cryosat_data.items():\n",
    "    # create a normalised colormap\n",
    "    norm = mcolors.Normalize(vmin=cs_data['elevation'].min(), vmax=cs_data['elevation'].max())\n",
    "\n",
    "    ax = austfonna_basin_3.to_crs(epsg=3413).plot(figsize=(10, 10), alpha=0.5)\n",
    "    cs_data.to_crs(epsg=3413).plot(ax=ax, column='elevation', markersize=.1, norm=norm)\n",
    "    \n",
    "    ax.set_xlabel('EPSG 3413 x (m)')\n",
    "    ax.set_ylabel('EPSG 3413 y (m)')\n",
    "    ax.set_title(f\"{ds_name} Elevation over Austfonna Basin 3 \\n(RGI v6.0, EPSG 3413)\")\n",
    "    ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=3413, attribution=False, zoom=10)\n",
    "    cbar = plt.gcf().colorbar(plt.cm.ScalarMappable(norm=norm), ax=ax, shrink=.5)\n",
    "    cbar.set_label('Elevation (m)')"
   ]
  }
 ],
 "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
}