To download these notebooks, please click "Download ZIP" to the top right:
Last active
September 9, 2024 04:01
-
-
Save robbibt/00aa3f3a1b8b5b5ea3f021b72f2af594 to your computer and use it in GitHub Desktop.
AEO 2024 DEA workshop material
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Burnt area mapping using Sentinel-2 data <img align=\"right\" src=\"../Supplementary_data/dea_logo.jpg\">\n", | |
| "\n", | |
| "* **[Sign up to the DEA Sandbox](https://app.sandbox.dea.ga.gov.au/)** to run this notebook interactively from a browser\n", | |
| "* **Compatibility:** Notebook compatible with the `NCI` and `DEA Sandbox` environments\n", | |
| "* **Products used:** \n", | |
| "[ga_s2am_ard_3](https://explorer.dea.ga.gov.au/ga_s2am_ard_3), \n", | |
| "[ga_s2bm_ard_3](https://explorer.dea.ga.gov.au/ga_s2bm_ard_3)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Background\n", | |
| "\n", | |
| "### Normalized Burn Ratio\n", | |
| "\n", | |
| "The Normalized Burn Ratio (NBR) is an index that uses the differences in the way healthy green vegetation and burnt vegetation reflect light to find burnt area. \n", | |
| "It is calculated using the following Sentinel-2 bands: Near Infrared/Band 8 and Shortwave Infrared/Band 12. The equation is defined below: \n", | |
| "\n", | |
| "\\begin{equation}\n", | |
| "NBR = \\frac{(NIR - SWIR)}{(NIR + SWIR)}\n", | |
| "\\end{equation}\n", | |
| "\n", | |
| "NBR returns values between -1 and 1. \n", | |
| "**Healthy green vegetation will have a high NBR value while burnt vegetation will have a low value**. \n", | |
| "Areas of dry, brown vegetation or bare soil will also return lower NBR values than green vegetation. \n", | |
| "\n", | |
| "### Delta Normalized Burn Ratio\n", | |
| "\n", | |
| "Change in Normalized Burn Ratio - also called Delta Normalized Burn Ratio (dNBR) - is calculated by subtracting the post-fire NBR value from the baseline NBR value as defined in this equation:\n", | |
| "\n", | |
| "\\begin{equation}\n", | |
| "dNBR = NBR_{baseline} - NBR_{post fire}\n", | |
| "\\end{equation}\n", | |
| "\n", | |
| "The dNBR value can be more useful than the NBR alone to determine what is burnt as it shows change from the baseline state. \n", | |
| "**A burnt area will have a positive dNBR value** while an unburnt area will have a negative dNBR value or a value close to zero.\n", | |
| "\n", | |
| "dNBR can also be used to describe burn severity (although this notebook does not look at burn severity). \n", | |
| "A higher severity fire will burn more vegetation, resulting in a higher dNBR. More information on NBR, dNBR and using it to measure burn severity can be found [on the UN-SPIDER knowledge portal](http://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/in-detail/normalized-burn-ratio).\n", | |
| "\n", | |
| "### Defining Burnt From Unburnt Areas\n", | |
| "\n", | |
| "[Rahman et al. 2018](https://doi.org/10.1109/IGARSS.2018.8518449) found a dNBR threshold value of +0.1 appropriate for differentiating burnt from unburnt areas when using Sentinel-2. \n", | |
| "However, some exploration with different threshold levels may be needed to get a good result in areas with different vegetation types. \n", | |
| "\n", | |
| "In the example presented in this notebook, which covers part of the Clyde Mountain fire in the area north of Batemans Bay, the fire occurred in heavily forested area, which returns a very strong dNBR result. \n", | |
| "Using +0.1 as a threshold here results in many false positives being picked up in the unburnt urban and forest areas where vegetation drying has occurred prior to the fire. \n", | |
| "A much more conservative threshold here of +0.3 produces a better result. \n", | |
| "Keep in mind the limitations of remote sensing and that in an ideal situation ground truth data collected in the location of interest would be used to assist in selecting a threshold.\n", | |
| "\n", | |
| "Some care should also be taken when interpreting results as a number of possible false positives can return a positive dNBR result:\n", | |
| "\n", | |
| "* A lot of smoke in the post burn image can interfere with the dNBR value\n", | |
| "* Areas that have been cleared of vegetation by some other means (logging, harvesting, and landslides) towards the end of the baseline period may incorrectly show up as burnt\n", | |
| "* Drying out of bright green vegetation such as grasses. \n", | |
| "If a fire event has been preceded by a rapid drying out of vegetation this can result in a low positive dNBR value in areas that have not burnt.\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Description\n", | |
| "\n", | |
| "This notebook calculates the change in Normalized Burn Ratio between a baseline composite image of the pre-fire condition of a selected area and a post-fire event image, in order to find burnt area extent. \n", | |
| " \n", | |
| "The user can change the location over which this notebook is run and specify a different date between which pre and post fire condition will be compared. \n", | |
| "The length of time over which the baseline composite image will be generated can be specified as 3, 6 or 12 months. \n", | |
| "The code in this notebook will automatically generate the composite image over the set length of time using Sentinel-2 data..\n", | |
| "\n", | |
| "The notebook contains the following steps:\n", | |
| "\n", | |
| "1. Select a location for the analysis\n", | |
| "2. Define fire event date and length of composite image\n", | |
| "3. Load all baseline data\n", | |
| "4. Generate Normalized Burn Ratio for baseline period\n", | |
| "5. Load post-fire data\n", | |
| "6. Generate Normalized Burn Ratio for post fire image\n", | |
| "7. Calculate Delta Normalized Burn Ratio\n", | |
| "8. Apply threshold to Delta Normalized Burn Ratio\n", | |
| "9. Calculate the area burnt\n", | |
| "10. Export results as a GeoTIFF\n", | |
| "\n", | |
| "***" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Getting started\n", | |
| "\n", | |
| "To run this analysis, run all the cells in the notebook, starting with the \"Load packages\" cell. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Load packages\n", | |
| "Import Python packages that are used for the analysis." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import datacube\n", | |
| "from datacube.utils import cog\n", | |
| "from datetime import datetime\n", | |
| "from datetime import timedelta\n", | |
| "import pandas as pd\n", | |
| "import geopandas as gpd\n", | |
| "import xarray as xr\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "import sys\n", | |
| "sys.path.insert(1, '../Tools/')\n", | |
| "from dea_tools.datahandling import load_ard\n", | |
| "from dea_tools.plotting import rgb, display_map, plot_variable_images\n", | |
| "from dea_tools.bandindices import calculate_indices\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Connect to the datacube\n", | |
| "\n", | |
| "Connect to the datacube so we can access DEA data.\n", | |
| "The `app` parameter is a unique name for the analysis which is based on the notebook file name." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "dc = datacube.Datacube(app=\"Burnt_area_mapping\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Select location\n", | |
| "\n", | |
| "The selected latitude and longitude will be displayed as a red box on the map below the next cell. \n", | |
| "This map can be used to find coordinates of other places, simply scroll and click on any point on the map to display the latitude and longitude of that location." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<div style=\"width:100%;\"><div style=\"position:relative;width:100%;height:0;padding-bottom:60%;\"><span style=\"color:#565656\">Make this Notebook Trusted to load map: File -> Trust Notebook</span><iframe srcdoc=\"<!DOCTYPE html>\n", | |
| "<html>\n", | |
| "<head>\n", | |
| " \n", | |
| " <meta http-equiv="content-type" content="text/html; charset=UTF-8" />\n", | |
| " \n", | |
| " <script>\n", | |
| " L_NO_TOUCH = false;\n", | |
| " L_DISABLE_3D = false;\n", | |
| " </script>\n", | |
| " \n", | |
| " <style>html, body {width: 100%;height: 100%;margin: 0;padding: 0;}</style>\n", | |
| " <style>#map {position:absolute;top:0;bottom:0;right:0;left:0;}</style>\n", | |
| " <script src="https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.js"></script>\n", | |
| " <script src="https://code.jquery.com/jquery-3.7.1.min.js"></script>\n", | |
| " <script src="https://cdn.jsdelivr.net/npm/bootstrap@5.2.2/dist/js/bootstrap.bundle.min.js"></script>\n", | |
| " <script src="https://cdnjs.cloudflare.com/ajax/libs/Leaflet.awesome-markers/2.0.2/leaflet.awesome-markers.js"></script>\n", | |
| " <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.css"/>\n", | |
| " <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/bootstrap@5.2.2/dist/css/bootstrap.min.css"/>\n", | |
| " <link rel="stylesheet" href="https://netdna.bootstrapcdn.com/bootstrap/3.0.0/css/bootstrap.min.css"/>\n", | |
| " <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/@fortawesome/fontawesome-free@6.2.0/css/all.min.css"/>\n", | |
| " <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/Leaflet.awesome-markers/2.0.2/leaflet.awesome-markers.css"/>\n", | |
| " <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/python-visualization/folium/folium/templates/leaflet.awesome.rotate.min.css"/>\n", | |
| " \n", | |
| " <meta name="viewport" content="width=device-width,\n", | |
| " initial-scale=1.0, maximum-scale=1.0, user-scalable=no" />\n", | |
| " <style>\n", | |
| " #map_45155595842c264a3aee6291495f9d96 {\n", | |
| " position: relative;\n", | |
| " width: 100.0%;\n", | |
| " height: 100.0%;\n", | |
| " left: 0.0%;\n", | |
| " top: 0.0%;\n", | |
| " }\n", | |
| " .leaflet-container { font-size: 1rem; }\n", | |
| " </style>\n", | |
| " \n", | |
| "</head>\n", | |
| "<body>\n", | |
| " \n", | |
| " \n", | |
| " <div class="folium-map" id="map_45155595842c264a3aee6291495f9d96" ></div>\n", | |
| " \n", | |
| "</body>\n", | |
| "<script>\n", | |
| " \n", | |
| " \n", | |
| " var map_45155595842c264a3aee6291495f9d96 = L.map(\n", | |
| " "map_45155595842c264a3aee6291495f9d96",\n", | |
| " {\n", | |
| " center: [-35.85, 136.95],\n", | |
| " crs: L.CRS.EPSG3857,\n", | |
| " zoom: 10,\n", | |
| " zoomControl: true,\n", | |
| " preferCanvas: false,\n", | |
| " }\n", | |
| " );\n", | |
| "\n", | |
| " \n", | |
| "\n", | |
| " \n", | |
| " \n", | |
| " var tile_layer_3500d9ef1c455baccf3acb14c7ae4516 = L.tileLayer(\n", | |
| " "http://mt1.google.com/vt/lyrs=y\\u0026z={z}\\u0026x={x}\\u0026y={y}",\n", | |
| " {"attribution": "Google", "detectRetina": false, "maxNativeZoom": 18, "maxZoom": 18, "minZoom": 0, "noWrap": false, "opacity": 1, "subdomains": "abc", "tms": false}\n", | |
| " );\n", | |
| " \n", | |
| " \n", | |
| " tile_layer_3500d9ef1c455baccf3acb14c7ae4516.addTo(map_45155595842c264a3aee6291495f9d96);\n", | |
| " \n", | |
| " \n", | |
| " var poly_line_5f1699d695032ba702682bd11e5cbba6 = L.polyline(\n", | |
| " [[-36.0, 136.8], [-36.0, 137.1], [-35.7, 137.1], [-35.7, 136.8], [-36.0, 136.8]],\n", | |
| " {"bubblingMouseEvents": true, "color": "red", "dashArray": null, "dashOffset": null, "fill": false, "fillColor": "red", "fillOpacity": 0.2, "fillRule": "evenodd", "lineCap": "round", "lineJoin": "round", "noClip": false, "opacity": 0.8, "smoothFactor": 1.0, "stroke": true, "weight": 3}\n", | |
| " ).addTo(map_45155595842c264a3aee6291495f9d96);\n", | |
| " \n", | |
| " \n", | |
| " var lat_lng_popup_291fa675bbc52c02c4bb2659b01a045c = L.popup();\n", | |
| " function latLngPop(e) {\n", | |
| " lat_lng_popup_291fa675bbc52c02c4bb2659b01a045c\n", | |
| " .setLatLng(e.latlng)\n", | |
| " .setContent("Latitude: " + e.latlng.lat.toFixed(4) +\n", | |
| " "<br>Longitude: " + e.latlng.lng.toFixed(4))\n", | |
| " .openOn(map_45155595842c264a3aee6291495f9d96);\n", | |
| " }\n", | |
| " map_45155595842c264a3aee6291495f9d96.on('click', latLngPop);\n", | |
| " \n", | |
| "</script>\n", | |
| "</html>\" style=\"position:absolute;width:100%;height:100%;left:0;top:0;border:none !important;\" allowfullscreen webkitallowfullscreen mozallowfullscreen></iframe></div></div>" | |
| ], | |
| "text/plain": [ | |
| "<folium.folium.Map at 0x7f3f380d4cd0>" | |
| ] | |
| }, | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "# Set the latitude and longitude\n", | |
| "ymax, xmin = -35.7, 136.8\n", | |
| "ymin, xmax = -36.0, 137.1\n", | |
| "\n", | |
| "from_date = \"2019-12-01\"\n", | |
| "to_date = \"2020-02-01\"\n", | |
| "\n", | |
| "# Set the buffer to load around the central coordinates\n", | |
| "buffer = 0.07\n", | |
| "\n", | |
| "# Compute the bounding box for the study area\n", | |
| "study_area_lat = (ymin, ymax)\n", | |
| "study_area_lon = (xmin, xmax)\n", | |
| "\n", | |
| "display_map(x=study_area_lon, y=study_area_lat, margin=-0.2)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Define fire event date and length of composite image\n", | |
| "\n", | |
| "Delta Normalized Burn Ratio produces the best result when using a post-fire image that was collected before much re-growth has occured. \n", | |
| "However, images collected while the fire is still active can be obscured by smoke and not show the full burn extent. \n", | |
| "As a result some adjustment of the fire event date entered may be needed to get the best result.\n", | |
| "\n", | |
| "The length of the baseline period can be automatically set to `3, 6 or 12 months`" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# Fire event date\n", | |
| "fire_date = '2019-12-01'\n", | |
| "\n", | |
| "# Length of baseline period\n", | |
| "baseline_length = '12 months'" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "#### Automaticaly define date range for baseline composite image" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "start_date_pre: 2018-12-01\n", | |
| "end_date_pre: 2019-11-30\n", | |
| "fire_date: 2019-12-01\n", | |
| "start_date_post: 2019-12-02\n", | |
| "end_date_post: 2020-02-29\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "# Define dates for loading data\n", | |
| "if baseline_length == '12 months':\n", | |
| " time_step = timedelta(days=365)\n", | |
| "if baseline_length == '6 months':\n", | |
| " time_step = timedelta(days=182.5)\n", | |
| "if baseline_length == '3 months':\n", | |
| " time_step = timedelta(days=91)\n", | |
| "\n", | |
| "# Calculate the start and end date for baseline data load\n", | |
| "start_date_pre = datetime.strftime(\n", | |
| " ((datetime.strptime(fire_date, '%Y-%m-%d')) - time_step), '%Y-%m-%d')\n", | |
| "end_date_pre = datetime.strftime(\n", | |
| " ((datetime.strptime(fire_date, '%Y-%m-%d')) - timedelta(days=1)),\n", | |
| " '%Y-%m-%d')\n", | |
| "\n", | |
| "# Calculate end date for post fire data load\n", | |
| "start_date_post = datetime.strftime(\n", | |
| " ((datetime.strptime(fire_date, '%Y-%m-%d')) + timedelta(days=1)),\n", | |
| " '%Y-%m-%d')\n", | |
| "end_date_post = datetime.strftime(\n", | |
| " ((datetime.strptime(fire_date, '%Y-%m-%d')) + timedelta(days=90)),\n", | |
| " '%Y-%m-%d')\n", | |
| "\n", | |
| "\n", | |
| "# Print dates\n", | |
| "print(f'start_date_pre: {start_date_pre}')\n", | |
| "print(f'end_date_pre: {end_date_pre}')\n", | |
| "print(f'fire_date: {fire_date}')\n", | |
| "print(f'start_date_post: {start_date_post}')\n", | |
| "print(f'end_date_post: {end_date_post}')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "#### Define load parameters" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "resolution = (-20, 20)\n", | |
| "measurements = ['nbart_blue', 'nbart_green', 'nbart_red',\n", | |
| " 'nbart_nir_1', 'nbart_swir_3']\n", | |
| "min_gooddata = 0.8\n", | |
| "output_crs = 'EPSG:3577'" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Load all baseline data" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Finding datasets\n", | |
| " ga_s2am_ard_3\n", | |
| " ga_s2bm_ard_3\n", | |
| "Counting good quality pixels for each time step using fmask\n", | |
| "Filtering to 17 out of 73 time steps with at least 80.0% good quality pixels\n", | |
| "Applying fmask pixel quality/cloud mask\n", | |
| "Loading 17 time steps\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "# Load all data in baseline period available from ARD data\n", | |
| "baseline = load_ard(dc=dc,\n", | |
| " products=['ga_s2am_ard_3', 'ga_s2bm_ard_3'],\n", | |
| " x=study_area_lon,\n", | |
| " y=study_area_lat,\n", | |
| " time=(start_date_pre, end_date_pre),\n", | |
| " measurements=measurements,\n", | |
| " min_gooddata=min_gooddata,\n", | |
| " output_crs=output_crs,\n", | |
| " resolution=resolution,\n", | |
| " group_by='solar_day')\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Generate Normalized Burn Ratio for baseline period" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# Calculate NBR for the baseline images\n", | |
| "baseline = calculate_indices(baseline, \n", | |
| " index='NBR', \n", | |
| " collection='ga_s2_3', \n", | |
| " drop=False)\n", | |
| "\n", | |
| "# Compute median using all observations in the dataset along the time axis\n", | |
| "baseline_image = baseline.median(dim='time')\n", | |
| "\n", | |
| "# Select NBR\n", | |
| "baseline_NBR = baseline_image.NBR" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Plot the baseline NBR data side-by-side with an RGB plot of the study area:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "tags": [] | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { |
View raw
(Sorry about that, but we can’t show files that are this big right now.)
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment