{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Sentinel-2 processing to get a burnt area map (Sen2Cor, SNAP, snappy)\n", "\n", "At the end of May 2018 a wildfire began in the North of New Mexico next to Cimarron in Ute Park area. The fire burnt a minimum of 36,000 acres (14,569 hectares) of drought-parched grassland and timber (source: [Express](http://web.archive.org/web/20180605150346/https://www.express.co.uk/news/world/969726/ute-park-fire-new-mexico-where-wildfire-spread-ute-park-map-NFSA)).\n", "\n", "\n", "In this tutorial we use `eodag` to recover a stack of data on the same tile and processing it with *SNAP* to get a burned area mask that we will display on a dynamic map to evaluate the extent of damage.\n", "\n", "To be able to follow this tutorial, you will need to install:\n", "\n", "* [Sen2Cor](http://step.esa.int/main/snap-supported-plugins/sen2cor/): the installation instructions of the version used in this notebook (2.5.5) are provided in [this online PDF file](http://step.esa.int/thirdparties/sen2cor/2.5.5/docs/S2-PDGS-MPC-L2A-SRN-V2.5.5.pdf). To use *Sen2Cor* and its `L2A_Process` command you need to have its `bin` folder in your system *PATH*.\n", "* [SNAP](http://step.esa.int/main/download/snap-download/) which we will use through its Python interface called `snappy`.\n", "\n", " > **Warning:** DO NOT INSTALL the package `snappy` you can get from PyPi with `python -m pip install snappy`. It's not the one we want to use in this tutorial.\n", " \n", " > In order to set up the right `snappy` package follow the [instructions provided by SNAP](https://senbox.atlassian.net/wiki/spaces/SNAP/pages/50855941/Configure+Python+to+use+the+SNAP-Python+snappy+interface). Also, make sure that the version of Python you are using to run this notebook is supported by `snappy` (this notebook was created with Python 3.6 and SNAP 8.0).\n", " \n", " > **Note:** If you need support to install or configure SNAP or SNAPPY, please refer to https://forum.step.esa.int/c/snap.\n", "\n", "* The following Python packages: [Folium](https://python-visualization.github.io/folium/installing.html), [imageio](https://imageio.readthedocs.io/en/stable/installation.html), [Matplotlib](https://matplotlib.org/3.3.3/users/installing.html), [NumPy](https://numpy.org/install/), [Pillow](https://pillow.readthedocs.io/en/stable/installation.html) and [rasterio](https://rasterio.readthedocs.io/en/latest/installation.html).\n", "\n", "## Configuration\n", "\n", "Let's start by setting your personal credentials to access [PEPS service](https://peps.cnes.fr) by filling your username and password below:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "os.environ[\"EODAG__PEPS__AUTH__CREDENTIALS__USERNAME\"] = \"PLEASE_CHANGE_ME\"\n", "os.environ[\"EODAG__PEPS__AUTH__CREDENTIALS__PASSWORD\"] = \"PLEASE_CHANGE_ME\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you don't have the `bin` folders of *Sen2Cor* and *SNAP* in your system *PATH* uncomment the following lines of code, adapt the paths to your installations and run it to check whether the paths were correctly prepended to your *PATH*." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Add absolute path to Sen2Cor bin folder to make the L2A_Process command available\n", "#os.environ[\"PATH\"] = \"PLEASE_CHANGE_ME\" + \":\" + os.environ[\"PATH\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you haven't yet configured `snappy` to be available from your Python interpreter you can uncomment the following code and adapt the path to your installation." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import sys\n", "# Add the path to the snappy package which should be \"/home/yourusername/.snap/snap-python\"\n", "sys.path.insert(0, \"/home/maxime/.snap/snap-python/\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check that the Python packages required to run this notebook are available:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import folium\n", "import imageio\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import rasterio\n", "\n", "import snappy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a workspace directory where all our files and configuration will live. We also start the [EODataAccessGateway](../../api_reference/core.rst#eodag.api.core.EODataAccessGateway) instance that we will use to search and download data from *PEPS*. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2021-01-11 23:38:06,984-15s eodag.config [INFO ] Loading user configuration from: /home/maxime/TRAVAIL/06_EODAG/eodag/examples/eodag_workspace_fire_snappy/eodag_conf.yml\n", "INFO:eodag.config:Loading user configuration from: /home/maxime/TRAVAIL/06_EODAG/eodag/examples/eodag_workspace_fire_snappy/eodag_conf.yml\n", "2021-01-11 23:38:07,263-15s eodag.core [INFO ] Locations configuration loaded from /home/maxime/.config/eodag/locations.yml\n", "INFO:eodag.core:Locations configuration loaded from /home/maxime/.config/eodag/locations.yml\n" ] } ], "source": [ "from eodag.api.core import EODataAccessGateway\n", "from eodag import setup_logging\n", "\n", "setup_logging(verbose=2)\n", "\n", "# Create the workspace folder.\n", "workspace = 'eodag_workspace_fire_snappy'\n", "if not os.path.isdir(workspace):\n", " os.mkdir(workspace)\n", "\n", "# Save the PEPS configuration file.\n", "yaml_content = \"\"\"\n", "peps:\n", " download:\n", " outputs_prefix: \"{}\"\n", " extract: true\n", "\"\"\".format(os.path.abspath(workspace))\n", "\n", "with open(os.path.join(workspace, 'eodag_conf.yml'), \"w\") as f_yml:\n", " f_yml.write(yaml_content.strip())\n", " \n", "dag = EODataAccessGateway(os.path.join(workspace, 'eodag_conf.yml'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Search and Download with `eodag`\n", "\n", "We define the type of product we want to work on and the search area." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [], "source": [ "product_type = 'S2_MSI_L1C'\n", "extent = {\n", " 'lonmin': -105.095901,\n", " 'lonmax': -104.977112,\n", " 'latmin': 36.500253,\n", " 'latmax': 36.559015\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to frame the fire period by studying two images before and after, at the end of May and at the beginning of June." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2021-01-11 23:38:08,839-15s eodag.core [INFO ] Searching product type 'S2_MSI_L1C' on provider: peps\n", "INFO:eodag.core:Searching product type 'S2_MSI_L1C' on provider: peps\n", "2021-01-11 23:38:08,840-15s eodag.plugins.search.qssearch [INFO ] Sending count request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-05-21&completionDate=2018-05-25&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=1&page=1\n", "INFO:eodag.plugins.search.qssearch:Sending count request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-05-21&completionDate=2018-05-25&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=1&page=1\n", "2021-01-11 23:38:09,196-15s eodag.plugins.search.qssearch [INFO ] Sending search request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-05-21&completionDate=2018-05-25&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=20&page=1\n", "INFO:eodag.plugins.search.qssearch:Sending search request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-05-21&completionDate=2018-05-25&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=20&page=1\n", "2021-01-11 23:38:09,572-15s eodag.core [INFO ] Found 2 result(s) on provider 'peps'\n", "INFO:eodag.core:Found 2 result(s) on provider 'peps'\n" ] }, { "data": { "text/plain": [ "[EOProduct(id=S2A_MSIL1C_20180524T174051_N0206_R098_T13SEA_20180524T225950, provider=peps), EOProduct(id=S2A_MSIL1C_20180524T174051_N0206_R098_T13SDA_20180524T225950, provider=peps)]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "products_before, _ = dag.search(\n", " productType=product_type,\n", " start='2018-05-21', \n", " end='2018-05-25',\n", " geom=extent,\n", ")\n", "products_before" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2021-01-11 23:38:09,581-15s eodag.core [INFO ] Searching product type 'S2_MSI_L1C' on provider: peps\n", "INFO:eodag.core:Searching product type 'S2_MSI_L1C' on provider: peps\n", "2021-01-11 23:38:09,582-15s eodag.plugins.search.qssearch [INFO ] Sending count request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-06-06&completionDate=2018-06-10&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=1&page=1\n", "INFO:eodag.plugins.search.qssearch:Sending count request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-06-06&completionDate=2018-06-10&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=1&page=1\n", "2021-01-11 23:38:09,870-15s eodag.plugins.search.qssearch [INFO ] Sending search request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-06-06&completionDate=2018-06-10&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=20&page=1\n", "INFO:eodag.plugins.search.qssearch:Sending search request: https://peps.cnes.fr/resto/api/collections/S2ST/search.json?startDate=2018-06-06&completionDate=2018-06-10&geometry=POLYGON ((-105.0959 36.5003, -105.0959 36.5590, -104.9771 36.5590, -104.9771 36.5003, -105.0959 36.5003))&productType=S2MSI1C&maxRecords=20&page=1\n", "2021-01-11 23:38:10,327-15s eodag.core [INFO ] Found 2 result(s) on provider 'peps'\n", "INFO:eodag.core:Found 2 result(s) on provider 'peps'\n" ] }, { "data": { "text/plain": [ "[EOProduct(id=S2B_MSIL1C_20180608T173859_N0206_R098_T13SDA_20180608T211610, provider=peps), EOProduct(id=S2B_MSIL1C_20180608T173859_N0206_R098_T13SEA_20180608T211610, provider=peps)]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "products_after, _ = dag.search(\n", " productType=product_type,\n", " start='2018-06-06',\n", " end='2018-06-10',\n", " geom=extent,\n", ")\n", "products_after" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have now two sets of products which overlaps our region of interest, one before the fire and one after. We would like to work with the products that cover the most the burned area. eodag provides a handy utility that can group the products obtained from multiple searches by their extent, i.e. their bounding box. We will use this utility to find out which tiles we want to work with:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[EOProduct(id=S2A_MSIL1C_20180524T174051_N0206_R098_T13SEA_20180524T225950, provider=peps), EOProduct(id=S2B_MSIL1C_20180608T173859_N0206_R098_T13SEA_20180608T211610, provider=peps)],\n", " [EOProduct(id=S2A_MSIL1C_20180524T174051_N0206_R098_T13SDA_20180524T225950, provider=peps), EOProduct(id=S2B_MSIL1C_20180608T173859_N0206_R098_T13SDA_20180608T211610, provider=peps)]]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sorted_products = dag.group_by_extent([products_before, products_after])\n", "sorted_products" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(sorted_products)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The utility constructed two different groups which means that we can choose between two extents. As can be seen from their names, the tiles available are T13SEA and T13SDA. We display the two groups on a map with a tooltip (move your mouse over the polygons) to display their title:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "