Project 3 - Satellite imagery with rasterio et al.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 3 - Satellite imagery with rasterio _et al._\n",
    "\n",
    "![GDAL logo](img/gdal-logo.png)\n",
    "\n",
    "**rasterio** is a Python library to read and write geospatial gridded raster datasets such as satellite imagery and terrain models. It wraps the raster capabilities of GDAL, a very mature and well-established software made by the  Open Source Geospatial Foundation. It is often used in combination with **fiona** (for vector data, wrapping GDAL/OGR), **Shapely** (for planar geometries, wrapping GEOS) and **pyproj** (for cartographic projections and coordinate transformations, wrapping PROJ).\n",
    "\n",
    "- **Documentation**:\n",
    "    - https://rasterio.readthedocs.io/\n",
    "    - https://fiona.readthedocs.io/\n",
    "    - https://shapely.readthedocs.io/\n",
    "    - https://pyproj4.github.io/pyproj/\n",
    "- **Source code**:\n",
    "    - https://github.com/mapbox/rasterio/\n",
    "    - https://github.com/Toblerity/Fiona/\n",
    "    - https://github.com/Toblerity/Shapely\n",
    "    - https://github.com/pyproj4/pyproj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining raster images\n",
    "\n",
    "We will make use of **Sentinel 2** images over the city of Madrid and its surroundings.\n",
    "\n",
    "> Sentinel-2 is an Earth observation mission from the Copernicus Programme that systematically acquires optical imagery at high spatial resolution (10 m to 60 m) over land and coastal waters. The mission is a constellation with two twin satellites, Sentinel-2A and Sentinel-2B.\n",
    "\n",
    "> One of the benefits of the Copernicus Programme is that the data and information produced in the framework of Copernicus are made available free-of-charge to all its users and the public, thus allowing downstream services to be developed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sentinel data can be obtained from Python using <a href=\"https://pypi.org/project/sentinelsat/\"><code>sentinelsat</code></a>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from IPython.display import GeoJSON\n",
    "\n",
    "from sentinelsat import read_geojson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we will look for images covering downtown Madrid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roi_geojson = read_geojson(\"search_polygon.geojson\")\n",
    "GeoJSON(roi_geojson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you have a sentinel account, you can use this script\n",
    "# %run download_data.py\n",
    "# Alternatively, you can download it using the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This takes approximately 2 minutes, depending on Internet connection\n",
    "import urllib.request\n",
    "\n",
    "urllib.request.urlretrieve(\n",
    "    (\n",
    "        \"https://satellite-imagery-storage.s3.fr-par.scw.cloud/\"\n",
    "        \"S2B_MSIL2A_20210930T105749_N0301_R094_T30TVK_20210930T125620_MANUAL.zip\"\n",
    "    ),\n",
    "    \"data/S2B_MSIL2A_20210930T105749_N0301_R094_T30TVK_20210930T125620.zip\"\n",
    ")\n",
    "\n",
    "product_title = \"S2B_MSIL2A_20210930T105749_N0301_R094_T30TVK_20210930T125620\"\n",
    "product_filename = \"S2B_MSIL2A_20210930T105749_N0301_R094_T30TVK_20210930T125620.SAFE\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, we unzip the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import shutil\n",
    "\n",
    "product_file_path = Path(\"data\") / (product_title + \".zip\")\n",
    "\n",
    "shutil.unpack_archive(product_file_path, \"./data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Reading raster data\n",
    "\n",
    "The True Color Image (TCI) corresponding to our data looks more or less like this:\n",
    "\n",
    "![Sentinel Madrid](img/sentinel-madrid.png)\n",
    "\n",
    "However, **we will not load it entirely in memory**, because the data is much bigger than our available RAM. Instead, we will use a window to load the data partially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "\n",
    "true_color_path = list((Path(\"data\") / product_filename).glob(\"**/*_TCI_10m.jp2\"))[0]\n",
    "\n",
    "src = rasterio.open(true_color_path)\n",
    "src"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This returns an open dataset, acting like a pointer to the actual data on disk. Therefore, we are still not loading it.\n",
    "\n",
    "We can access to several metadata attributes of this dataset, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.crs  # Coordinate reference system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.bounds  # Expressed in the CRS above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.count  # Number of bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.shape  # Shape of the underlying arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.colorinterp  # Color interpretation of each band"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vector data\n",
    "\n",
    "We are interested in downtown Madrid, therefore we will create a rectangle enclosing the Region of Interest (ROI). Using Shapely, we can load our GeoJSON data as raw coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import shape\n",
    "\n",
    "roi = shape(roi_geojson)\n",
    "roi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we have to transform these coordinates into the CRS of the raster data. For that, we need another Python package called pyproj:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyproj import Transformer\n",
    "from shapely.ops import transform\n",
    "\n",
    "def reproject_vector(vector, to_epsg, from_epsg=4326):\n",
    "    # https://shapely.readthedocs.io/en/latest/manual.html#other-transformations\n",
    "    project = Transformer.from_crs(from_epsg, to_epsg, always_xy=True).transform\n",
    "    return transform(project, roi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = reproject_vector(roi, src.crs.to_epsg()).bounds\n",
    "bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use these bounds to create a rasterio `Window` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "win = src.window(*bounds)\n",
    "win"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, we can use this to load and visualize a partial subset of our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show\n",
    "\n",
    "show(src.read(window=win), transform=src.transform)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, we can load just one of the bands specifying an integer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "show(src.read(1, window=win), transform=src.transform, cmap=\"Reds\", title=\"Red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or display them all at the same time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, (ax_r, ax_g, ax_b) = plt.subplots(ncols=3, figsize=(21, 7))\n",
    "\n",
    "show(src.read(1, window=win), ax=ax_r, cmap='Reds', title='Red')\n",
    "show(src.read(2, window=win), ax=ax_g, cmap='Greens', title='Green')\n",
    "show(src.read(3, window=win), ax=ax_b, cmap='Blues', title='Blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading external bands\n",
    "\n",
    "According to the [Sentinel 2 official definitions](https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/definitions):\n",
    "\n",
    "> The TCI is an RGB image built from the B02 (Blue), B03 (Green), and B04 (Red) Bands.\n",
    "\n",
    "Therefore, we could load these bands separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b02 = rasterio.open(list((Path(\"data\") / product_filename).glob(\"**/*_B02_10m.jp2\"))[0])\n",
    "b03 = rasterio.open(list((Path(\"data\") / product_filename).glob(\"**/*_B03_10m.jp2\"))[0])\n",
    "b04 = rasterio.open(list((Path(\"data\") / product_filename).glob(\"**/*_B04_10m.jp2\"))[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the range of these files is much wider:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b04.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(b04.read(window=win), transform=src.transform, cmap='Reds', title='Red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The saturation level of 255 digital counts correspond to a level of 3558 for L1C products or 2000 for L2A products (0.3558 and 0.2 in reflectance value respectively.\n",
    "\n",
    "Therefore, we need to rescale the data manually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(\n",
    "    (b04.read(window=win).clip(max=3558) / 3558 * 255).clip(min=1).astype(np.uint8),\n",
    "    transform=src.transform,\n",
    "    cmap='Reds',\n",
    "    title='Red',\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}