Burned aread analysis for Charter Call 693 - Wildfire in Goseong - South Korea

Objective

This notebook does the burned area analysis using Sentinel-2 data for the Charter Call 693 - Wildfire in Goseong - South Korea

[1]:

import snappy

%load_ext autoreload
%autoreload 2

import sys
import os
sys.path.append(os.getcwd())
from py_snap_helpers import op_help, get_operator_default_parameters, GraphProcessor

from geopandas import GeoDataFrame
import pandas as pd

import cgi
import cioppy

from shapely.geometry import box

from shapely.wkt import loads

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors

import math
%matplotlib inline
import lxml.etree as etree
import gdal
from ipywidgets import HTML
sys.path.append('/opt/OTB/lib/python')
sys.path.append('/opt/OTB/lib/libfftw3.so.3')
os.environ['OTB_APPLICATION_PATH'] = '/opt/OTB/lib/otb/applications'
os.environ['LD_LIBRARY_PATH'] = '/opt/OTB/lib'
os.environ['ITK_AUTOLOAD_PATH'] = '/opt/OTB/lib/otb/applications'
os.environ['GDAL_DATA'] = '/opt/anaconda/share/gdal/'
import otbApplication
from os.path import exists
from osgeo.gdalconst import GA_ReadOnly
from struct import unpack
from PIL import Image
from PIL import ImageDraw
from StringIO import StringIO
from base64 import b64encode

from ipyleaflet import *
from shapely.ops import cascaded_union

Charter activation for the wildfire in Goseong - South Korea

This section queries the Charter calls catalogue using the freetext search to retrieve the event area and time of interest

[2]:

charter_endpoint = 'https://catalog.terradue.com/cos2/series/activations/description'

charter_search_params = dict()

charter_search_params['q'] = 'Goseong'

[3]:

ciop = cioppy.Cioppy()

activation = ciop.search(end_point=charter_endpoint,
                        params=charter_search_params,
                        output_fields='identifier,self,wkt,startdate,link:results,title',
                        model='GeoTime')

[4]:

activation

[4]:

[{'identifier': 'Call-693',
  'link:results': 'https://catalog.terradue.com/cos2/series/call-693-areas/search',
  'self': 'https://catalog.terradue.com//cos2/series/activations/search?format=atom&uid=Call-693',
  'startdate': '2019-04-05T00:49:00.0000000Z',
  'title': '[Act-602/Call-693] Fire in South Korea',
  'wkt': 'POINT(128.50810625 38.24185)'}]

[5]:

call_areas = GeoDataFrame(ciop.search(end_point=activation[0]['link:results'],
                        params=[],
                        output_fields='identifier,self,wkt,link:related',
                        model='GeoTime'))

There are several call areas: a point and two polygons.

[6]:

call_areas

[6]:
identifier link:related self wkt
0 Call-693-Area_0 https://catalog.terradue.com//cos2/series/call... POINT(128.50810625 38.24185)
1 Call-693-Area_1 https://catalog.terradue.com//cos2/series/call... POLYGON((128.418 38.2536,128.4275 38.1604,128....
2 Call-693-Area_2 https://catalog.terradue.com//cos2/series/call... POLYGON((128.49 38.26,128.49 38.2,128.61 38.2,... 
[7]:

call_areas['wkt'] = call_areas['wkt'].apply(loads)

[8]:

call_areas

[8]:
identifier link:related self wkt
0 Call-693-Area_0 https://catalog.terradue.com//cos2/series/call... POINT (128.50810625 38.24185)
1 Call-693-Area_1 https://catalog.terradue.com//cos2/series/call... POLYGON ((128.418 38.2536, 128.4275 38.1604, 1...
2 Call-693-Area_2 https://catalog.terradue.com//cos2/series/call... POLYGON ((128.49 38.26, 128.49 38.2, 128.61 38...

Create a map to visualize the areas of interest

[9]:

for index, row in call_areas.iterrows():

    if row['wkt'].geom_type == 'Point':
        m = Map(center=(row['wkt'].y,
                row['wkt'].x),
                zoom=10)

        break

marker = Marker(location=(row['wkt'].y,
                          row['wkt'].x), draggable=False)

d = {'title': activation[0]['title'],
     'date' : activation[0]['startdate']}

html_value="""
        <div>
        <ul class='list-group'>
           {title}
           <br>
           {date}
        </ul></div>""".format(**d)


html_widget_slave = HTML(
            value=html_value,
    placeholder='',
    description='',
    )


marker.popup = html_widget_slave



m.add_layer(marker);

m

[10]:

call_areas['wkt'].values

[10]:

array([<shapely.geometry.point.Point object at 0x7fef731fc690>,
       <shapely.geometry.polygon.Polygon object at 0x7fef731fc6d0>,
       <shapely.geometry.polygon.Polygon object at 0x7fef731fc590>], dtype=object)

[11]:

aoi = cascaded_union(call_areas['wkt'].values)

[12]:

aoi.wkt

[12]:

'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))'

[13]:

min_lon, min_lat, max_lon, max_lat =  list(aoi.bounds)

Add the area of interest to the map

[14]:

p = Polygon(locations=np.asarray([t[::-1] for t in list(aoi.exterior.coords)]).tolist(), color="red", fill_color="green")

layer_group = LayerGroup(layers=(m.layers), name='Call AOI')

layer_group.add_layer(p)

m.add_layer(layer_group)

Search for Sentinel-2 acquisitions

Search for Sentinle-2 acquistions over the call AOI and over a period of two weeks centered on the activation call:

[15]:

sentinel2_endpoint = 'https://catalog.terradue.com/sentinel2/description'

[16]:

search_params = dict()

search_params['geom'] = aoi.wkt
search_params['pt'] = 'S2MSI2A'
search_params['start'] = (pd.to_datetime(activation[0]['startdate']) - np.timedelta64(7, 'D')).isoformat() + 'Z'
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'
search_params['do'] = 'terradue'
search_params['cc'] = '20['

[17]:

search_params

[17]:

{'cc': '20[',
 'do': 'terradue',
 '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))',
 'pt': 'S2MSI2A',
 'start': '2019-03-29T00:49:00Z',
 'stop': '2019-04-12T00:49:00Z'}

[18]:

sentinel2_search = GeoDataFrame(ciop.search(end_point=sentinel2_endpoint,
                        params=search_params,
                        output_fields='identifier,self,wkt,startdate,enddate,enclosure,orbitDirection,cc',
                        model='EOP'))

[19]:

sentinel2_search

[19]:
cc enclosure enddate identifier orbitDirection self startdate wkt
0 14.461197 https://store.terradue.com/download/sentinel2/... 2019-04-08T02:16:09.0240000Z S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-08T02:16:09.0240000Z POLYGON((127.862827804024 37.8539506483739,129...
1 16.589038 https://store.terradue.com/download/sentinel2/... 2019-04-05T02:06:59.0240000Z S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-05T02:06:59.0240000Z POLYGON((128.234033401221 37.855568274556,129....
2 4.844317 https://store.terradue.com/download/sentinel2/... 2019-04-03T02:16:51.0240000Z S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-03T02:16:51.0240000Z POLYGON((127.862827804024 37.8539506483739,129... 
[20]:

sentinel2_search['wkt'] = sentinel2_search['wkt'].apply(loads)

[21]:

sentinel2_search['identifier'].values

[21]:

array(['S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111',
       'S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128',
       'S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016'], dtype=object)

[22]:

sentinel2_search

[22]:
cc enclosure enddate identifier orbitDirection self startdate wkt
0 14.461197 https://store.terradue.com/download/sentinel2/... 2019-04-08T02:16:09.0240000Z S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-08T02:16:09.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1...
1 16.589038 https://store.terradue.com/download/sentinel2/... 2019-04-05T02:06:59.0240000Z S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-05T02:06:59.0240000Z POLYGON ((128.234033401221 37.855568274556, 12...
2 4.844317 https://store.terradue.com/download/sentinel2/... 2019-04-03T02:16:51.0240000Z S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-03T02:16:51.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1... 
[23]:

def analyse(row, aoi):

    aoi_intersection = (aoi.intersection(row['wkt']).area / aoi.area) * 100

    series = dict([('aoi_intersection', aoi_intersection)])

    series['utm_zone'] = row['identifier'][39:41]
    series['latitude_band'] = row['identifier'][41]
    series['grid_square']  = row['identifier'][42:44]


    return pd.Series(series)

[24]:

sentinel2_search = sentinel2_search.merge(sentinel2_search.apply(lambda row: analyse(row, aoi), axis=1),
              left_index=True,
              right_index=True)

[25]:

sentinel2_search

[25]:
cc enclosure enddate identifier orbitDirection self startdate wkt aoi_intersection grid_square latitude_band utm_zone
0 14.461197 https://store.terradue.com/download/sentinel2/... 2019-04-08T02:16:09.0240000Z S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-08T02:16:09.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 DH S 52
1 16.589038 https://store.terradue.com/download/sentinel2/... 2019-04-05T02:06:59.0240000Z S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-05T02:06:59.0240000Z POLYGON ((128.234033401221 37.855568274556, 12... 100.0 DH S 52
2 4.844317 https://store.terradue.com/download/sentinel2/... 2019-04-03T02:16:51.0240000Z S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-03T02:16:51.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 DH S 52

Stage in the selected Sentinel-2 data

[26]:

target_dir = '/workspace/data/'

def stage_in(row):


    try:
        local_path = ciop.copy(row['enclosure'], extract=True, target=target_dir)

    except:
        local_path = os.path.join(target_dir, row['identifier'] )

    row['local_path'] = os.path.join(target_dir, row['identifier'])

    return row

[27]:

sentinel2_search = sentinel2_search.apply(lambda row: stage_in(row), axis=1)

[28]:

sentinel2_search

[28]:
cc enclosure enddate identifier orbitDirection self startdate wkt aoi_intersection grid_square latitude_band utm_zone local_path
0 14.461197 https://store.terradue.com/download/sentinel2/... 2019-04-08T02:16:09.0240000Z S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-08T02:16:09.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 DH S 52 /workspace/data/S2B_MSIL2A_20190408T021609_N02...
1 16.589038 https://store.terradue.com/download/sentinel2/... 2019-04-05T02:06:59.0240000Z S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-05T02:06:59.0240000Z POLYGON ((128.234033401221 37.855568274556, 12... 100.0 DH S 52 /workspace/data/S2B_MSIL2A_20190405T020659_N02...
2 4.844317 https://store.terradue.com/download/sentinel2/... 2019-04-03T02:16:51.0240000Z S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_2... DESCENDING https://catalog.terradue.com/sentinel2/search?... 2019-04-03T02:16:51.0240000Z POLYGON ((127.862827804024 37.8539506483739, 1... 100.0 DH S 52 /workspace/data/S2A_MSIL2A_20190403T021651_N02...

Near Infrared and Short-Wave Infrared RGB Composite

Using the Sentinel-2 bands B12, B11 and B8A, generate the RGB composite

[29]:

bands = ['B12', 'B11', 'B8A']

[30]:

def get_band_path(row, band):

    ns = {'xfdu': 'urn:ccsds:schema:xfdu:1',
          'safe': 'http://www.esa.int/safe/sentinel/1.1',
          'gml': 'http://www.opengis.net/gml'}

    path_manifest = os.path.join(row['local_path'],
                                 row['identifier'] + '.SAFE',
                                'manifest.safe')

    root = etree.parse(path_manifest)

    bands = [band]

    for index, band in enumerate(bands):

        sub_path = os.path.join(row['local_path'],
                                row['identifier'] + '.SAFE',
                                root.xpath('//dataObjectSection/dataObject/byteStream/fileLocation[contains(@href,("%s%s")) and contains(@href,("%s")) ]' % (row['latitude_band'],
                                row['grid_square'],
                                band),
                                  namespaces=ns)[0].attrib['href'][2:])

    return sub_path

Use OTB for the contrast enhancement:

[31]:

def contrast_enhancement(in_tif, out_tif, hfact=2.0):

    ContrastEnhancement = otbApplication.Registry.CreateApplication("ContrastEnhancement")

    ContrastEnhancement.SetParameterString("in", in_tif)
    ContrastEnhancement.SetParameterString("out", out_tif)
    ContrastEnhancement.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_uint8)
    ContrastEnhancement.SetParameterFloat("nodata", 0.0)
    ContrastEnhancement.SetParameterFloat("hfact", hfact)
    ContrastEnhancement.SetParameterInt("bins", 256)
    ContrastEnhancement.SetParameterInt("spatial.local.w", 500)
    ContrastEnhancement.SetParameterInt("spatial.local.h", 500)
    ContrastEnhancement.SetParameterString("mode","lum")

    ContrastEnhancement.ExecuteAndWriteOutput()


[32]:

rgb_composites = []

for index, row in sentinel2_search.iterrows():

    vrt_bands = []

    for j, band in enumerate(bands):

        vrt_bands.append(get_band_path(row, band))

    vrt = '{0}.vrt'.format(row['identifier'])
    ds = gdal.BuildVRT(vrt,
                       vrt_bands,
                       srcNodata=0,
                       xRes=10,
                       yRes=10,
                       separate=True)
    ds.FlushCache()

    tif =  '{0}.tif'.format(row['identifier'])
    gdal.Translate(tif,
                       vrt,
                       projWin=[min_lon, max_lat, max_lon, min_lat],
                       projWinSRS='EPSG:4326',
                       outputType=gdal.GDT_Byte,
                       scaleParams=[[0,10000,0,255]])

    tif_e =  '{0}_e.tif'.format(row['identifier'])

    contrast_enhancement(tif, tif_e)

    rgb_composites.append(tif_e)

Create a new map with the RGB composites:

[33]:

m = Map(center=(aoi.centroid.y,
                aoi.centroid.x),
                zoom=13)

m.add_control(LayersControl())

[34]:

rgb_pngs = []

layers = []

for index, composite in enumerate(rgb_composites):

    im = Image.open(composite)

    f = StringIO()

    im.save(f, 'png')
    data = b64encode(f.getvalue())

    rgb_pngs.append('data:image/png;base64,' + data)

    image = ImageOverlay(
        url=rgb_pngs[index],
        bounds=((min_lat, min_lon), (max_lat, max_lon))
        )

    layer_group = LayerGroup(layers=(), name='RGB composite #{0}'.format(index))

    layer_group.add_layer(image)

    m.add_layer(layer_group)

    layers.append(layer_group)

[35]:

m

Pre-processing for the burned area analysis

Use a SNAP graph to calculate the NDWI (to descriminate the water) and the NBR:

[36]:

def pre_processing(**kwargs):

    options = dict()

    operators = ['Read',
                 'Resample',
                 'BandMaths',
                 'Reproject',
                 'Subset',
                 'Write']

    for operator in operators:

        print 'Getting default values for Operator {}'.format(operator)
        parameters = get_operator_default_parameters(operator)

        options[operator] = parameters

    for key, value in kwargs.items():

        print 'Updating Operator {}'.format(key)
        options[key.replace('_', '-')].update(value)

    mygraph = GraphProcessor()

    for index, operator in enumerate(operators):

        print 'Adding Operator {} to graph'.format(operator)
        if index == 0:
            source_node_id = ''

        else:
            source_node_id = operators[index - 1]

        mygraph.add_node(operator,
                         operator,
                         options[operator], source_node_id)

    mygraph.view_graph()

    mygraph.run()

[37]:

pre_processed = []

resample = dict()
resample['referenceBandName'] = 'B2'

reproject = dict()
reproject['crs'] = 'EPSG:4326'

subset = dict()
subset['geoRegion'] = box(*aoi.bounds).wkt
subset['copyMetadata'] = 'true'

bands = '''<targetBands>
    <targetBand>
      <name>NDWI</name>
      <type>float32</type>
      <expression>(B3 - B8) / (B3 + B8)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>NBR</name>
      <type>float32</type>
      <expression>(B8 - B12) / (B8 + B12)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>'''

band_maths = dict()
band_maths['targetBandDescriptors'] = bands

for index, row in sentinel2_search.iterrows():

    print os.path.join(row['local_path'], row['identifier'] + '.SAFE', 'MTD_MSIL2A.xml')

    read = dict()
    read['file'] = os.path.join(row['local_path'], row['identifier'] + '.SAFE', 'MTD_MSIL2A.xml') #, 'manifest.safe')
    read['formatName'] = 'SENTINEL-2-MSI-MultiRes-UTM52N'

    write = dict()
    write['file'] = 'pre_{}'.format(row['identifier'])

    pre_processed.append('pre_{}'.format(row['identifier']))

    pre_processing(Read=read,
                 Resample=resample,
                 Reproject=reproject,
                 Subset=subset,
                 BandMaths=band_maths,
                 Write=write)

/workspace/data/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.SAFE/MTD_MSIL2A.xml
Getting default values for Operator Read
Getting default values for Operator Resample
Getting default values for Operator BandMaths
Getting default values for Operator Reproject
Getting default values for Operator Subset
Getting default values for Operator Write
Updating Operator Subset
Updating Operator Read
Updating Operator Resample
Updating Operator BandMaths
Updating Operator Write
Updating Operator Reproject
Adding Operator Read to graph
Adding Operator Resample to graph
Adding Operator BandMaths to graph
Adding Operator Reproject to graph
Adding Operator Subset to graph
Adding Operator Write to graph
<graph>
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName>SENTINEL-2-MSI-MultiRes-UTM52N</formatName>
      <file>/workspace/data/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111/S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.SAFE/MTD_MSIL2A.xml</file>
    </parameters>
  </node>
  <node id="Resample">
    <operator>Resample</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetHeight/>
      <flagDownsamplingMethod>First</flagDownsamplingMethod>
      <targetWidth/>
      <resampleOnPyramidLevels>true</resampleOnPyramidLevels>
      <downsamplingMethod>First</downsamplingMethod>
      <upsamplingMethod>Nearest</upsamplingMethod>
      <targetResolution/>
      <referenceBandName>B2</referenceBandName>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Resample"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
    <targetBand>
      <name>NDWI</name>
      <type>float32</type>
      <expression>(B3 - B8) / (B3 + B8)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>NBR</name>
      <type>float32</type>
      <expression>(B8 - B12) / (B8 + B12)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Reproject">
    <operator>Reproject</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <crs>EPSG:4326</crs>
      <orthorectify>false</orthorectify>
      <pixelSizeX/>
      <orientation>0</orientation>
      <elevationModelName/>
      <includeTiePointGrids>true</includeTiePointGrids>
      <pixelSizeY/>
      <resamplingName>Nearest</resamplingName>
      <height/>
      <width/>
      <wktFile/>
      <easting/>
      <referencePixelX/>
      <referencePixelY/>
      <addDeltaBands>false</addDeltaBands>
      <noDataValue/>
      <tileSizeX/>
      <tileSizeY/>
      <northing/>
    </parameters>
  </node>
  <node id="Subset">
    <operator>Subset</operator>
    <sources>
      <sourceProduct refid="Reproject"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <tiePointGridNames/>
      <bandNames/>
      <subSamplingY>1</subSamplingY>
      <subSamplingX>1</subSamplingX>
      <region/>
      <copyMetadata>true</copyMetadata>
      <fullSwath>false</fullSwath>
      <geoRegion>POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))</geoRegion>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="Subset"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <clearCacheAfterRowWrite>false</clearCacheAfterRowWrite>
      <formatName>BEAM-DIMAP</formatName>
      <file>pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111</file>
      <deleteOutputOnFailure>true</deleteOutputOnFailure>
      <writeEntireTileRows>true</writeEntireTileRows>
    </parameters>
  </node>
</graph>

Processing the graph
Process PID: 3584
('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')
Done.
/workspace/data/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128.SAFE/MTD_MSIL2A.xml
Getting default values for Operator Read
Getting default values for Operator Resample
Getting default values for Operator BandMaths
Getting default values for Operator Reproject
Getting default values for Operator Subset
Getting default values for Operator Write
Updating Operator Subset
Updating Operator Read
Updating Operator Resample
Updating Operator BandMaths
Updating Operator Write
Updating Operator Reproject
Adding Operator Read to graph
Adding Operator Resample to graph
Adding Operator BandMaths to graph
Adding Operator Reproject to graph
Adding Operator Subset to graph
Adding Operator Write to graph
<graph>
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName>SENTINEL-2-MSI-MultiRes-UTM52N</formatName>
      <file>/workspace/data/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128/S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128.SAFE/MTD_MSIL2A.xml</file>
    </parameters>
  </node>
  <node id="Resample">
    <operator>Resample</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetHeight/>
      <flagDownsamplingMethod>First</flagDownsamplingMethod>
      <targetWidth/>
      <resampleOnPyramidLevels>true</resampleOnPyramidLevels>
      <downsamplingMethod>First</downsamplingMethod>
      <upsamplingMethod>Nearest</upsamplingMethod>
      <targetResolution/>
      <referenceBandName>B2</referenceBandName>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Resample"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
    <targetBand>
      <name>NDWI</name>
      <type>float32</type>
      <expression>(B3 - B8) / (B3 + B8)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>NBR</name>
      <type>float32</type>
      <expression>(B8 - B12) / (B8 + B12)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Reproject">
    <operator>Reproject</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <crs>EPSG:4326</crs>
      <orthorectify>false</orthorectify>
      <pixelSizeX/>
      <orientation>0</orientation>
      <elevationModelName/>
      <includeTiePointGrids>true</includeTiePointGrids>
      <pixelSizeY/>
      <resamplingName>Nearest</resamplingName>
      <height/>
      <width/>
      <wktFile/>
      <easting/>
      <referencePixelX/>
      <referencePixelY/>
      <addDeltaBands>false</addDeltaBands>
      <noDataValue/>
      <tileSizeX/>
      <tileSizeY/>
      <northing/>
    </parameters>
  </node>
  <node id="Subset">
    <operator>Subset</operator>
    <sources>
      <sourceProduct refid="Reproject"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <tiePointGridNames/>
      <bandNames/>
      <subSamplingY>1</subSamplingY>
      <subSamplingX>1</subSamplingX>
      <region/>
      <copyMetadata>true</copyMetadata>
      <fullSwath>false</fullSwath>
      <geoRegion>POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))</geoRegion>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="Subset"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <clearCacheAfterRowWrite>false</clearCacheAfterRowWrite>
      <formatName>BEAM-DIMAP</formatName>
      <file>pre_S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128</file>
      <deleteOutputOnFailure>true</deleteOutputOnFailure>
      <writeEntireTileRows>true</writeEntireTileRows>
    </parameters>
  </node>
</graph>

Processing the graph
Process PID: 3673
('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')
Done.
/workspace/data/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.SAFE/MTD_MSIL2A.xml
Getting default values for Operator Read
Getting default values for Operator Resample
Getting default values for Operator BandMaths
Getting default values for Operator Reproject
Getting default values for Operator Subset
Getting default values for Operator Write
Updating Operator Subset
Updating Operator Read
Updating Operator Resample
Updating Operator BandMaths
Updating Operator Write
Updating Operator Reproject
Adding Operator Read to graph
Adding Operator Resample to graph
Adding Operator BandMaths to graph
Adding Operator Reproject to graph
Adding Operator Subset to graph
Adding Operator Write to graph
<graph>
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName>SENTINEL-2-MSI-MultiRes-UTM52N</formatName>
      <file>/workspace/data/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016/S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.SAFE/MTD_MSIL2A.xml</file>
    </parameters>
  </node>
  <node id="Resample">
    <operator>Resample</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetHeight/>
      <flagDownsamplingMethod>First</flagDownsamplingMethod>
      <targetWidth/>
      <resampleOnPyramidLevels>true</resampleOnPyramidLevels>
      <downsamplingMethod>First</downsamplingMethod>
      <upsamplingMethod>Nearest</upsamplingMethod>
      <targetResolution/>
      <referenceBandName>B2</referenceBandName>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Resample"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
    <targetBand>
      <name>NDWI</name>
      <type>float32</type>
      <expression>(B3 - B8) / (B3 + B8)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>NBR</name>
      <type>float32</type>
      <expression>(B8 - B12) / (B8 + B12)</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Reproject">
    <operator>Reproject</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <crs>EPSG:4326</crs>
      <orthorectify>false</orthorectify>
      <pixelSizeX/>
      <orientation>0</orientation>
      <elevationModelName/>
      <includeTiePointGrids>true</includeTiePointGrids>
      <pixelSizeY/>
      <resamplingName>Nearest</resamplingName>
      <height/>
      <width/>
      <wktFile/>
      <easting/>
      <referencePixelX/>
      <referencePixelY/>
      <addDeltaBands>false</addDeltaBands>
      <noDataValue/>
      <tileSizeX/>
      <tileSizeY/>
      <northing/>
    </parameters>
  </node>
  <node id="Subset">
    <operator>Subset</operator>
    <sources>
      <sourceProduct refid="Reproject"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <tiePointGridNames/>
      <bandNames/>
      <subSamplingY>1</subSamplingY>
      <subSamplingX>1</subSamplingX>
      <region/>
      <copyMetadata>true</copyMetadata>
      <fullSwath>false</fullSwath>
      <geoRegion>POLYGON ((128.6102 38.1604, 128.6102 38.3037, 128.4164 38.3037, 128.4164 38.1604, 128.6102 38.1604))</geoRegion>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="Subset"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <clearCacheAfterRowWrite>false</clearCacheAfterRowWrite>
      <formatName>BEAM-DIMAP</formatName>
      <file>pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016</file>
      <deleteOutputOnFailure>true</deleteOutputOnFailure>
      <writeEntireTileRows>true</writeEntireTileRows>
    </parameters>
  </node>
</graph>

Processing the graph
Process PID: 3762
('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')
Done.

Collocation

Collocate the Sentinel-2 processed results using as Master the pre-fire and as Slave the post-fire Sentinel-2 acquisitions

[38]:

pre_processed

[38]:

['pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111',
 'pre_S2B_MSIL2A_20190405T020659_N0211_R103_T52SDH_20190405T043128',
 'pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016']

[56]:

master_path = pre_processed[2] + '.dim'
slave_path = pre_processed[0] + '.dim'

[57]:

mygraph = GraphProcessor()

[58]:

operator = 'Read'

node_id = 'Read_M'

source_node_id = ''

parameters = get_operator_default_parameters(operator)

parameters['file'] = master_path

mygraph.add_node(node_id,
                 operator,
                         parameters, source_node_id)


[59]:

operator = 'Read'

node_id = 'Read_S'

source_node_id = ''

parameters = get_operator_default_parameters(operator)

parameters['file'] = slave_path

mygraph.add_node(node_id,
                 operator,
                 parameters, source_node_id)

[60]:

operator = 'Collocate'

parameters = get_operator_default_parameters(operator)

parameters['masterComponentPattern'] = 'PRE_FIRE_${ORIGINAL_NAME}'
parameters['slaveComponentPattern'] = 'POST_FIRE_${ORIGINAL_NAME}'

source_node_id = dict()

source_node_id['master'] = 'Read_M'

source_node_id['slave'] = 'Read_S'


node_id = 'Collocate'

mygraph.add_node(operator,
                         operator,
                         parameters, source_node_id)

[61]:


operator = 'Write'

node_id = 'Write'

source_node_id = 'Collocate'



parameters = get_operator_default_parameters(operator)

parameters['file'] = 'collocated'
parameters['formatName'] = 'BEAM-DIMAP'

mygraph.add_node(                         node_id,
                 operator,
                         parameters, source_node_id)

[62]:

mygraph.view_graph()

<graph>
  <version>1.0</version>
  <node id="Read_M">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName/>
      <file>pre_S2A_MSIL2A_20190403T021651_N0211_R003_T52SDH_20190404T105016.dim</file>
    </parameters>
  </node>
  <node id="Read_S">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName/>
      <file>pre_S2B_MSIL2A_20190408T021609_N0211_R003_T52SDH_20190408T045111.dim</file>
    </parameters>
  </node>
  <node id="Collocate">
    <operator>Collocate</operator>
    <sources>
      <slave>Read_S</slave>
      <master>Read_M</master>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <renameSlaveComponents>true</renameSlaveComponents>
      <slaveComponentPattern>POST_FIRE_${ORIGINAL_NAME}</slaveComponentPattern>
      <targetProductName>_collocated</targetProductName>
      <renameMasterComponents>true</renameMasterComponents>
      <targetProductType>COLLOCATED</targetProductType>
      <masterComponentPattern>PRE_FIRE_${ORIGINAL_NAME}</masterComponentPattern>
      <resamplingType>NEAREST_NEIGHBOUR</resamplingType>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="Collocate"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <deleteOutputOnFailure>true</deleteOutputOnFailure>
      <clearCacheAfterRowWrite>false</clearCacheAfterRowWrite>
      <formatName>BEAM-DIMAP</formatName>
      <writeEntireTileRows>true</writeEntireTileRows>
      <file>collocated</file>
    </parameters>
  </node>
</graph>


[63]:

mygraph.run()

Processing the graph
Process PID: 4000
('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')
Done.

dNBR, RBR and burned area severity

[64]:

def burned_area(**kwargs):

    options = dict()

    operators = ['Read',
                 'BandMaths',
                 'Write']

    for operator in operators:

        print 'Getting default values for Operator {}'.format(operator)
        parameters = get_operator_default_parameters(operator)

        options[operator] = parameters

    for key, value in kwargs.items():

        print 'Updating Operator {}'.format(key)
        options[key.replace('_', '-')].update(value)

    mygraph = GraphProcessor()

    for index, operator in enumerate(operators):

        print 'Adding Operator {} to graph'.format(operator)
        if index == 0:
            source_node_id = ''

        else:
            source_node_id = operators[index - 1]

        mygraph.add_node(operator,
                         operator,
                         options[operator], source_node_id)

    mygraph.view_graph()

    mygraph.run()

[65]:

collocated_input = 'collocated.dim'

read = dict()
read['file'] = collocated_input

bands = '''<targetBands>
    <targetBand>
      <name>dNBR</name>
      <type>float32</type>
      <expression>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 </expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>RBR</name>
      <type>float32</type>
      <expression>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</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>'''

band_maths = dict()
band_maths['targetBandDescriptors'] = bands

write = dict()
write['file'] = 'burned_area'
write['formatName'] = 'GeoTIFF'

burned_area(Read=read,
            BandMaths=band_maths,
            Write=write)

Getting default values for Operator Read
Getting default values for Operator BandMaths
Getting default values for Operator Write
Updating Operator Read
Updating Operator Write
Updating Operator BandMaths
Adding Operator Read to graph
Adding Operator BandMaths to graph
Adding Operator Write to graph
<graph>
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <formatName/>
      <file>collocated.dim</file>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
    <targetBand>
      <name>dNBR</name>
      <type>float32</type>
      <expression>PRE_FIRE_NDWI &gt;= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) &gt; 0.27 ? PRE_FIRE_NBR - POST_FIRE_NBR : 0 </expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    <targetBand>
      <name>RBR</name>
      <type>float32</type>
      <expression>PRE_FIRE_NDWI &gt;= 0.0 ? 0 : ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) &gt; 0.27 ? ((PRE_FIRE_NBR - POST_FIRE_NBR) / (PRE_FIRE_NBR + 1.001)) : 0</expression>
      <description/>
      <unit/>
      <noDataValue>NaN</noDataValue>
    </targetBand>
    </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <clearCacheAfterRowWrite>false</clearCacheAfterRowWrite>
      <formatName>GeoTIFF</formatName>
      <file>burned_area</file>
      <deleteOutputOnFailure>true</deleteOutputOnFailure>
      <writeEntireTileRows>true</writeEntireTileRows>
    </parameters>
  </node>
</graph>

Processing the graph
Process PID: 4083
('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')
Done.

Burned area severity

[66]:

def raster2rgb(raster_file, color_table, out_file_name, raster_band=1, discrete=True):

    #Reading the band
    data_types ={'Byte':'B','UInt16':'H','Int16':'h','UInt32':'I','Int32':'i','Float32':'f','Float64':'d'}
    if exists(raster_file) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + raster_file + '\'')

    dataset = gdal.Open(raster_file, GA_ReadOnly )

    if dataset == None:
        raise Exception("Unable to read the data file")

    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()

    band = dataset.GetRasterBand(raster_band)
    values = band.ReadRaster( 0, 0, band.XSize, band.YSize, band.XSize, band.YSize, band.DataType )
    values = unpack(data_types[gdal.GetDataTypeName(band.DataType)]*band.XSize*band.YSize,values)

    #Preparing the color table and the output file
    classification_values = color_table.keys()
    classification_values.sort()

    base = Image.new( 'RGBA', (band.XSize,band.YSize) )
    base_draw = ImageDraw.Draw(base)
    alpha_mask = Image.new('L', (band.XSize,band.YSize), 255)
    alpha_draw = ImageDraw.Draw(alpha_mask)

    #Reading the value and setting the output color for each pixel
    for pos in range(len(values)):
        y = pos/band.XSize
        x = pos - y * band.XSize
        for index in range(len(classification_values)):

            if values[pos] <= classification_values[index] or index == len(classification_values)-1:
                if discrete == True:
                    if index == 0:
                        index = 1
                    elif index == len(classification_values)-1 and values[pos] >= classification_values[index]:
                        index = index + 1
                    color = color_table[classification_values[index-1]]
                    base_draw.point((x,y), (color[0],color[1],color[2]))
                    alpha_draw.point((x,y),color[3])
                else:
                    if index == 0:
                        r = color_table[classification_values[0]][0]
                        g = color_table[classification_values[0]][1]
                        b = color_table[classification_values[0]][2]
                        a = color_table[classification_values[0]][3]
                    elif index == len(classification_values)-1 and values[pos] >= classification_values[index]:
                        r = color_table[classification_values[index]][0]
                        g = color_table[classification_values[index]][1]
                        b = color_table[classification_values[index]][2]
                        a = color_table[classification_values[index]][3]
                    else:
                        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])
                        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])
                        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])
                        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])

                    base_draw.point((x,y), (int(r),int(g),int(b)))
                    alpha_draw.point((x,y),int(a))

                break
    #Adding transparency and saving the output image
    color_layer = Image.new('RGBA', base.size, (255, 255, 255, 0))
    base = Image.composite(color_layer, base, alpha_mask)
    base.save(out_file_name)

    # update geolocation
    ds_rgb = gdal.Open(out_file_name,1)
    ds_rgb.SetGeoTransform(geoTransform)
    ds_rgb.SetProjection(proj)

    ds_rgb.FlushCache()

    ds_rgb = None

[67]:

severity_palette = {-1: [159, 159, 159, 0], # grey
                    -0.5: [43, 25, 223, 0], # blue
                    -0.251: [139, 221, 231, 0], # cyan
                    -0.101: [97, 169, 45, 255], # unburned, green
                    0.099: [250, 254, 76, 0], # yellow
                    0.269: [228, 173, 55, 0], # orange
                    0.439: [202, 59, 18, 0],  # red
                    0.659: [85, 15, 112, 0]} # purple

[68]:

raster2rgb('burned_area.tif',
           severity_palette,
           'burned_area.rgb.tif',
           raster_band=1,
           discrete=True)

View the results on a map

[69]:

im = Image.open('burned_area.rgb.tif')

f = StringIO()

im.save(f, 'png')
data = b64encode(f.getvalue())

burned_area = 'data:image/png;base64,' + data

[70]:

m = Map(center=(aoi.centroid.y,
                aoi.centroid.x),
                zoom=13)

m.add_control(LayersControl())

[71]:

image = ImageOverlay(
    url=burned_area,
    bounds=((min_lat, min_lon), (max_lat, max_lon))
)

layer_group = LayerGroup(layers=(m.layers), name='Burned area severity')

layer_group.add_layer(image)

m.add_layer(layer_group)

# add the RGB layers
for layer in layers:

    m.add_layer(layer)




[72]:

m

This work is licenced under a Attribution-ShareAlike 4.0 International License (CC BY-SA 4.0)

YOU ARE FREE TO:

  • Share - copy and redistribute the material in any medium or format.
  • Adapt - remix, transform, and built upon the material for any purpose, even commercially.

UNDER THE FOLLOWING TERMS:

  • 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.
  • ShareAlike - If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.
[ ]:



[ ]: