{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Comparison of dNBR for three wild fires in 2018 and 2019\n", "\n", "This notebook contains the calculations necessary for mapping the burn severity of the previously selected products.\n", "This burn severity is measured as the Normalized Burn Ratio (NBR).\n", "The development of this burn ratio over time is described using the ${\\Delta}$NBR, which is defined as:\n", "\n", "$$\n", "{\\Delta}\\text{NBR} = \\text{NBR}_\\text{pre-fire} - \\text{NBR}_\\text{post-fire}\n", "$$\n", "\n", "Instead of only calculating values by overlaying different bands in one product, the ${\\Delta}$NBR adds a time dimension to determine changes that have occurred in the time span between the two chosen products.\n", "\n", "## Methodology\n", "\n", "The NBR of the products downloaded in the previous notebook is calculated using the notebook [02 Spectral Indices](\"02 Spectral Indices.ipynb\").\n", "These NBR values are plotted side-by-side to get a visual impression of the changes that may have occurred between two given dates.\n", "The ${\\Delta}$NBR is calculated and, if showing signs of burned areas, compared to data given by the responsible German ministries.\n", "\n", "## Calculcations\n", "### Setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import geopandas as gpd\n", "\n", "from sentinel_helpers import geodataframe_on_map, scihub_band_paths\n", "nbr_path = Path('output/spectral_indices/')\n", "product_path = Path('input/forest_fires/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Treuenbrietzen (Postdam-Mittelmark) 2018\n", "\n", "The first case is the wild fire in August 2018 that led to the evacuation of three villages (source).\n", "We load the pre-calculated NBR values and geometries we fetched from OpenStreetMap:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "potsdam_mittelmark_nbr_2018 = list(sorted(nbr_path.glob('*2018*NBR*.tif')))\n", "potsdam_mittelmark_nbr_2018" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `postdam_mittelmark_geom` is a collection of three points:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
place_idosm_typeosm_iddisplay_nameplace_rankcategorytypeimportanceicongeometry
03285415node387079362Frohnsdorf, Treuenbrietzen, Potsdam-Mittelmark...19placevillage0.495https://nominatim.openstreetmap.org/ui/mapicon...POINT (12.90217 52.05466)
1554316node226935349Klausdorf, Treuenbrietzen, Potsdam-Mittelmark,...19placevillage0.495https://nominatim.openstreetmap.org/ui/mapicon...POINT (12.94213 52.04879)
2303964452node282202396Tiefenbrunnen, Treuenbrietzen, Potsdam-Mittelm...22placeisolated_dwelling0.420NonePOINT (12.94439 52.03532)
\n", "
" ], "text/plain": [ " place_id osm_type osm_id \\\n", "0 3285415 node 387079362 \n", "1 554316 node 226935349 \n", "2 303964452 node 282202396 \n", "\n", " display_name place_rank category \\\n", "0 Frohnsdorf, Treuenbrietzen, Potsdam-Mittelmark... 19 place \n", "1 Klausdorf, Treuenbrietzen, Potsdam-Mittelmark,... 19 place \n", "2 Tiefenbrunnen, Treuenbrietzen, Potsdam-Mittelm... 22 place \n", "\n", " type importance \\\n", "0 village 0.495 \n", "1 village 0.495 \n", "2 isolated_dwelling 0.420 \n", "\n", " icon \\\n", "0 https://nominatim.openstreetmap.org/ui/mapicon... \n", "1 https://nominatim.openstreetmap.org/ui/mapicon... \n", "2 None \n", "\n", " geometry \n", "0 POINT (12.90217 52.05466) \n", "1 POINT (12.94213 52.04879) \n", "2 POINT (12.94439 52.03532) " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "potsdam_mittelmark_geom = gpd.read_file(product_path / 'evacuated_2018.json')\n", "# geodataframe_on_map(postdam_mittelmark_geom)\n", "potsdam_mittelmark_geom" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### NBR plots\n", "\n", "We plot the NBR values for all dates side by side to get a first visual impression of changes that have occurred." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sentinel_helpers import scihub_band_paths, scihub_band_date, RasterReaderList\n", "import rasterio as r\n", "import rasterio.plot as rplt\n", "\n", "# we can save some time by reading only parts of the product we are interested in\n", "from rasterio.features import geometry_window\n", "\n", "def plot_nbrs(products, geom):\n", " with RasterReaderList(products) as readers:\n", " fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(20, 20))\n", "\n", " # we need to reproject from WGS84 so the geometry can be correctly plotted on the map\n", " _geom = geom.to_crs(readers[0].crs)\n", "\n", " # we don't need the entire NBR, we only use a slice\n", " window = geometry_window(readers[0], _geom.buffer(5000))\n", " window_transform = readers[0].window_transform(window)\n", "\n", " nbrs = np.array([reader.read(masked=True, window=window) for reader in readers])\n", " \n", " # ensure that subplot colors are chosen on a fixed scale\n", " vmin = nbrs.min()\n", " vmax = nbrs.max()\n", " \n", " for reader, nbr, ax in zip(readers, nbrs, axes):\n", " product_dt = scihub_band_date(reader)\n", " ax.tick_params(axis='x', labelrotation=90)\n", " rplt.show(nbr,\n", " ax=ax,\n", " transform=window_transform,\n", " cmap='Greys',\n", " title=product_dt.strftime('%Y-%m-%d'),\n", " vmin=vmin,\n", " vmax=vmax)\n", " _geom.plot(ax=ax, facecolor='none', edgecolor='red')\n", " \n", " # increase horizontal whitespace between subplots\n", " plt.subplots_adjust(wspace=0.32)\n", " \n", " # add colorbar using the last image\n", " img = axes[-1].get_images()[0]\n", " fig.colorbar(img, ax=axes, shrink=0.2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Red dots mark the location of the evacuated villages, the background shows the calculated NBR values, where areas with a low NBR value are brighter and areas with a high NBR value are darker." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "ename": "IndexError", "evalue": "list index out of range", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplot_nbrs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpotsdam_mittelmark_nbr_2018\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpotsdam_mittelmark_geom\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m\u001b[0m in \u001b[0;36mplot_nbrs\u001b[0;34m(products, geom)\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[0;31m# we need to reproject from WGS84 so the geometry can be correctly plotted on the map\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m \u001b[0m_geom\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgeom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_crs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreaders\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcrs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;31m# we don't need the entire NBR, we only use a slice\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mIndexError\u001b[0m: list index out of range" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABIkAAARiCAYAAAAgMacZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAolklEQVR4nO3dX4jl53nY8e9TKYLmT+MQb0IqyUQtchxd2MXZKKY0rdPQRnIvRCAXckJMTUCYRiGXFoUmF75pLgrBxIkQRgjfRBeNSZSixBRK4oLrViuwZSvGYatQa6uApTik4ECF7LcXZ6ZMNqOzZ6WZ3TnjzwcG9pzzauadl5nfY74+58ystQIAAADgW9vfudkbAAAAAODmE4kAAAAAEIkAAAAAEIkAAAAASCQCAAAAIJEIAAAAgHaIRDPz+Mx8dWa++DqPz8x8dGYuz8xzM/Puk98mAGeVOQHANuYEwP7Y5ZlET1T3bXn8/urug4+Hqt9689sCYI88kTkBwOt7InMCYC9cMxKttT5dfW3LkgeqT6yNz1ZvmZkfOKkNAnC2mRMAbGNOAOyPk3hPoturF4/cvnJwHwCUOQHAduYEwBlx6wl8jjnmvnXswpmH2jyFtO/4ju/4kXe84x0n8OUBzpdnn332lbXWhZu9jxNkTgCcIHPCnADY5s3MiZOIRFeqO4/cvqN66biFa63HqseqLl68uC5dunQCXx7gfJmZ/3Wz93DCzAmAE2ROmBMA27yZOXESLzd7qvrAwV8leE/1V2utPz+BzwvA+WBOALCNOQFwRlzzmUQz89vVe6u3zsyV6lerb6taaz1aPV29r7pc/XX1wdPaLABnjzkBwDbmBMD+uGYkWmu9/xqPr+oXT2xHAOwVcwKAbcwJgP1xEi83AwAAAGDPiUQAAAAAiEQAAAAAiEQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAA0I6RaGbum5kvz8zlmXnkmMe/e2Z+f2Y+PzPPz8wHT36rAJxV5gQA25gTAPvhmpFoZm6pPlbdX91TvX9m7rlq2S9Wf7LWelf13uo/zMxtJ7xXAM4gcwKAbcwJgP2xyzOJ7q0ur7VeWGu9Wj1ZPXDVmlV918xM9Z3V16rXTnSnAJxV5gQA25gTAHtil0h0e/XikdtXDu476jeqH65eqr5Q/fJa65tXf6KZeWhmLs3MpZdffvkNbhmAM8acAGAbcwJgT+wSieaY+9ZVt3+q+lz196t/VP3GzPy9v/UfrfXYWuviWuvihQsXrnOrAJxR5gQA25gTAHtil0h0pbrzyO072hT+oz5YfXJtXK7+rHrHyWwRgDPOnABgG3MCYE/sEomeqe6embsO3jzuweqpq9Z8pfrJqpn5/uqHqhdOcqMAnFnmBADbmBMAe+LWay1Ya702Mw9Xn6puqR5faz0/Mx86ePzR6iPVEzPzhTZPJ/3wWuuVU9w3AGeEOQHANuYEwP64ZiSqWms9XT191X2PHvn3S9W/PNmtAbAvzAkAtjEnAPbDLi83AwAAAOCcE4kAAAAAEIkAAAAAEIkAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAAaMdINDP3zcyXZ+byzDzyOmveOzOfm5nnZ+aPT3abAJxl5gQA25gTAPvh1mstmJlbqo9V/6K6Uj0zM0+ttf7kyJq3VL9Z3bfW+srMfN8p7ReAM8acAGAbcwJgf+zyTKJ7q8trrRfWWq9WT1YPXLXmZ6tPrrW+UrXW+urJbhOAM8ycAGAbcwJgT+wSiW6vXjxy+8rBfUe9vfqemfmjmXl2Zj5w3CeamYdm5tLMXHr55Zff2I4BOGvMCQC2MScA9sQukWiOuW9ddfvW6keqf1X9VPXvZubtf+s/WuuxtdbFtdbFCxcuXPdmATiTzAkAtjEnAPbENd+TqE3pv/PI7Tuql45Z88pa6+vV12fm09W7qj89kV0CcJaZEwBsY04A7Ildnkn0THX3zNw1M7dVD1ZPXbXm96ofn5lbZ+bbqx+rvnSyWwXgjDInANjGnADYE9d8JtFa67WZebj6VHVL9fha6/mZ+dDB44+utb40M39YPVd9s/r4WuuLp7lxAM4GcwKAbcwJgP0xa139cuAb4+LFi+vSpUs35WsDnGUz8+xa6+LN3sfNZk4AHM+c2DAnAI73ZubELi83AwAAAOCcE4kAAAAAEIkAAAAAEIkAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAoB0j0czcNzNfnpnLM/PIlnU/OjPfmJmfObktAnDWmRMAbGNOAOyHa0aimbml+lh1f3VP9f6Zued11v1a9amT3iQAZ5c5AcA25gTA/tjlmUT3VpfXWi+stV6tnqweOGbdL1W/U331BPcHwNlnTgCwjTkBsCd2iUS3Vy8euX3l4L7/b2Zur366enTbJ5qZh2bm0sxcevnll693rwCcTeYEANuYEwB7YpdINMfct666/evVh9da39j2idZaj621Lq61Ll64cGHHLQJwxpkTAGxjTgDsiVt3WHOluvPI7Tuql65ac7F6cmaq3lq9b2ZeW2v97klsEoAzzZwAYBtzAmBP7BKJnqnunpm7qv9dPVj97NEFa627Dv89M09U/8kFHeBbhjkBwDbmBMCeuGYkWmu9NjMPt/krA7dUj6+1np+ZDx08vvV1wwCcb+YEANuYEwD7Y5dnErXWerp6+qr7jr2Yr7X+9ZvfFgD7xJwAYBtzAmA/7PLG1QAAAACccyIRAAAAACIRAAAAACIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAACJRAAAAAAkEgEAAACQSAQAAABAIhEAAAAAiUQAAAAAJBIBAAAAkEgEAAAAQCIRAAAAAO0YiWbmvpn58sxcnplHjnn852bmuYOPz8zMu05+qwCcVeYEANuYEwD74ZqRaGZuqT5W3V/dU71/Zu65atmfVf9srfXO6iPVYye9UQDOJnMCgG3MCYD9scszie6tLq+1XlhrvVo9WT1wdMFa6zNrrb88uPnZ6o6T3SYAZ5g5AcA25gTAntglEt1evXjk9pWD+17PL1R/cNwDM/PQzFyamUsvv/zy7rsE4CwzJwDYxpwA2BO7RKI55r517MKZn2hzUf/wcY+vtR5ba11ca128cOHC7rsE4CwzJwDYxpwA2BO37rDmSnXnkdt3VC9dvWhm3ll9vLp/rfUXJ7M9APaAOQHANuYEwJ7Y5ZlEz1R3z8xdM3Nb9WD11NEFM/O26pPVz6+1/vTktwnAGWZOALCNOQGwJ675TKK11msz83D1qeqW6vG11vMz86GDxx+tfqX63uo3Z6bqtbXWxdPbNgBnhTkBwDbmBMD+mLWOfTnwqbt48eK6dOnSTfnaAGfZzDzrfxibEwCvx5zYMCcAjvdm5sQuLzcDAAAA4JwTiQAAAAAQiQAAAAAQiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABIJAIAAAAgkQgAAACARCIAAAAAEokAAAAASCQCAAAAIJEIAAAAgEQiAAAAABKJAAAAAEgkAgAAACCRCAAAAIBEIgAAAAASiQAAAABox0g0M/fNzJdn5vLMPHLM4zMzHz14/LmZeffJbxWAs8qcAGAbcwJgP1wzEs3MLdXHqvure6r3z8w9Vy27v7r74OOh6rdOeJ8AnFHmBADbmBMA+2OXZxLdW11ea72w1nq1erJ64Ko1D1SfWBufrd4yMz9wwnsF4GwyJwDYxpwA2BO7RKLbqxeP3L5ycN/1rgHgfDInANjGnADYE7fusGaOuW+9gTXNzENtnj5a9X9n5os7fP3z7q3VKzd7EzeZM9hwDhvOoX7oZm/gOpkTp8vvhDM45Bw2nIM5YU78TX4nnMEh57DhHN7EnNglEl2p7jxy+47qpTewprXWY9VjVTNzaa118bp2ew45B2dwyDlsOIfNGdzsPVwnc+IUOQdncMg5bDgHcyJz4m9wDs7gkHPYcA5vbk7s8nKzZ6q7Z+aumbmterB66qo1T1UfOPirBO+p/mqt9edvdFMA7BVzAoBtzAmAPXHNZxKttV6bmYerT1W3VI+vtZ6fmQ8dPP5o9XT1vupy9dfVB09vywCcJeYEANuYEwD7Y5eXm7XWerrNhfvofY8e+feqfvE6v/Zj17n+vHIOzuCQc9hwDnt4BubEqXIOzuCQc9hwDnt4BubEqXIOzuCQc9hwDm/iDGZzPQYAAADgW9ku70kEAAAAwDl36pFoZu6bmS/PzOWZeeSYx2dmPnrw+HMz8+7T3tONtsMZ/NzB9/7czHxmZt51M/Z52q51DkfW/ejMfGNmfuZG7u9G2eUcZua9M/O5mXl+Zv74Ru/xtO3wO/HdM/P7M/P5gzM4l+9LMDOPz8xXX+/P934rXB/LnChz4pA5sWFOmBNlRhxlTpgTh8yJDXPCnKhTnBNrrVP7aPPGdP+z+gfVbdXnq3uuWvO+6g+qqd5T/ffT3NON/tjxDP5x9T0H/77/vJ3BrudwZN1/afOa9Z+52fu+ST8Pb6n+pHrbwe3vu9n7vgln8G+rXzv494Xqa9VtN3vvp3AW/7R6d/XF13n8XF8fr+Pn4Vyfgzmx+zkcWWdOmBPnfk6YEdf183Cuz8Kc2P0cjqwzJ8wJc+INXhtP+5lE91aX11ovrLVerZ6sHrhqzQPVJ9bGZ6u3zMwPnPK+bqRrnsFa6zNrrb88uPnZ6o4bvMcbYZefhapfqn6n+uqN3NwNtMs5/Gz1ybXWV6rWWuftLHY5g1V918xM9Z1tLuqv3dhtnr611qfbfG+v57xfH8ucKHPikDmxYU6YE5UZcYQ5YU4cMic2zAlzojq9OXHakej26sUjt68c3He9a/bZ9X5/v9Cm9p031zyHmbm9+unq0c6vXX4e3l59z8z80cw8OzMfuGG7uzF2OYPfqH64eqn6QvXLa61v3pjtnSnn/fpY5kSZE4fMiQ1zwpzY1Xm/Nh4yJ8yJQ+bEhjlhTuzqDV0bbz217WzMMfdd/efUdlmzz3b+/mbmJ9pc1P/Jqe7o5tjlHH69+vBa6xub4Hsu7XIOt1Y/Uv1k9Xer/zYzn11r/elpb+4G2eUMfqr6XPXPq39Y/eeZ+a9rrf9zyns7a8779bHMiTInDpkTG+aEObGr835tPGROmBOHzIkNc8Kc2NUbujaediS6Ut155PYdbUre9a7ZZzt9fzPzzurj1f1rrb+4QXu7kXY5h4vVkwcX9LdW75uZ19Zav3tDdnhj7Po78cpa6+vV12fm09W7qvNyUd/lDD5Y/fu1eTHt5Zn5s+od1f+4MVs8M8779bHMiTInDpkTG+aEObGr835tPGROmBOHzIkNc8Kc2NUbujae9svNnqnunpm7Zua26sHqqavWPFV94OCdt99T/dVa689PeV830jXPYGbeVn2y+vlzVHevds1zWGvdtdb6wbXWD1b/sfo35+yCXrv9Tvxe9eMzc+vMfHv1Y9WXbvA+T9MuZ/CVNv/PRzPz/dUPVS/c0F2eDef9+ljmRJkTh8yJDXPCnNjVeb82HjInzIlD5sSGOWFO7OoNXRtP9ZlEa63XZubh6lNt3oH88bXW8zPzoYPHH23zrvPvqy5Xf92m+J0bO57Br1TfW/3mQfV+ba118Wbt+TTseA7n3i7nsNb60sz8YfVc9c3q42utY/+s4T7a8WfhI9UTM/OFNk+T/PBa65WbtulTMjO/Xb23euvMXKl+tfq2+ta4PpY5UebEIXNiw5wwJw6ZERvmhDlxyJzYMCfMiUOnNSdm8+wrAAAAAL6VnfbLzQAAAADYAyIRAAAAACIRAAAAACIRAAAAAIlEAAAAACQSAQAAAJBIBAAAAEAiEQAAAADV/wO0caPkz8390QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_nbrs(potsdam_mittelmark_nbr_2018, potsdam_mittelmark_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The NBR can be interpreted as follows:\n", "\n", "> Healthy vegetation shows a very high reflectance in the NIR[^nir], and low reflectance in the SWIR[^swir] portion of the spectrum - the opposite of what is seen in areas devastated by fire. Recently burnt areas demonstrate low reflectance in the NIR and high reflectance in the SWIR, i.e. the difference between the spectral responses of healthy vegetation and burnt areas reach their peak in the NIR and the SWIR regions of the spectrum.\n", ">\n", "> ![image.png](https://un-spider.org/sites/default/files/Spectral_responses.jpg) \n", "> Figure 2. Comparison of the spectral response of healthy vegetation and burned areas. Source: U.S. Forest service.\n", "> \n", "> To benefit from the magnitude of spectral difference, NBR uses the ratio between NIR and SWIR bands, according to the formula shown below. A high NBR value indicates healthy vegetation while a low value indicates bare ground and recently burnt areas. Non-burnt areas are normally attributed to values close to zero.\n", ">\n", "> $\\text{NBR} = \\frac{\\text{NIR} - \\text{SWIR}}{\\text{NIR} + \\text{SWIR}}$\n", "\n", "Source: [UN-Spider Knowledge Portal](https://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/in-detail/normalized-burn-ratio)\n", "\n", "[^nir]: Near Infra-Red\n", "[^swir]: Short-wave Infra-Red\n", "\n", "### dNBR calculation\n", "\n", "In order to calculate ${\\Delta}$NBR, the NBR after the fire is sutracted from the NBR before the fire.\n", "Reminder:\n", "\n", "$$\n", "{\\Delta}\\text{NBR} = \\text{NBR}_\\text{pre} - \\text{NBR}_\\text{post}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def calculate_dnbr(pre_nbr, post_nbr, geom=[]):\n", " with RasterReaderList([pre_nbr, post_nbr]) as readers:\n", " if len(geom):\n", " # if a geometry is passed, perform all calculations only in the\n", " # surroundings if this geometry\n", " _geom = geom.to_crs(readers[0].crs)\n", " window = geometry_window(readers[0], _geom.buffer(5000))\n", " window_transform = readers[0].window_transform(window)\n", " else:\n", " window = window_transform = None\n", " \n", " pre_fire = readers[0].read(masked=True, window=window)\n", " post_fire = readers[1].read(masked=True, window=window)\n", "\n", " # we need to mask invalid pixels in any of the input files for the resulting file\n", " dnbr = pre_fire - post_fire\n", " dnbr.mask = pre_fire.mask | post_fire.mask\n", " \n", " return (dnbr, window, window_transform)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> dNBR values can vary from case to case, and so, if possible, interpretation in specific instances should also be carried out through field assessment; in order to obtain the best results. However, the United States Geological Survey (USGS) proposed a classification table to interpret the burn severity, which can be seen below (Table 1).\n", ">\n", "> ![image.png](https://un-spider.org/sites/default/files/table+legend.PNG)\n", "> Table 1. Burn severity levels obtained calculating dNBR, proposed by USGS.\n", "\n", "Source: [UN-Spider Knowledge Portal](https://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/in-detail/normalized-burn-ratio)\n", "\n", "We construct a custom color scheme for matplotlib to match the values given in the table above and define a function that helps us plot the $\\Delta$NBR values using different color scales:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib.colors import BoundaryNorm, ListedColormap\n", "\n", "# define discrete color scale based on table above\n", "dnbr_cmap = ListedColormap(['#778835', '#a7c050', '#0be344', '#f8fc11', '#f8b140', '#f8671a', '#a600d4'])\n", "boundaries = [-0.5, -0.25, -0.1, 0.1, 0.27, 0.44, 0.66, 1.3]\n", "dnbr_norm = BoundaryNorm(boundaries, dnbr_cmap.N, clip=True)\n", "\n", "def plot_dnbr(dnbr, dnbr_crs, dnbr_transform, geom):\n", " fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 8))\n", " _geom = geom.to_crs(dnbr_crs)\n", "\n", " # plot black and white image\n", " rplt.show(dnbr, ax=ax1, transform=dnbr_transform, cmap='Greys', vmin=boundaries[0], vmax=boundaries[-1])\n", " ax1.tick_params(axis='x', labelrotation=90)\n", " _geom.plot(ax=ax1, facecolor='none', edgecolor='red')\n", " \n", " # plot image using colormap from above\n", " rplt.show(dnbr, ax=ax2, transform=dnbr_transform, cmap=dnbr_cmap, norm=dnbr_norm)\n", " ax2.tick_params(axis='x', labelrotation=90)\n", " _geom.plot(ax=ax2, facecolor='none', edgecolor='black')\n", " \n", " fig.colorbar(ax1.get_images()[0], ax=ax1, shrink=0.2)\n", " fig.colorbar(ax2.get_images()[0], ax=ax2, shrink=0.2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with r.open(potsdam_mittelmark_nbr_2018[0]) as src:\n", " # we only need the metadata which is now available in `src`\n", " pass\n", "\n", "dnbr, window, window_transform = calculate_dnbr(potsdam_mittelmark_nbr_2018[1], potsdam_mittelmark_nbr_2018[2], potsdam_mittelmark_geom)\n", "plot_dnbr(dnbr, src.crs, window_transform, potsdam_mittelmark_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The left image shows a dark spot between the three locations marking a decrease in NBR values after the start of the fire on August 23.\n", "is very likely the result of a loss of vegetation due to the fire.\n", "However, we can also see other examples of vegetation loss in the bottom left and top right corner of both plots which, according to the burn severity scale, are classified as \"high severity\" burns.\n", "\n", "Given the regular shape of these vegetation loss, we may assume that they are the result of crop harvests or that they have other causes unrelated to the fire.\n", "This interpretation is **based on assumptions** that need to be verified by experts.\n", "\n", "The area selected is verified against the _Waldbrandstatistik_:\n", "\n", "> Die höchste mediale Aufmerksamkeit erreichte der am 13.08.2018 entstandene Waldbrand bei Treuenbrietzen (Obf. Dippmannsdorf). Ausgelöst an mehreren Zündstellen und durch Weltkriegsmunition sowie böige Winde verstärkt, erfasste der Brand schnell über 300 ha Wald. Mehrere Dörfer wurden evakuiert. Der Baumbestand auf der Fläche, zumeist Kiefern, ist fast komplett vernichtet. \n", "\n", "Source: [Waldbrandstatistik 2018, p. 2](https://forst.brandenburg.de/sixcms/media.php/9/wbra2018.pdf)\n", "\n", "Please note that the date given in the Waldbrandstatistik is different to the one [given by Deutsche Welle](https://www.dw.com/en/forest-fires-near-berlin-could-burn-for-days/a-45203042).\n", "It is likely to be the result of an uncorrected typing error (August 13 vs. August 23), which is corrected in the following pages of the report.\n", "\n", "We can extract the area between the three evacuated villages and calculate the burned area.\n", "The coordinates for the box used for the extract can be read from this plot that marks its axes using array coordinates:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import box\n", "from math import floor\n", "\n", "potsdam_mittelmark_extract = dnbr[0,500:800,420:810]\n", "\n", "img = plt.imshow(potsdam_mittelmark_extract, cmap='Greys')\n", "plt.colorbar(img)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now calculate all pixels in this extract that are at least a low burn severity.\n", "Each pixel represents an area of 100m².\n", "One ha is 10 000m².\n", "The conversion of numbers of pixels to ha is therefore\n", "\n", "$$\n", "\\frac{n_\\text{pixels} \\cdot 100}{10\\ 000} = \\frac{n_\\text{pixels}}{100}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "potsdam_mittelmark_burned_area = (potsdam_mittelmark_extract >= 0.1).sum() / 100\n", "potsdam_mittelmark_burned_area" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The area given by the [Waldbrandstatistik 2019, p. 3](https://forst.brandenburg.de/sixcms/media.php/9/wbra2018.pdf) is exactly 300ha.\n", "The absolute error is therefore:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "abs((potsdam_mittelmark_burned_area / 300) - 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An error of 0 means the calculated area matches the official figure exactly.\n", "The number above means that the calculated burned area is within 7.8% of the official figure.\n", "\n", "### Jüterbog 2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jueterbog_geom = gpd.read_file(product_path / 'jueterbog_2019.json')\n", "jueterbog_nbr_2019 = list(sorted(nbr_path.glob('T33UUT*201906*NBR*.tif')))\n", "jueterbog_nbr_2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# geodataframe_on_map(jueterbog_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### NBR comparison\n", "\n", "The three NBR values are again plotted side by side.\n", "The red line marks the boundary of the town of Jüterbog as retrieved from OpenStreetMap." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_nbrs(jueterbog_nbr_2019, jueterbog_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot on the right contains grey spots covering the area, which are clouds that block the view to the ground.\n", "The ${\\Delta}$NBR is calculated using the left and middle images.\n", "\n", "#### dNBR calculation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with r.open(jueterbog_nbr_2019[0]) as src:\n", " # we only need the metadata which is now available in `src`\n", " pass\n", "\n", "dnbr, window, window_transform = calculate_dnbr(jueterbog_nbr_2019[0], jueterbog_nbr_2019[1], jueterbog_geom)\n", "plot_dnbr(dnbr, src.crs, window_transform, jueterbog_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Towards the top of the area, many areas show an increase of NBR between the first and second capture date.\n", "Like in Treuenbrietzen, regularly shaped losses indicate non-fire-related value changes.\n", "\n", "Media reports note that the fire broke out former military site which is now part of a nature reserve.\n", "Using OpenStreetmap, we can find the exact boundaries of the nature reserve and use it to restrict the area in which we look at NBR values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sentinel_helpers import search_osm\n", "nature_reserve_geom = search_osm('NSG Forst Zinna-Jüterbog-Keilberg')\n", "# geodataframe_on_map(nature_reserve_geom)\n", "nature_reserve_geom" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_dnbr(dnbr, src.crs, window_transform, nature_reserve_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the former military site encloses an area with a high ${\\Delta}$NBR value.\n", "We can create a raster mask from its shape and use it to select only the area it encloses." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from rasterio.features import geometry_mask\n", "\n", "nature_reserve_reprojected = nature_reserve_geom.to_crs(src.crs).iloc[0]['geometry']\n", "nature_reserve_mask = geometry_mask([nature_reserve_geom.to_crs(src.crs).iloc[0]['geometry']],\n", " out_shape=(window.height, window.width),\n", " transform=window_transform)\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(nature_reserve_mask)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "burned_in_nature_reserve = dnbr.copy()\n", "burned_in_nature_reserve.mask = burned_in_nature_reserve.mask | nature_reserve_mask\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(burned_in_nature_reserve[0], cmap=dnbr_cmap, norm=dnbr_norm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After masking the array with the shape of the nature reserve, we only keep the pixels within that nature reserve that have at marks of at least a low severity burn." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "burned_in_nature_reserve.mask = burned_in_nature_reserve.mask | (burned_in_nature_reserve.data < 0.1)\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(burned_in_nature_reserve[0], cmap=dnbr_cmap, norm=dnbr_norm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The size of the burned area in ha is the sum of all unmasked (i.e. opaque) pixels divided by 100:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "opaque_pixels = ~burned_in_nature_reserve.mask\n", "burned_area_in_ha = opaque_pixels.sum() / 100\n", "burned_area_in_ha" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fire area for June 2019 is given as 746ha by the Waldbrandstatistik 2019 (p. 17).\n", "The absolute error of the detection above is given by the following figure:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "abs((burned_area_in_ha / 746) - 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please note that fire area and area with burn marks are two different things, but it is the most reliable available information available that is close to ground truth.\n", "Compared with the Waldbrandstatistik, the area we calculated has an error of approximately 25%." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lübtheen 2019\n", "\n", "The fire in Lübtheen in July 2019 is interesting because there was a flyover of Sentinel-2 while the fire was active and burning, which\n", "ESA used to publish an animation of the burning fire:\n", "\n", "> ![ESA Lübtheen Flyover Image](http://www.esa.int/var/esa/storage/images/esa_multimedia/images/2019/07/german_wildfire/19466843-1-eng-GB/German_wildfire_pillars.gif)\n", "\n", "Image Source: [ESA, 2019](http://www.esa.int/ESA_Multimedia/Images/2019/07/German_wildfire)\n", "\n", "The animation above contains a true-color rendering and a frame that highlights wavelengths in the infrared spectrum to create a faux-color image highlighting the burn-scar and decreasing the visibility of smoke.\n", "\n", "Using the same data we can get a glimpse of the site before, during and after the fire:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "luebtheen_geom = gpd.read_file(product_path / 'luebtheen_2019.json')\n", "luebtheen_esa_pre_fire = Path('input/forest_fires/S2A_MSIL2A_20190629T103031_N0212_R108_T32UPE_20190629T135351.zip')\n", "luebtheen_esa_active_fire = Path('input/forest_fires/S2B_MSIL2A_20190701T102029_N0212_R065_T32UPE_20190701T134657.zip')\n", "luebtheen_esa_post_fire = Path('input/forest_fires/S2A_MSIL2A_20190726T102031_N0213_R065_T32UPE_20190726T125507.zip')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import box\n", "from sentinel_helpers import RasterReaderList\n", "\n", "def preview_true_color(products, geometry):\n", " tci_paths = map(lambda p: scihub_band_paths(p, ['TCI'], '60m')[0], products)\n", " \n", " with RasterReaderList(tci_paths) as readers:\n", " fig, axes = plt.subplots(ncols=len(readers), figsize=(20,10))\n", " _geom = geometry.to_crs(readers[0].crs)\n", " \n", " window = geometry_window(readers[0], _geom.buffer(5000))\n", " window_transform = readers[0].window_transform(window)\n", " \n", " for src, ax in zip(readers, axes):\n", " capture_date = scihub_band_date(src.name)\n", " rplt.show(src.read(window=window),\n", " transform=window_transform,\n", " ax=ax,\n", " title=capture_date.strftime('%Y-%m-%d'))\n", " _geom.plot(ax=ax, facecolor='none', edgecolor='red')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "preview_true_color([luebtheen_esa_pre_fire, luebtheen_esa_active_fire, luebtheen_esa_post_fire], luebtheen_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The area highlighted in the ESA animation is positioned in the lower-right corner of the superimposed boundary of Lübtheen.\n", "\n", "Looking at these true-color image we can see that on June 29 and July 26 we get a very clear view of the area.\n", "We can also see that there is a lot of vegetation change that is likely not related to the fire, visible when comparing the top and bottom-left parts of the images.\n", "\n", "This vegetation change is likely also reflected in the NBR and $\\Delta$NBR values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "luebtheen_nbr_2019 = list(sorted(nbr_path.glob('T32UPE_2019*.tif')))\n", "luebtheen_nbr_2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_nbrs(luebtheen_nbr_2019, luebtheen_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The NBR rendering of the active fire again contains clouds as almost monotonous gray patches scattered across the image.\n", "\n", "The images overall appear noisier than the previous two, making it harder to make out a clear spot to isolate." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with r.open(luebtheen_nbr_2019[0]) as src:\n", " # we only need the metadata which is now available in `src`\n", " pass\n", "\n", "dnbr, window, window_transform = calculate_dnbr(luebtheen_nbr_2019[0], luebtheen_nbr_2019[2], luebtheen_geom)\n", "plot_dnbr(dnbr, src.crs, window_transform, luebtheen_geom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The image contains a lot of high $\\Delta$NBR values, many of which are unlikely to have been caused by events related to the fire.\n", "\n", "The [Bundesanstalt für Immobilienaufgaben](https://web.archive.org/web/20210226162714/https://www.bundesimmobilien.de/bundesforstbetrieb-trave-ueberwacht-jetzt-das-gelaende-9457f54f122de1a3) published a report after the fire in which they note that it, too, broke out on a former military site. Using a the [OpenStreetMap search interface](https://www.openstreetmap.org/search?query=L%C3%BCbtheen#map=11/53.2894/11.0779) we can find the name of the former military site and retrieve its geometry using the Nominatim API as before:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "luebtheen_military_site = search_osm('Lübtheener Heide (ehem. Truppenübungsplatz)')\n", "luebtheen_military_site" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_dnbr(dnbr, src.crs, window_transform, luebtheen_military_site)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The boundary of the military site exceeds the area for which the $\\Delta$NBR was calculated to the right.\n", "This part is discarded in further calculations because it is several kilometers away from the apparent fire site.\n", "\n", "Following the method above, we isolate the burned pixels within the military site's boundaries and compare their area with the figure mentioned by official sources." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "luebtheen_military_site_mask = geometry_mask([luebtheen_military_site.to_crs(src.crs).iloc[0]['geometry']],\n", " out_shape=(window.height, window.width),\n", " transform=window_transform)\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(luebtheen_military_site_mask)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "burned_in_military_site = dnbr.copy()\n", "burned_in_military_site.mask = burned_in_military_site.mask | luebtheen_military_site_mask\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(burned_in_military_site[0], cmap=dnbr_cmap, norm=dnbr_norm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# keep only the pixels with a burn serverity of at least 0.1\n", "burned_in_military_site.mask = burned_in_military_site.mask | (burned_in_military_site.data < 0.1)\n", "\n", "plt.figure(figsize=(10, 10))\n", "plt.imshow(burned_in_military_site[0], cmap=dnbr_cmap, norm=dnbr_norm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate the area of all unmasked pixels\n", "opaque_pixels = ~burned_in_military_site.mask\n", "burned_area_in_ha = opaque_pixels.sum() / 100\n", "burned_area_in_ha" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The _Bundesanstalt für Immobilien_ puts the affected burned area at 944ha.\n", "This results in the following absolute error:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "abs((burned_area_in_ha / 944) - 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.8.6" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }