{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Burned aread analysis for Charter Call 693 - Wildfire in Goseong - South Korea" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Objective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook does the burned area analysis using Sentinel-2 data for the Charter Call 693 - Wildfire in Goseong - South Korea " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import snappy\n", "\n", "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import sys\n", "import os\n", "sys.path.append(os.getcwd())\n", "from py_snap_helpers import op_help, get_operator_default_parameters, GraphProcessor\n", "\n", "from geopandas import GeoDataFrame\n", "import pandas as pd\n", "\n", "import cgi\n", "import cioppy\n", "\n", "from shapely.geometry import box\n", "\n", "from shapely.wkt import loads\n", "\n", "import numpy as np\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colors as colors\n", "\n", "import math\n", "%matplotlib inline\n", "import lxml.etree as etree\n", "import gdal\n", "from ipywidgets import HTML\n", "sys.path.append('/opt/OTB/lib/python')\n", "sys.path.append('/opt/OTB/lib/libfftw3.so.3')\n", "os.environ['OTB_APPLICATION_PATH'] = '/opt/OTB/lib/otb/applications'\n", "os.environ['LD_LIBRARY_PATH'] = '/opt/OTB/lib'\n", "os.environ['ITK_AUTOLOAD_PATH'] = '/opt/OTB/lib/otb/applications'\n", "os.environ['GDAL_DATA'] = '/opt/anaconda/share/gdal/'\n", "import otbApplication\n", "from os.path import exists\n", "from osgeo.gdalconst import GA_ReadOnly\n", "from struct import unpack\n", "from PIL import Image\n", "from PIL import ImageDraw\n", "from StringIO import StringIO\n", "from base64 import b64encode\n", "\n", "from ipyleaflet import *\n", "from shapely.ops import cascaded_union" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Charter activation for the wildfire in Goseong - South Korea" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This section queries the Charter calls catalogue using the freetext search to retrieve the event area and time of interest" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "charter_endpoint = 'https://catalog.terradue.com/cos2/series/activations/description'\n", "\n", "charter_search_params = dict()\n", "\n", "charter_search_params['q'] = 'Goseong'" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "ciop = cioppy.Cioppy()\n", "\n", "activation = ciop.search(end_point=charter_endpoint, \n", " params=charter_search_params,\n", " output_fields='identifier,self,wkt,startdate,link:results,title', \n", " model='GeoTime')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{'identifier': 'Call-693',\n", " 'link:results': 'https://catalog.terradue.com/cos2/series/call-693-areas/search',\n", " 'self': 'https://catalog.terradue.com//cos2/series/activations/search?format=atom&uid=Call-693',\n", " 'startdate': '2019-04-05T00:49:00.0000000Z',\n", " 'title': '[Act-602/Call-693] Fire in South Korea',\n", " 'wkt': 'POINT(128.50810625 38.24185)'}]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "activation" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "call_areas = GeoDataFrame(ciop.search(end_point=activation[0]['link:results'], \n", " params=[],\n", " output_fields='identifier,self,wkt,link:related', \n", " model='GeoTime'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are several call areas: a point and two polygons. " ] }, { "cell_type": "code", "execution_count": 6, "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", "
identifierlink:relatedselfwkt
0Call-693-Area_0https://catalog.terradue.com//cos2/series/call...POINT(128.50810625 38.24185)
1Call-693-Area_1https://catalog.terradue.com//cos2/series/call...POLYGON((128.418 38.2536,128.4275 38.1604,128....
2Call-693-Area_2https://catalog.terradue.com//cos2/series/call...POLYGON((128.49 38.26,128.49 38.2,128.61 38.2,...
\n", "
" ], "text/plain": [ " identifier link:related \\\n", "0 Call-693-Area_0 \n", "1 Call-693-Area_1 \n", "2 Call-693-Area_2 \n", "\n", " self \\\n", "0 https://catalog.terradue.com//cos2/series/call... \n", "1 https://catalog.terradue.com//cos2/series/call... \n", "2 https://catalog.terradue.com//cos2/series/call... \n", "\n", " wkt \n", "0 POINT(128.50810625 38.24185) \n", "1 POLYGON((128.418 38.2536,128.4275 38.1604,128.... \n", "2 POLYGON((128.49 38.26,128.49 38.2,128.61 38.2,... " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "call_areas" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "call_areas['wkt'] = call_areas['wkt'].apply(loads)" ] }, { "cell_type": "code", "execution_count": 8, "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", "
identifierlink:relatedselfwkt
0Call-693-Area_0https://catalog.terradue.com//cos2/series/call...POINT (128.50810625 38.24185)
1Call-693-Area_1https://catalog.terradue.com//cos2/series/call...POLYGON ((128.418 38.2536, 128.4275 38.1604, 1...
2Call-693-Area_2https://catalog.terradue.com//cos2/series/call...POLYGON ((128.49 38.26, 128.49 38.2, 128.61 38...
\n", "
" ], "text/plain": [ " identifier link:related \\\n", "0 Call-693-Area_0 \n", "1 Call-693-Area_1 \n", "2 Call-693-Area_2 \n", "\n", " self \\\n", "0 https://catalog.terradue.com//cos2/series/call... \n", "1 https://catalog.terradue.com//cos2/series/call... \n", "2 https://catalog.terradue.com//cos2/series/call... \n", "\n", " wkt \n", "0 POINT (128.50810625 38.24185) \n", "1 POLYGON ((128.418 38.2536, 128.4275 38.1604, 1... \n", "2 POLYGON ((128.49 38.26, 128.49 38.2, 128.61 38... " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "call_areas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a map to visualize the areas of interest" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d9d56a01a2f04bb5bbae4ff15059d9ea", "version_major": 2, "version_minor": 0 }, "text/plain": [ "TWFwKGJhc2VtYXA9eyd1cmwnOiAnaHR0cHM6Ly97c30udGlsZS5vcGVuc3RyZWV0bWFwLm9yZy97en0ve3h9L3t5fS5wbmcnLCAnbWF4X3pvb20nOiAxOSwgJ2F0dHJpYnV0aW9uJzogJ01hcCDigKY=\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for index, row in call_areas.iterrows():\n", "\n", " if row['wkt'].geom_type == 'Point':\n", " m = Map(center=(row['wkt'].y, \n", " row['wkt'].x), \n", " zoom=10)\n", " \n", " break\n", "\n", "marker = Marker(location=(row['wkt'].y,\n", " row['wkt'].x), draggable=False)\n", "\n", "d = {'title': activation[0]['title'], \n", " 'date' : activation[0]['startdate']}\n", " \n", "html_value=\"\"\"\n", "
\n", "
\"\"\".format(**d)\n", " \n", " \n", "html_widget_slave = HTML(\n", " value=html_value,\n", " placeholder='',\n", " description='',\n", " )\n", " \n", " \n", "marker.popup = html_widget_slave\n", "\n", "\n", "\n", "m.add_layer(marker);\n", " \n", "m" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([,\n", " ,\n", " ], dtype=object)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "call_areas['wkt'].values" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "aoi = cascaded_union(call_areas['wkt'].values)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'POLYGON ((128.5680275793651 38.26, 128.61 38.26, 128.61 38.2, 128.5942576964478 38.2, 128.6102 38.1604, 128.4275 38.1604, 128.418 38.2536, 128.4164 38.2971, 128.4259 38.2995, 128.5369 38.3037, 128.5680275793651 38.26))'" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aoi.wkt" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "min_lon, min_lat, max_lon, max_lat = list(aoi.bounds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add the area of interest to the map" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "p = Polygon(locations=np.asarray([t[::-1] for t in list(aoi.exterior.coords)]).tolist(), color=\"red\", fill_color=\"green\")\n", "\n", "layer_group = LayerGroup(layers=(m.layers), name='Call AOI')\n", "\n", "layer_group.add_layer(p)\n", "\n", "m.add_layer(layer_group)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Search for Sentinel-2 acquisitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Search for Sentinle-2 acquistions over the call AOI and over a period of two weeks centered on the activation call:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "sentinel2_endpoint = 'https://catalog.terradue.com/sentinel2/description'" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "search_params = dict()\n", "\n", "search_params['geom'] = aoi.wkt\n", "search_params['pt'] = 'S2MSI2A'\n", "search_params['start'] = (pd.to_datetime(activation[0]['startdate']) - np.timedelta64(7, 'D')).isoformat() + 'Z'\n", "search_params['stop'] = (pd.to_datetime(activation[0]['startdate']) + np.timedelta64(7, 'D')).isoformat() + 'Z'#(pd.to_datetime(activation[0]['startdate']) + np.timedelta64(1, 'D')).isoformat() + 'Z'\n", "search_params['do'] = 'terradue'\n", "search_params['cc'] = '20['" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'cc': '20[',\n", " 'do': 'terradue',\n", " 'geom': 'POLYGON ((128.5680275793651 38.26, 128.61 38.26, 128.61 38.2, 128.5942576964478 38.2, 128.6102 38.1604, 128.4275 38.1604, 128.418 38.2536, 128.4164 38.2971, 128.4259 38.2995, 128.5369 38.3037, 128.5680275793651 38.26))',\n", " 'pt': 'S2MSI2A',\n", " 'start': '2019-03-29T00:49:00Z',\n", " 'stop': '2019-04-12T00:49:00Z'}" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "search_params" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "sentinel2_search = GeoDataFrame(ciop.search(end_point=sentinel2_endpoint, \n", " params=search_params,\n", " output_fields='identifier,self,wkt,startdate,enddate,enclosure,orbitDirection,cc', \n", " model='EOP'))" ] }, { "cell_type": "code", "execution_count": 19, "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", "
ccenclosureenddateidentifierorbitDirectionselfstartdatewkt
014.461197https://store.terradue.com/download/sentinel2/...2019-04-08T02:16:09.0240000ZS2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-08T02:16:09.0240000ZPOLYGON((127.862827804024 37.8539506483739,129...
116.589038https://store.terradue.com/download/sentinel2/...2019-04-05T02:06:59.0240000ZS2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-05T02:06:59.0240000ZPOLYGON((128.234033401221 37.855568274556,129....
24.844317https://store.terradue.com/download/sentinel2/...2019-04-03T02:16:51.0240000ZS2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-03T02:16:51.0240000ZPOLYGON((127.862827804024 37.8539506483739,129...
\n", "
" ], "text/plain": [ " cc enclosure \\\n", "0 14.461197 https://store.terradue.com/download/sentinel2/... \n", "1 16.589038 https://store.terradue.com/download/sentinel2/... \n", "2 4.844317 https://store.terradue.com/download/sentinel2/... \n", "\n", " enddate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " identifier orbitDirection \\\n", "0 S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING \n", "1 S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING \n", "2 S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING \n", "\n", " self \\\n", "0 https://catalog.terradue.com/sentinel2/search?... \n", "1 https://catalog.terradue.com/sentinel2/search?... \n", "2 https://catalog.terradue.com/sentinel2/search?... \n", "\n", " startdate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " wkt \n", "0 POLYGON((127.862827804024 37.8539506483739,129... \n", "1 POLYGON((128.234033401221 37.855568274556,129.... \n", "2 POLYGON((127.862827804024 37.8539506483739,129... " ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel2_search" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "sentinel2_search['wkt'] = sentinel2_search['wkt'].apply(loads)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111',\n", " 'S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128',\n", " 'S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016'], dtype=object)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel2_search['identifier'].values" ] }, { "cell_type": "code", "execution_count": 22, "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", "
ccenclosureenddateidentifierorbitDirectionselfstartdatewkt
014.461197https://store.terradue.com/download/sentinel2/...2019-04-08T02:16:09.0240000ZS2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-08T02:16:09.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...
116.589038https://store.terradue.com/download/sentinel2/...2019-04-05T02:06:59.0240000ZS2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-05T02:06:59.0240000ZPOLYGON ((128.234033401221 37.855568274556, 12...
24.844317https://store.terradue.com/download/sentinel2/...2019-04-03T02:16:51.0240000ZS2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-03T02:16:51.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...
\n", "
" ], "text/plain": [ " cc enclosure \\\n", "0 14.461197 https://store.terradue.com/download/sentinel2/... \n", "1 16.589038 https://store.terradue.com/download/sentinel2/... \n", "2 4.844317 https://store.terradue.com/download/sentinel2/... \n", "\n", " enddate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " identifier orbitDirection \\\n", "0 S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING \n", "1 S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING \n", "2 S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING \n", "\n", " self \\\n", "0 https://catalog.terradue.com/sentinel2/search?... \n", "1 https://catalog.terradue.com/sentinel2/search?... \n", "2 https://catalog.terradue.com/sentinel2/search?... \n", "\n", " startdate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " wkt \n", "0 POLYGON ((127.862827804024 37.8539506483739, 1... \n", "1 POLYGON ((128.234033401221 37.855568274556, 12... \n", "2 POLYGON ((127.862827804024 37.8539506483739, 1... " ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel2_search" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def analyse(row, aoi):\n", " \n", " aoi_intersection = (aoi.intersection(row['wkt']).area / aoi.area) * 100\n", " \n", " series = dict([('aoi_intersection', aoi_intersection)])\n", " \n", " series['utm_zone'] = row['identifier'][39:41]\n", " series['latitude_band'] = row['identifier'][41]\n", " series['grid_square'] = row['identifier'][42:44]\n", "\n", " \n", " return pd.Series(series)\n", "\n", " " ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "sentinel2_search = sentinel2_search.merge(sentinel2_search.apply(lambda row: analyse(row, aoi), axis=1), \n", " left_index=True,\n", " right_index=True)" ] }, { "cell_type": "code", "execution_count": 25, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ccenclosureenddateidentifierorbitDirectionselfstartdatewktaoi_intersectiongrid_squarelatitude_bandutm_zone
014.461197https://store.terradue.com/download/sentinel2/...2019-04-08T02:16:09.0240000ZS2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-08T02:16:09.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...100.0DHS52
116.589038https://store.terradue.com/download/sentinel2/...2019-04-05T02:06:59.0240000ZS2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-05T02:06:59.0240000ZPOLYGON ((128.234033401221 37.855568274556, 12...100.0DHS52
24.844317https://store.terradue.com/download/sentinel2/...2019-04-03T02:16:51.0240000ZS2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-03T02:16:51.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...100.0DHS52
\n", "
" ], "text/plain": [ " cc enclosure \\\n", "0 14.461197 https://store.terradue.com/download/sentinel2/... \n", "1 16.589038 https://store.terradue.com/download/sentinel2/... \n", "2 4.844317 https://store.terradue.com/download/sentinel2/... \n", "\n", " enddate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " identifier orbitDirection \\\n", "0 S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING \n", "1 S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING \n", "2 S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING \n", "\n", " self \\\n", "0 https://catalog.terradue.com/sentinel2/search?... \n", "1 https://catalog.terradue.com/sentinel2/search?... \n", "2 https://catalog.terradue.com/sentinel2/search?... \n", "\n", " startdate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " wkt aoi_intersection \\\n", "0 POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 \n", "1 POLYGON ((128.234033401221 37.855568274556, 12... 100.0 \n", "2 POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 \n", "\n", " grid_square latitude_band utm_zone \n", "0 DH S 52 \n", "1 DH S 52 \n", "2 DH S 52 " ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel2_search" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Stage in the selected Sentinel-2 data" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "target_dir = '/workspace/data/'\n", "\n", "def stage_in(row):\n", " \n", " \n", " try: \n", " local_path = ciop.copy(row['enclosure'], extract=True, target=target_dir)\n", " \n", " except:\n", " local_path = os.path.join(target_dir, row['identifier'] )\n", "\n", " row['local_path'] = os.path.join(target_dir, row['identifier'])\n", " \n", " return row \n", " " ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "sentinel2_search = sentinel2_search.apply(lambda row: stage_in(row), axis=1)" ] }, { "cell_type": "code", "execution_count": 28, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ccenclosureenddateidentifierorbitDirectionselfstartdatewktaoi_intersectiongrid_squarelatitude_bandutm_zonelocal_path
014.461197https://store.terradue.com/download/sentinel2/...2019-04-08T02:16:09.0240000ZS2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-08T02:16:09.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...100.0DHS52/workspace/data/S2B_MSIL2A_20190408T021609_N02...
116.589038https://store.terradue.com/download/sentinel2/...2019-04-05T02:06:59.0240000ZS2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-05T02:06:59.0240000ZPOLYGON ((128.234033401221 37.855568274556, 12...100.0DHS52/workspace/data/S2B_MSIL2A_20190405T020659_N02...
24.844317https://store.terradue.com/download/sentinel2/...2019-04-03T02:16:51.0240000ZS2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2...DESCENDINGhttps://catalog.terradue.com/sentinel2/search?...2019-04-03T02:16:51.0240000ZPOLYGON ((127.862827804024 37.8539506483739, 1...100.0DHS52/workspace/data/S2A_MSIL2A_20190403T021651_N02...
\n", "
" ], "text/plain": [ " cc enclosure \\\n", "0 14.461197 https://store.terradue.com/download/sentinel2/... \n", "1 16.589038 https://store.terradue.com/download/sentinel2/... \n", "2 4.844317 https://store.terradue.com/download/sentinel2/... \n", "\n", " enddate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " identifier orbitDirection \\\n", "0 S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING \n", "1 S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING \n", "2 S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING \n", "\n", " self \\\n", "0 https://catalog.terradue.com/sentinel2/search?... \n", "1 https://catalog.terradue.com/sentinel2/search?... \n", "2 https://catalog.terradue.com/sentinel2/search?... \n", "\n", " startdate \\\n", "0 2019-04-08T02:16:09.0240000Z \n", "1 2019-04-05T02:06:59.0240000Z \n", "2 2019-04-03T02:16:51.0240000Z \n", "\n", " wkt aoi_intersection \\\n", "0 POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 \n", "1 POLYGON ((128.234033401221 37.855568274556, 12... 100.0 \n", "2 POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 \n", "\n", " grid_square latitude_band utm_zone \\\n", "0 DH S 52 \n", "1 DH S 52 \n", "2 DH S 52 \n", "\n", " local_path \n", "0 /workspace/data/S2B_MSIL2A_20190408T021609_N02... \n", "1 /workspace/data/S2B_MSIL2A_20190405T020659_N02... \n", "2 /workspace/data/S2A_MSIL2A_20190403T021651_N02... " ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel2_search" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Near Infrared and Short-Wave Infrared RGB Composite\n", "\n", " Using the Sentinel-2 bands B12, B11 and B8A, generate the RGB composite" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "bands = ['B12', 'B11', 'B8A']" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "def get_band_path(row, band):\n", " \n", " ns = {'xfdu': 'urn:ccsds:schema:xfdu:1',\n", " 'safe': 'http://www.esa.int/safe/sentinel/1.1',\n", " 'gml': 'http://www.opengis.net/gml'}\n", " \n", " path_manifest = os.path.join(row['local_path'],\n", " row['identifier'] + '.SAFE', \n", " 'manifest.safe')\n", " \n", " root = etree.parse(path_manifest)\n", " \n", " bands = [band]\n", "\n", " for index, band in enumerate(bands):\n", "\n", " sub_path = os.path.join(row['local_path'],\n", " row['identifier'] + '.SAFE',\n", " root.xpath('//dataObjectSection/dataObject/byteStream/fileLocation[contains(@href,(\"%s%s\")) and contains(@href,(\"%s\")) ]' % (row['latitude_band'],\n", " row['grid_square'], \n", " band), \n", " namespaces=ns)[0].attrib['href'][2:])\n", " \n", " return sub_path" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use OTB for the contrast enhancement:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "def contrast_enhancement(in_tif, out_tif, hfact=2.0):\n", "\n", " ContrastEnhancement = otbApplication.Registry.CreateApplication(\"ContrastEnhancement\")\n", "\n", " ContrastEnhancement.SetParameterString(\"in\", in_tif)\n", " ContrastEnhancement.SetParameterString(\"out\", out_tif)\n", " ContrastEnhancement.SetParameterOutputImagePixelType(\"out\", otbApplication.ImagePixelType_uint8)\n", " ContrastEnhancement.SetParameterFloat(\"nodata\", 0.0)\n", " ContrastEnhancement.SetParameterFloat(\"hfact\", hfact)\n", " ContrastEnhancement.SetParameterInt(\"bins\", 256)\n", " ContrastEnhancement.SetParameterInt(\"spatial.local.w\", 500)\n", " ContrastEnhancement.SetParameterInt(\"spatial.local.h\", 500)\n", " ContrastEnhancement.SetParameterString(\"mode\",\"lum\")\n", "\n", " ContrastEnhancement.ExecuteAndWriteOutput()\n" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "rgb_composites = []\n", "\n", "for index, row in sentinel2_search.iterrows():\n", " \n", " vrt_bands = []\n", " \n", " for j, band in enumerate(bands):\n", " \n", " vrt_bands.append(get_band_path(row, band))\n", " \n", " vrt = '{0}.vrt'.format(row['identifier'])\n", " ds = gdal.BuildVRT(vrt,\n", " vrt_bands,\n", " srcNodata=0,\n", " xRes=10, \n", " yRes=10,\n", " separate=True)\n", " ds.FlushCache()\n", " \n", " tif = '{0}.tif'.format(row['identifier'])\n", " gdal.Translate(tif,\n", " vrt,\n", " projWin=[min_lon, max_lat, max_lon, min_lat],\n", " projWinSRS='EPSG:4326',\n", " outputType=gdal.GDT_Byte, \n", " scaleParams=[[0,10000,0,255]])\n", " \n", " tif_e = '{0}_e.tif'.format(row['identifier'])\n", " \n", " contrast_enhancement(tif, tif_e)\n", "\n", " rgb_composites.append(tif_e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a new map with the RGB composites:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "m = Map(center=(aoi.centroid.y, \n", " aoi.centroid.x), \n", " zoom=13)\n", "\n", "m.add_control(LayersControl())" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "rgb_pngs = []\n", "\n", "layers = []\n", "\n", "for index, composite in enumerate(rgb_composites):\n", " \n", " im = Image.open(composite)\n", " \n", " f = StringIO()\n", "\n", " im.save(f, 'png')\n", " data = b64encode(f.getvalue())\n", "\n", " rgb_pngs.append('data:image/png;base64,' + data)\n", " \n", " image = ImageOverlay(\n", " url=rgb_pngs[index],\n", " bounds=((min_lat, min_lon), (max_lat, max_lon))\n", " )\n", "\n", " layer_group = LayerGroup(layers=(), name='RGB composite #{0}'.format(index))\n", "\n", " layer_group.add_layer(image)\n", "\n", " m.add_layer(layer_group)\n", "\n", " layers.append(layer_group)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ef3b23c1ce8b464aa0af3ef139772176", "version_major": 2, "version_minor": 0 }, "text/plain": [ "TWFwKGJhc2VtYXA9eyd1cmwnOiAnaHR0cHM6Ly97c30udGlsZS5vcGVuc3RyZWV0bWFwLm9yZy97en0ve3h9L3t5fS5wbmcnLCAnbWF4X3pvb20nOiAxOSwgJ2F0dHJpYnV0aW9uJzogJ01hcCDigKY=\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pre-processing for the burned area analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use a SNAP graph to calculate the NDWI (to descriminate the water) and the NBR:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "def pre_processing(**kwargs):\n", " \n", " options = dict()\n", " \n", " operators = ['Read', \n", " 'Resample',\n", " 'BandMaths',\n", " 'Reproject',\n", " 'Subset',\n", " 'Write']\n", " \n", " for operator in operators:\n", " \n", " print 'Getting default values for Operator {}'.format(operator)\n", " parameters = get_operator_default_parameters(operator)\n", " \n", " options[operator] = parameters\n", "\n", " for key, value in kwargs.items():\n", " \n", " print 'Updating Operator {}'.format(key)\n", " options[key.replace('_', '-')].update(value)\n", " \n", " mygraph = GraphProcessor()\n", " \n", " for index, operator in enumerate(operators):\n", " \n", " print 'Adding Operator {} to graph'.format(operator)\n", " if index == 0: \n", " source_node_id = ''\n", " \n", " else:\n", " source_node_id = operators[index - 1]\n", " \n", " mygraph.add_node(operator,\n", " operator, \n", " options[operator], source_node_id)\n", " \n", " mygraph.view_graph()\n", " \n", " mygraph.run()\n", " " ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/workspace/data/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.SAFE/MTD_MSIL2A.xml\n", "Getting default values for Operator Read\n", "Getting default values for Operator Resample\n", "Getting default values for Operator BandMaths\n", "Getting default values for Operator Reproject\n", "Getting default values for Operator Subset\n", "Getting default values for Operator Write\n", "Updating Operator Subset\n", "Updating Operator Read\n", "Updating Operator Resample\n", "Updating Operator BandMaths\n", "Updating Operator Write\n", "Updating Operator Reproject\n", "Adding Operator Read to graph\n", "Adding Operator Resample to graph\n", "Adding Operator BandMaths to graph\n", "Adding Operator Reproject to graph\n", "Adding Operator Subset to graph\n", "Adding Operator Write to graph\n", "\n", " 1.0\n", " \n", " Read\n", " \n", " \n", " SENTINEL-2-MSI-MultiRes-UTM52N\n", " /workspace/data/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.SAFE/MTD_MSIL2A.xml\n", " \n", " \n", " \n", " Resample\n", " \n", " \n", " \n", " \n", " \n", " First\n", " \n", " true\n", " First\n", " Nearest\n", " \n", " B2\n", " \n", " \n", " \n", " BandMaths\n", " \n", " \n", " \n", " \n", " \n", " \n", " NDWI\n", " float32\n", " (B3 - B8) / (B3 + B8)\n", " \n", " \n", " NaN\n", " \n", " \n", " NBR\n", " float32\n", " (B8 - B12) / (B8 + B12)\n", " \n", " \n", " NaN\n", " \n", " \n", " \n", " \n", " \n", " \n", " Reproject\n", " \n", " \n", " \n", " \n", " EPSG:4326\n", " false\n", " \n", " 0\n", " \n", " true\n", " \n", " Nearest\n", " \n", " \n", " \n", " \n", " \n", " \n", " false\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " Subset\n", " \n", " \n", " \n", " \n", " \n", " \n", " 1\n", " 1\n", " \n", " true\n", " false\n", " POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))\n", " \n", " \n", " \n", " Write\n", " \n", " \n", " \n", " \n", " false\n", " BEAM-DIMAP\n", " pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111\n", " true\n", " true\n", " \n", " \n", "\n", "\n", "Processing the graph\n", "Process PID: 3584\n", "('Executing processing graph\\n.15%.30%.45%.60%.75%.90% done.\\n', 'INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters\\nINFO: org.esa.s2tbx.dataio.s2.ortho.S2OrthoProductReaderPlugIn: Building product reader - EPSG:32652\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [metadata_level] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nINFO: org.hsqldb.persist.Logger: dataFileCache open start\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: missing file /workspace/data/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.SAFE/GRANULE/L2A_T52SDH_A010897_20190408T021607/QI_DATA/L2A_T52SDH_20190408T021609_DDV_20m.jp2\\n\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: no image files found for band quality_dense_dark_vegetation\\n\\n')\n", "Done.\n", "/workspace/data/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128.SAFE/MTD_MSIL2A.xml\n", "Getting default values for Operator Read\n", "Getting default values for Operator Resample\n", "Getting default values for Operator BandMaths\n", "Getting default values for Operator Reproject\n", "Getting default values for Operator Subset\n", "Getting default values for Operator Write\n", "Updating Operator Subset\n", "Updating Operator Read\n", "Updating Operator Resample\n", "Updating Operator BandMaths\n", "Updating Operator Write\n", "Updating Operator Reproject\n", "Adding Operator Read to graph\n", "Adding Operator Resample to graph\n", "Adding Operator BandMaths to graph\n", "Adding Operator Reproject to graph\n", "Adding Operator Subset to graph\n", "Adding Operator Write to graph\n", "\n", " 1.0\n", " \n", " Read\n", " \n", " \n", " SENTINEL-2-MSI-MultiRes-UTM52N\n", " /workspace/data/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128.SAFE/MTD_MSIL2A.xml\n", " \n", " \n", " \n", " Resample\n", " \n", " \n", " \n", " \n", " \n", " First\n", " \n", " true\n", " First\n", " Nearest\n", " \n", " B2\n", " \n", " \n", " \n", " BandMaths\n", " \n", " \n", " \n", " \n", " \n", " \n", " NDWI\n", " float32\n", " (B3 - B8) / (B3 + B8)\n", " \n", " \n", " NaN\n", " \n", " \n", " NBR\n", " float32\n", " (B8 - B12) / (B8 + B12)\n", " \n", " \n", " NaN\n", " \n", " \n", " \n", " \n", " \n", " \n", " Reproject\n", " \n", " \n", " \n", " \n", " EPSG:4326\n", " false\n", " \n", " 0\n", " \n", " true\n", " \n", " Nearest\n", " \n", " \n", " \n", " \n", " \n", " \n", " false\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " Subset\n", " \n", " \n", " \n", " \n", " \n", " \n", " 1\n", " 1\n", " \n", " true\n", " false\n", " POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))\n", " \n", " \n", " \n", " Write\n", " \n", " \n", " \n", " \n", " false\n", " BEAM-DIMAP\n", " pre_S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128\n", " true\n", " true\n", " \n", " \n", "\n", "\n", "Processing the graph\n", "Process PID: 3673\n", "('Executing processing graph\\n.15%.30%.45%.60%.75%.90% done.\\n', 'INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters\\nINFO: org.esa.s2tbx.dataio.s2.ortho.S2OrthoProductReaderPlugIn: Building product reader - EPSG:32652\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [metadata_level] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nINFO: org.hsqldb.persist.Logger: dataFileCache open start\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: missing file /workspace/data/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128.SAFE/GRANULE/L2A_T52SDH_A010854_20190405T020656/QI_DATA/L2A_T52SDH_20190405T020659_DDV_20m.jp2\\n\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: no image files found for band quality_dense_dark_vegetation\\n\\n')\n", "Done.\n", "/workspace/data/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.SAFE/MTD_MSIL2A.xml\n", "Getting default values for Operator Read\n", "Getting default values for Operator Resample\n", "Getting default values for Operator BandMaths\n", "Getting default values for Operator Reproject\n", "Getting default values for Operator Subset\n", "Getting default values for Operator Write\n", "Updating Operator Subset\n", "Updating Operator Read\n", "Updating Operator Resample\n", "Updating Operator BandMaths\n", "Updating Operator Write\n", "Updating Operator Reproject\n", "Adding Operator Read to graph\n", "Adding Operator Resample to graph\n", "Adding Operator BandMaths to graph\n", "Adding Operator Reproject to graph\n", "Adding Operator Subset to graph\n", "Adding Operator Write to graph\n", "\n", " 1.0\n", " \n", " Read\n", " \n", " \n", " SENTINEL-2-MSI-MultiRes-UTM52N\n", " /workspace/data/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.SAFE/MTD_MSIL2A.xml\n", " \n", " \n", " \n", " Resample\n", " \n", " \n", " \n", " \n", " \n", " First\n", " \n", " true\n", " First\n", " Nearest\n", " \n", " B2\n", " \n", " \n", " \n", " BandMaths\n", " \n", " \n", " \n", " \n", " \n", " \n", " NDWI\n", " float32\n", " (B3 - B8) / (B3 + B8)\n", " \n", " \n", " NaN\n", " \n", " \n", " NBR\n", " float32\n", " (B8 - B12) / (B8 + B12)\n", " \n", " \n", " NaN\n", " \n", " \n", " \n", " \n", " \n", " \n", " Reproject\n", " \n", " \n", " \n", " \n", " EPSG:4326\n", " false\n", " \n", " 0\n", " \n", " true\n", " \n", " Nearest\n", " \n", " \n", " \n", " \n", " \n", " \n", " false\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " Subset\n", " \n", " \n", " \n", " \n", " \n", " \n", " 1\n", " 1\n", " \n", " true\n", " false\n", " POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))\n", " \n", " \n", " \n", " Write\n", " \n", " \n", " \n", " \n", " false\n", " BEAM-DIMAP\n", " pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016\n", " true\n", " true\n", " \n", " \n", "\n", "\n", "Processing the graph\n", "Process PID: 3762\n", "('Executing processing graph\\n.15%.30%.45%.60%.75%.90% done.\\n', 'INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters\\nINFO: org.esa.s2tbx.dataio.s2.ortho.S2OrthoProductReaderPlugIn: Building product reader - EPSG:32652\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [metadata_level] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nWARNING: org.esa.s2tbx.dataio.metadata.GenericXmlMetadata: Metadata: the path to element [bandid] does not exist\\nINFO: org.hsqldb.persist.Logger: dataFileCache open start\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: missing file /workspace/data/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.SAFE/GRANULE/L2A_T52SDH_A019734_20190403T021647/QI_DATA/L2A_T52SDH_20190403T021651_DDV_20m.jp2\\n\\nWARNING: org.esa.s2tbx.dataio.s2.ortho.Sentinel2OrthoProductReader: Warning: no image files found for band quality_dense_dark_vegetation\\n\\n')\n", "Done.\n" ] } ], "source": [ "pre_processed = []\n", "\n", "resample = dict()\n", "resample['referenceBandName'] = 'B2'\n", "\n", "reproject = dict()\n", "reproject['crs'] = 'EPSG:4326'\n", "\n", "subset = dict()\n", "subset['geoRegion'] = box(*aoi.bounds).wkt\n", "subset['copyMetadata'] = 'true'\n", "\n", "bands = '''\n", " \n", " NDWI\n", " float32\n", " (B3 - B8) / (B3 + B8)\n", " \n", " \n", " NaN\n", " \n", " \n", " NBR\n", " float32\n", " (B8 - B12) / (B8 + B12)\n", " \n", " \n", " NaN\n", " \n", " '''\n", "\n", "band_maths = dict()\n", "band_maths['targetBandDescriptors'] = bands \n", "\n", "for index, row in sentinel2_search.iterrows():\n", " \n", " print os.path.join(row['local_path'], row['identifier'] + '.SAFE', 'MTD_MSIL2A.xml')\n", " \n", " read = dict()\n", " read['file'] = os.path.join(row['local_path'], row['identifier'] + '.SAFE', 'MTD_MSIL2A.xml') #, 'manifest.safe')\n", " read['formatName'] = 'SENTINEL-2-MSI-MultiRes-UTM52N'\n", " \n", " write = dict()\n", " write['file'] = 'pre_{}'.format(row['identifier'])\n", " \n", " pre_processed.append('pre_{}'.format(row['identifier']))\n", " \n", " pre_processing(Read=read, \n", " Resample=resample, \n", " Reproject=reproject, \n", " Subset=subset,\n", " BandMaths=band_maths,\n", " Write=write)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Collocation " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Collocate the Sentinel-2 processed results using as Master the pre-fire and as Slave the post-fire Sentinel-2 acquisitions" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111',\n", " 'pre_S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128',\n", " 'pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016']" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pre_processed" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "master_path = pre_processed[2] + '.dim'\n", "slave_path = pre_processed[0] + '.dim'" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "mygraph = GraphProcessor()" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "operator = 'Read'\n", "\n", "node_id = 'Read_M'\n", "\n", "source_node_id = ''\n", "\n", "parameters = get_operator_default_parameters(operator)\n", " \n", "parameters['file'] = master_path \n", " \n", "mygraph.add_node(node_id,\n", " operator,\n", " parameters, source_node_id)\n" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "operator = 'Read'\n", "\n", "node_id = 'Read_S'\n", "\n", "source_node_id = ''\n", "\n", "parameters = get_operator_default_parameters(operator)\n", " \n", "parameters['file'] = slave_path \n", " \n", "mygraph.add_node(node_id, \n", " operator,\n", " parameters, source_node_id)" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "operator = 'Collocate'\n", "\n", "parameters = get_operator_default_parameters(operator)\n", "\n", "parameters['masterComponentPattern'] = 'PRE_FIRE_${ORIGINAL_NAME}'\n", "parameters['slaveComponentPattern'] = 'POST_FIRE_${ORIGINAL_NAME}'\n", "\n", "source_node_id = dict()\n", "\n", "source_node_id['master'] = 'Read_M'\n", "\n", "source_node_id['slave'] = 'Read_S'\n", "\n", "\n", "node_id = 'Collocate'\n", "\n", "mygraph.add_node(operator,\n", " operator, \n", " parameters, source_node_id)" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "\n", "operator = 'Write'\n", "\n", "node_id = 'Write'\n", "\n", "source_node_id = 'Collocate'\n", "\n", "\n", "\n", "parameters = get_operator_default_parameters(operator)\n", "\n", "parameters['file'] = 'collocated'\n", "parameters['formatName'] = 'BEAM-DIMAP'\n", "\n", "mygraph.add_node( node_id, \n", " operator,\n", " parameters, source_node_id)" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " 1.0\n", " \n", " Read\n", " \n", " \n", " \n", " pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.dim\n", " \n", " \n", " \n", " Read\n", " \n", " \n", " \n", " pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.dim\n", " \n", " \n", " \n", " Collocate\n", " \n", " Read_S\n", " Read_M\n", " \n", " \n", " true\n", " POST_FIRE_${ORIGINAL_NAME}\n", " _collocated\n", " true\n", " COLLOCATED\n", " PRE_FIRE_${ORIGINAL_NAME}\n", " NEAREST_NEIGHBOUR\n", " \n", " \n", " \n", " Write\n", " \n", " \n", " \n", " \n", " true\n", " false\n", " BEAM-DIMAP\n", " true\n", " collocated\n", " \n", " \n", "\n", "\n" ] } ], "source": [ "mygraph.view_graph()" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Processing the graph\n", "Process PID: 4000\n", "('Executing processing graph\\n.12%.24%.36%.48%.60%.72%.84%. done.\\n', 'INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters\\nINFO: org.hsqldb.persist.Logger: dataFileCache open start\\n')\n", "Done.\n" ] } ], "source": [ "mygraph.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### dNBR, RBR and burned area severity" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "def burned_area(**kwargs):\n", " \n", " options = dict()\n", " \n", " operators = ['Read', \n", " 'BandMaths',\n", " 'Write']\n", " \n", " for operator in operators:\n", " \n", " print 'Getting default values for Operator {}'.format(operator)\n", " parameters = get_operator_default_parameters(operator)\n", " \n", " options[operator] = parameters\n", "\n", " for key, value in kwargs.items():\n", " \n", " print 'Updating Operator {}'.format(key)\n", " options[key.replace('_', '-')].update(value)\n", " \n", " mygraph = GraphProcessor()\n", " \n", " for index, operator in enumerate(operators):\n", " \n", " print 'Adding Operator {} to graph'.format(operator)\n", " if index == 0: \n", " source_node_id = ''\n", " \n", " else:\n", " source_node_id = operators[index - 1]\n", " \n", " mygraph.add_node(operator,\n", " operator, \n", " options[operator], source_node_id)\n", " \n", " mygraph.view_graph()\n", " \n", " mygraph.run()\n", " " ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Getting default values for Operator Read\n", "Getting default values for Operator BandMaths\n", "Getting default values for Operator Write\n", "Updating Operator Read\n", "Updating Operator Write\n", "Updating Operator BandMaths\n", "Adding Operator Read to graph\n", "Adding Operator BandMaths to graph\n", "Adding Operator Write to graph\n", "\n", " 1.0\n", " \n", " Read\n", " \n", " \n", " \n", " collocated.dim\n", " \n", " \n", " \n", " BandMaths\n", " \n", " \n", " \n", " \n", " \n", " \n", " dNBR\n", " float32\n", " PRE_FIRE_NDWI >= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) > 0.27 ? PRE_FIRE_NBR - POST_FIRE_NBR : 0 \n", " \n", " \n", " NaN\n", " \n", " \n", " RBR\n", " float32\n", " PRE_FIRE_NDWI >= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) > 0.27 ? ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) : 0\n", " \n", " \n", " NaN\n", " \n", " \n", " \n", " \n", " \n", " \n", " Write\n", " \n", " \n", " \n", " \n", " false\n", " GeoTIFF\n", " burned_area\n", " true\n", " true\n", " \n", " \n", "\n", "\n", "Processing the graph\n", "Process PID: 4083\n", "('Executing processing graph\\n.12%.24%.36%.48%.60%.72%.84%. done.\\n', 'INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters\\nINFO: org.hsqldb.persist.Logger: dataFileCache open start\\n')\n", "Done.\n" ] } ], "source": [ "collocated_input = 'collocated.dim'\n", "\n", "read = dict()\n", "read['file'] = collocated_input\n", "\n", "bands = '''\n", " \n", " dNBR\n", " float32\n", " PRE_FIRE_NDWI >= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) > 0.27 ? PRE_FIRE_NBR - POST_FIRE_NBR : 0 \n", " \n", " \n", " NaN\n", " \n", " \n", " RBR\n", " float32\n", " PRE_FIRE_NDWI >= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) > 0.27 ? ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) : 0\n", " \n", " \n", " NaN\n", " \n", " '''\n", "\n", "band_maths = dict()\n", "band_maths['targetBandDescriptors'] = bands \n", "\n", "write = dict()\n", "write['file'] = 'burned_area'\n", "write['formatName'] = 'GeoTIFF'\n", "\n", "burned_area(Read=read, \n", " BandMaths=band_maths,\n", " Write=write)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Burned area severity" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "def raster2rgb(raster_file, color_table, out_file_name, raster_band=1, discrete=True):\n", " \n", " #Reading the band\n", " data_types ={'Byte':'B','UInt16':'H','Int16':'h','UInt32':'I','Int32':'i','Float32':'f','Float64':'d'}\n", " if exists(raster_file) is False:\n", " raise Exception('[Errno 2] No such file or directory: \\'' + raster_file + '\\'') \n", " \n", " dataset = gdal.Open(raster_file, GA_ReadOnly )\n", " \n", " if dataset == None:\n", " raise Exception(\"Unable to read the data file\")\n", " \n", " geoTransform = dataset.GetGeoTransform()\n", " proj = dataset.GetProjection()\n", " \n", " band = dataset.GetRasterBand(raster_band)\n", " values = band.ReadRaster( 0, 0, band.XSize, band.YSize, band.XSize, band.YSize, band.DataType )\n", " values = unpack(data_types[gdal.GetDataTypeName(band.DataType)]*band.XSize*band.YSize,values)\n", " \n", " #Preparing the color table and the output file\n", " classification_values = color_table.keys()\n", " classification_values.sort()\n", " \n", " base = Image.new( 'RGBA', (band.XSize,band.YSize) )\n", " base_draw = ImageDraw.Draw(base)\n", " alpha_mask = Image.new('L', (band.XSize,band.YSize), 255)\n", " alpha_draw = ImageDraw.Draw(alpha_mask)\n", " \n", " #Reading the value and setting the output color for each pixel\n", " for pos in range(len(values)):\n", " y = pos/band.XSize\n", " x = pos - y * band.XSize\n", " for index in range(len(classification_values)):\n", "\n", " if values[pos] <= classification_values[index] or index == len(classification_values)-1:\n", " if discrete == True:\n", " if index == 0:\n", " index = 1\n", " elif index == len(classification_values)-1 and values[pos] >= classification_values[index]:\n", " index = index + 1\n", " color = color_table[classification_values[index-1]]\n", " base_draw.point((x,y), (color[0],color[1],color[2]))\n", " alpha_draw.point((x,y),color[3])\n", " else:\n", " if index == 0:\n", " r = color_table[classification_values[0]][0]\n", " g = color_table[classification_values[0]][1]\n", " b = color_table[classification_values[0]][2]\n", " a = color_table[classification_values[0]][3]\n", " elif index == len(classification_values)-1 and values[pos] >= classification_values[index]:\n", " r = color_table[classification_values[index]][0]\n", " g = color_table[classification_values[index]][1]\n", " b = color_table[classification_values[index]][2]\n", " a = color_table[classification_values[index]][3]\n", " else:\n", " r = color_table[classification_values[index-1]][0] + (values[pos] - classification_values[index-1])*(color_table[classification_values[index]][0] - color_table[classification_values[index-1]][0])/(classification_values[index]-classification_values[index-1]) \n", " g = color_table[classification_values[index-1]][1] + (values[pos] - classification_values[index-1])*(color_table[classification_values[index]][1] - color_table[classification_values[index-1]][1])/(classification_values[index]-classification_values[index-1]) \n", " b = color_table[classification_values[index-1]][2] + (values[pos] - classification_values[index-1])*(color_table[classification_values[index]][2] - color_table[classification_values[index-1]][2])/(classification_values[index]-classification_values[index-1]) \n", " a = color_table[classification_values[index-1]][3] + (values[pos] - classification_values[index-1])*(color_table[classification_values[index]][3] - color_table[classification_values[index-1]][3])/(classification_values[index]-classification_values[index-1]) \n", " \n", " base_draw.point((x,y), (int(r),int(g),int(b)))\n", " alpha_draw.point((x,y),int(a))\n", " \n", " break\n", " #Adding transparency and saving the output image \n", " color_layer = Image.new('RGBA', base.size, (255, 255, 255, 0))\n", " base = Image.composite(color_layer, base, alpha_mask)\n", " base.save(out_file_name)\n", "\n", " # update geolocation\n", " ds_rgb = gdal.Open(out_file_name,1)\n", " ds_rgb.SetGeoTransform(geoTransform)\n", " ds_rgb.SetProjection(proj)\n", " \n", " ds_rgb.FlushCache()\n", " \n", " ds_rgb = None" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "severity_palette = {-1: [159, 159, 159, 0], # grey\n", " -0.5: [43, 25, 223, 0], # blue\n", " -0.251: [139, 221, 231, 0], # cyan\n", " -0.101: [97, 169, 45, 255], # unburned, green \n", " 0.099: [250, 254, 76, 0], # yellow\n", " 0.269: [228, 173, 55, 0], # orange\n", " 0.439: [202, 59, 18, 0], # red\n", " 0.659: [85, 15, 112, 0]} # purple" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [], "source": [ "raster2rgb('burned_area.tif',\n", " severity_palette,\n", " 'burned_area.rgb.tif',\n", " raster_band=1,\n", " discrete=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### View the results on a map" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [], "source": [ "im = Image.open('burned_area.rgb.tif')\n", " \n", "f = StringIO()\n", "\n", "im.save(f, 'png')\n", "data = b64encode(f.getvalue())\n", " \n", "burned_area = 'data:image/png;base64,' + data" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [], "source": [ "m = Map(center=(aoi.centroid.y, \n", " aoi.centroid.x), \n", " zoom=13)\n", "\n", "m.add_control(LayersControl())" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [], "source": [ "image = ImageOverlay(\n", " url=burned_area,\n", " bounds=((min_lat, min_lon), (max_lat, max_lon))\n", ")\n", "\n", "layer_group = LayerGroup(layers=(m.layers), name='Burned area severity')\n", " \n", "layer_group.add_layer(image)\n", "\n", "m.add_layer(layer_group)\n", "\n", "# add the RGB layers\n", "for layer in layers:\n", " \n", " m.add_layer(layer)\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9a87960832fb4dc2a970c97c160b5974", "version_major": 2, "version_minor": 0 }, "text/plain": [ "TWFwKGJhc2VtYXA9eyd1cmwnOiAnaHR0cHM6Ly97c30udGlsZS5vcGVuc3RyZWV0bWFwLm9yZy97en0ve3h9L3t5fS5wbmcnLCAnbWF4X3pvb20nOiAxOSwgJ2F0dHJpYnV0aW9uJzogJ01hcCDigKY=\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "This work is licenced under a [Attribution-ShareAlike 4.0 International License (CC BY-SA 4.0)](http://creativecommons.org/licenses/by-sa/4.0/) \n", "\n", "YOU ARE FREE TO:\n", "\n", "* Share - copy and redistribute the material in any medium or format.\n", "* Adapt - remix, transform, and built upon the material for any purpose, even commercially.\n", "\n", "UNDER THE FOLLOWING TERMS:\n", "\n", "* Attribution - You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.\n", "* ShareAlike - If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:mypython2]", "language": "python", "name": "conda-env-mypython2-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.12" } }, "nbformat": 4, "nbformat_minor": 2 }