{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ab31574b",
   "metadata": {},
   "source": [
    "# STAC Metadata Visualizations\n",
    "\n",
    "This notebook illustrates a simple way to display footprints of discovered Items after searching a STAC API, and making simple plots using Pandas and Holoviews. Only the metadata is visualized in these examples through maps and plots. The actual STAC data (i.e., Item Assets) are not accessed.\n",
    "\n",
    "The libraries GeoPandas and hvplot are used for visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3188b5a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystac_client import Client\n",
    "\n",
    "# set pystac_client logger to DEBUG to see API calls\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logger = logging.getLogger('pystac_client')\n",
    "logger.setLevel(logging.DEBUG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d404982",
   "metadata": {},
   "source": [
    "Define the STAC API to use, along with any custom headers (such as for authentication)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d40ed5c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# STAC API root URL\n",
    "URL = 'https://planetarycomputer.microsoft.com/api/stac/v1'\n",
    "\n",
    "# custom headers\n",
    "headers = []\n",
    "\n",
    "cat = Client.open(URL, headers=headers)\n",
    "cat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e513e548",
   "metadata": {},
   "source": [
    "## Search\n",
    "\n",
    "Perform a spatio-temporal search of ASTER data for a small AOI in the northern part of The Netherlands between 2000 and 2010."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae7a3eca",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# AOI around Delfzijl, in northern Netherlands\n",
    "geom = {\n",
    "    \"type\": \"Polygon\",\n",
    "    \"coordinates\": [\n",
    "      [\n",
    "        [\n",
    "          6.42425537109375,\n",
    "          53.174765470134616\n",
    "        ],\n",
    "        [\n",
    "          7.344360351562499,\n",
    "          53.174765470134616\n",
    "        ],\n",
    "        [\n",
    "          7.344360351562499,\n",
    "          53.67393435835391\n",
    "        ],\n",
    "        [\n",
    "          6.42425537109375,\n",
    "          53.67393435835391\n",
    "        ],\n",
    "        [\n",
    "          6.42425537109375,\n",
    "          53.174765470134616\n",
    "        ]\n",
    "      ]\n",
    "    ]\n",
    "}\n",
    "\n",
    "# limit sets the # of items per page so we can see multiple pages getting fetched\n",
    "search = cat.search(\n",
    "    max_items = 50,\n",
    "    collections = \"aster-l1t\",\n",
    "    intersects = geom,\n",
    "    datetime = \"2000-01-01/2010-12-31\",\n",
    ")\n",
    "\n",
    "# retrieve the items as dictionaries, rather than Item objects\n",
    "items = list(search.items_as_dicts())\n",
    "len(items)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a23660f",
   "metadata": {},
   "source": [
    "## GeoPandas\n",
    "\n",
    "A GeoDataFrame is constructed from the AOI geometry, in order to make visualizations easy.\n",
    "\n",
    "The STAC Items, which are a GeoJSON FeatureCollection can be converted to a GeoDataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed0da77b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from shapely.geometry import shape\n",
    "\n",
    "# convert a list of STAC Items into a GeoDataFrame\n",
    "def items_to_geodataframe(items):\n",
    "    _items = []\n",
    "    for i in items:\n",
    "        _i = deepcopy(i)\n",
    "        _i['geometry'] = shape(_i['geometry'])\n",
    "        _items.append(_i)\n",
    "    gdf = gpd.GeoDataFrame(pd.json_normalize(_items))\n",
    "    for field in ['properties.datetime', 'properties.created', 'properties.updated']:\n",
    "        if field in gdf:\n",
    "            gdf[field] = pd.to_datetime(gdf[field])\n",
    "    gdf.set_index('properties.datetime', inplace=True)\n",
    "    return gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4bed793",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert geometry to a GeoDataFrame\n",
    "aoi_gdf = gpd.GeoDataFrame([{'geometry': shape(geom)}])\n",
    "aoi_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1fc8732",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert found items to a GeoDataFrame\n",
    "items_gdf = items_to_geodataframe(items)\n",
    "items_gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d570626",
   "metadata": {},
   "source": [
    "## Plot Geometries on a Map\n",
    "\n",
    "Holoviews is used to display geometries on a map by using `hvplot`. The The `*` Holoviews operator to overlay two plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fbfad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.pandas\n",
    "\n",
    "# plot polygons on a map with background tiles.\n",
    "def plot_polygons(data, *args, **kwargs):\n",
    "    return data.hvplot.polygons(*args, geo=True, projection=\"GOOGLE_MERCATOR\", xaxis=None, yaxis=None,\n",
    "                             frame_width=600, frame_height=600, fill_alpha=0,\n",
    "                             line_width=4, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72094092",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polygons(aoi_gdf, tiles=\"OSM\", line_color=\"red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7d6cfc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polygons(items_gdf, tiles=\"OSM\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a07e87b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polygons(items_gdf, tiles=\"OSM\") * plot_polygons(aoi_gdf, line_color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0604c84",
   "metadata": {},
   "source": [
    "## Line Plots\n",
    "\n",
    "Numeric STAC metadata can also be plotted, most often this will be plotted vs the Item `datetime`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c198e261",
   "metadata": {},
   "outputs": [],
   "source": [
    "items_df = pd.DataFrame(items_gdf)\n",
    "\n",
    "plot_fields = [\n",
    "    'properties.eo:cloud_cover',\n",
    "    'properties.view:sun_azimuth',\n",
    "    'properties.view:sun_elevation'\n",
    "]\n",
    "\n",
    "items_df[plot_fields].hvplot(height=500, width=800).opts(legend_position=\"top_right\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "6b6313dbab648ff537330b996f33bf845c0da10ea77ae70864d6ca8e2699c7ea"
  },
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}