01_Find_Load_Visualise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find, Load, and Visualise Earth Observation Imagery\n",
    "\n",
    "This notebook demonstrates how to find, load, and visualise Earth observation imagery using cloud native approaches, which work well on your desktop or in cloud environments!\n",
    "\n",
    "The notebook demonstrates how to find images from the Element 84 Earth Search STAC catalog, load them using `odc-stac`, and visualise them using extensions to xarray provided by `odc-geo`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up\n",
    "\n",
    "The first step is to set up the required Python libraries and local imports.\n",
    "\n",
    "* `odc.stac` and `pystac_client` are used to access the Earth Search STAC catalog\n",
    "* `numpy` is used to manipulate data\n",
    "* `odc.geo.xr` enables additional functionality for xarrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from odc.stac import configure_s3_access, load\n",
    "from pystac.client import Client\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second step is to start a Dask client.\n",
    "\n",
    "Dask supports local parallel processing and can help speed up computation times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client as DaskClient\n",
    "\n",
    "dask_client = DaskClient()\n",
    "dask_client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Find\n",
    "\n",
    "### 1.1 Connecting to the catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The catalog URL for the Earth Search STAC API\n",
    "catalog = \"https://earth-search.aws.element84.com/v1\"\n",
    "\n",
    "# Pystac_client's Client class is used to connect to the catalog\n",
    "stac_client = Client.open(catalog)\n",
    "\n",
    "# Configure settings for reading from Earth Search STAC\n",
    "configure_s3_access(\n",
    "    cloud_defaults=True,\n",
    "    aws_unsigned=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Selecting an area to query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from odc.geo import BoundingBox\n",
    "\n",
    "bbox = BoundingBox(\n",
    "    left=-113.2,\n",
    "    bottom=37.0,\n",
    "    right=-112.8,\n",
    "    top=37.3,\n",
    "    crs=\"EPSG:4326\",\n",
    ")\n",
    "bbox.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Set year and month to query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set a start and end date\n",
    "date_query = \"2025-04\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Choose collections and filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set Earth Search product ID as the STAC \"Collection\"\n",
    "collections_query = [\"sentinel-2-c1-l2a\"]\n",
    "\n",
    "# Set up a filter query. This is less than 20% cloud cover\n",
    "filter_query = {\"eo:cloud_cover\": {\"lt\": 20}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Running the query to indentify matching STAC items"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Query with filtering for cloud cover\n",
    "items = stac_client.search(\n",
    "    intersects=bbox.boundary(),\n",
    "    collections=collections_query,\n",
    "    datetime=date_query,\n",
    "    query=filter_query,\n",
    ").item_collection()\n",
    "\n",
    "print(f\"Found {len(items)} items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load\n",
    "\n",
    "### 2.1 Using odc-stac to load identified items\n",
    "\n",
    "This may take a few minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load our filtered data\n",
    "ds_filtered = load(\n",
    "    items,\n",
    "    bands=[\"red\", \"green\", \"blue\"],\n",
    "    crs=\"utm\",\n",
    "    chunks={},\n",
    "    resolution=20,\n",
    "    groupby=\"solar_day\",\n",
    "    intersects=bbox.boundary(),\n",
    ").compute()\n",
    "\n",
    "ds_filtered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Review loaded imagery\n",
    "\n",
    "Identify which image you want to export and note the date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To_array sets up a 3D array with the time dimension, which works directly\n",
    "# with the plot function to make an RGB image\n",
    "ds_filtered.to_array().plot.imshow(col=\"time\", col_wrap=3, robust=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Visualise\n",
    "\n",
    "### 3.1 Select best image\n",
    "\n",
    "Update the `best_image_date` parameter to match the date you identified in the previous step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_image_date = \"2025-04-06\"\n",
    "\n",
    "best_image = ds_filtered.sel(time=best_image_date).squeeze()\n",
    "\n",
    "best_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 View the selected image on an interactive map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualisation = best_image.odc.to_rgba()\n",
    "\n",
    "visualisation.odc.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Improve the image contrast\n",
    "\n",
    "Calculate the values corresponding to the 1st and 99th percentiles. \n",
    "These can be used in the `to_rgba()` function to stretch the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percentile_stretch = (1, 99)\n",
    "\n",
    "rgb_array = best_image.to_array().values\n",
    "\n",
    "stretch_vmin, stretch_vmax = np.nanpercentile(rgb_array, percentile_stretch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the percentile stretch values to the visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch_visualisation = best_image.odc.to_rgba(vmin=stretch_vmin, vmax=stretch_vmax)\n",
    "\n",
    "stretch_visualisation.odc.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Export to a cloud-optimised GeoTIFF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch_visualisation.odc.write_cog(\"sentinel2_example.tif\", overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}