Multi-Choice Dashboard

Note: This notebook is a dashboard, check out this link to preview: https://nbviewer.jupyter.org/github/GeoscienceAustralia/dea-notebooks/blob/master/06_Composite_generation/MultiChoice_Dashboard.ipynb

What does this notebook do?: This notebook interactively creates a composite image of Landsat 8 data. You can choose to create a geomedian or a tasseled cap index composite. If you choose the tasseled cap composite, you can choose to plot the percentage exceedance, the means, the wetness summary or the standard deviations. You can then save the data as a png image, netcdf and geotiff.

Background: Landsat 8 data is available from April 2013 onwards.

The geometric median (geomedian) is a spatial median composite of satellite data that can be used to look at the statistically “average” behavior of an area summarising the chosen time interval. The geometric median preserves relationships between spectral bands - so instead of getting an “average” behaviour for each band, you get an overall measure of the central tendency that works spatially, temporally and spectrally.

The Tasseled Cap Index (TCI) is a method of reducing 6 bands of satellite data (BLUE, GREEN, RED, NIR, SWIR1, SWIR2) to 3 bands (Brightness, Greenness, Wetness) using a Principal Components Analysis and Procrustes’ Rotation (Roberts et al 2018). This notebook uses the published coefficients of Crist 1985 as applied to Digital Earth Australia’s Landsat satellite data.

Tasseled cap percentage exceedance: Percent of observations within epoch where TCI is above chosen threshold (i.e. where pixel has been ‘wet’ ‘green or ‘bright’).

Before you run this notebook: This notebook uses the dea statistics module. You need to run “module load dea” in a terminal window and then launch jupyter notebooks in the same window so that your notebook can ‘see’ the stats module.

Date: June 2018

Authors: Bex Dunn, Vanessa Newey, Erin Telfer

Tags: Landsat8, geomedian, shapefile, ipywidgets, Geotiff, tasseled_cap, TCWstats, pct_exceedance, NetCDF, polygondrill, DrawnPolygon, SelectFileButton

Import modules from standard libraries, datacube and files

Select ‘Trust this notebook’ to import these modules and run the dashboard.

[3]:
#import standard libraries
import datetime #for widget date validation
import fiona
import numpy as np
import os
import pandas as pd
import rasterio
import sys
import xarray as xr
from ipywidgets import interact #for widgets
from IPython.display import display
import ipywidgets as widgets
import matplotlib.pyplot as plt
from matplotlib.pyplot import imshow
from mpl_toolkits.axes_grid1 import make_axes_locatable
import rasterio.features
import shapely.geometry


#modules for datacube and dea stats
import datacube
from datacube.helpers import ga_pq_fuser, write_geotiff
from datacube.storage import masking
from datacube.storage.masking import mask_to_dict
from datacube.storage.storage import write_dataset_to_netcdf
from datacube_stats.statistics import GeoMedian, TCWStats

from datacube.utils import geometry
# set datacube alias
dc = datacube.Datacube(app='dc-tcw and geomedian')

# Import external functions from dea-notebooks
sys.path.append(os.path.expanduser('~/dea-notebooks/10_Scripts/'))
import DEAPlotting, DEADataHandling
from FileDialogs import *

#ignore datacube warnings (needs to be last import statement)
import warnings
warnings.filterwarnings('ignore', module='datacube')

#set global variables for widgets
global GEOM
global NBART
global NBART_GM
global NBART_TC
global OUTPUT_STAT
global OUTPUT_STAT_NAME
global SHAPE_NAME
GEOM = None
OUTPUT_STAT=None
OUTPUT_STAT_NAME=None
NBART = None
NBART_GM = None
NBART_TC = None
SHAPE_NAME =None

Select a shape file for your area of analysis

  • Shape file can be selected by either choosing a file or by drawing a polygon below

  • Note: if area is larger than 20km x 20km or time period is long, notebook may run out of memory

[ ]:
#code to work with an input polygon
style = {'description_width': 'initial'}
shape_file = '/g/data1a/r78/rjd547/groundwater_activities/Burdekin/Burdekin_shapefiles/fels_sml.shp'
shape_file_text = widgets.Text(value=shape_file,placeholder='choose a file',
    description='path to shape file',
    style = {'description_width': 'initial'},
    disabled=False,
    layout=widgets.Layout(width='70%'))
def handle_shape_file(sender):
    shape_file=shape_file_text.value
shape_file_text.observe(handle_shape_file)
#display(shape_file_text)

Use the select file button to choose a shape file

[ ]:
f = SelectFileButton()
def handle_file_selector(sender):
    shape_file_text.value = f.file
f.on_click(handle_file_selector)
display(widgets.HBox([shape_file_text,f],disabled=False))

Use the check shapefile button to check if the chosen shapefile is valid

[10]:
button =widgets.Button(
    description='Check shapefile',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Click me',
)

def on_button_clicked(b):
    #check whether shapefile is a valid shapefile and exists
    if not os.path.isfile(shape_file_text.value):
        print('this is not a shapefile')
    else:
        try:
            with fiona.open(shape_file_text.value, 'r') as input:
            # check that the file type is a shapefile
                if input.driver == 'ESRI Shapefile':
                    global GEOM
                    global SHAPE_NAME
                    GEOM, SHAPE_NAME = DEADataHandling.open_polygon_from_shapefile(shape_file_text.value)
                    print ("Specified file is a shapefile!")
                else:
                    print ("NOT a shapefile!")
        except:
            print ("NOT a shapefile!")

    #shapes = fiona.open(shape_file)
button.on_click(on_button_clicked)
display(button)

Alternatively, draw a polygon on the map below

After the polygon is drawn, select ‘DrawnPolygon’ below figure to use the polygon.

[11]:
from ipyleaflet import (
    Map,
    Marker,
    TileLayer, ImageOverlay,
    Polyline, Polygon, Rectangle, Circle, CircleMarker,
    GeoJSON,
    DrawControl,
    basemaps
)


from traitlets import link
import ipywidgets as widgets
from datacube.utils import geometry
#center = [-12.543839666237682, 131.08474731445315]
center =[-20.0183,145.3134]
zoom = 10


m = Map(center=center, zoom=zoom,layout=dict(width='400px', height='400px'))

#Let's draw a second map

m2 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery,layout=m.layout)


# Now create the DrawControls and add it to the Maps using add_control.
# We also register handlers for draw events. These will fire when a polygon is drawn created, edited or deleted
# (there are the actions). The geo_json argument is the serialized geometry of the drawn path,
# along with its embedded style.

draw_control = DrawControl(polyline={})
draw_control2 = DrawControl(polygon={'shapeOptions': {'color': '#0000FF'}}, polyline={})
def handle_draw(self, action, geo_json):
    if action == 'created':
        m2.add_layer(GeoJSON(data=draw_control.last_draw))
        draw_control2.last_draw =draw_control.last_draw
    if action == 'deleted':
        while len(m2.layers)>1:
            m2.remove_layer(m2.layers[1])

def handle_draw2(self, action, geo_json):

    if action == 'created':
        m.add_layer(GeoJSON(data=draw_control2.last_draw))
        draw_control.last_draw =draw_control2.last_draw
    if action == 'deleted':
        while len(m.layers)>1:
            m.remove_layer(m.layers[1])

#add handlers to draw controls
draw_control.on_draw(handle_draw)
draw_control2.on_draw(handle_draw2)

#add draw controls to maps
m.add_control(draw_control)
m2.add_control(draw_control2)

#We can use link to synchronize traitlets of the two maps:

map_center_link = link((m, 'center'), (m2, 'center'))
map_zoom_link = link((m, 'zoom'), (m2, 'zoom'))

#use_map_geom_button = widgets.Button(description="Use Drawn Geom")

use_map_geom_button = widgets.ToggleButtons(
    options=['Shapefile', 'DrawnPolygon'],
    description='Geometry Used for Datacube Query:',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltips=['Use Shapefile Specified for the DC query', 'Use Polygon Drawn on the Map for the DC query'],
#     icons=['check'] * 3
)

def on_map_button_click(b):
    global GEOM
    global SHAPE_NAME
    #Leaflet return geometry in EPSG:4326
    geom_crs = geometry.CRS('EPSG:4326')
    if use_map_geom_button.value == "Shapefile":
        GEOM, SHAPE_NAME = DEADataHandling.open_polygon_from_shapefile(shape_file_text.value)
    else:
        try:
            GEOM=geometry.Geometry(draw_control.last_draw['geometry'], geom_crs)
            SHAPE_NAME = "drawn_polygon"
        except:
            print("Draw a polygon first.")

use_map_geom_button.on_trait_change(on_map_button_click)


from ipywidgets import HBox, VBox
h_box= HBox([m,m2])

VBox([h_box,use_map_geom_button])
/g/data/v10/public/modules/dea-env/20180405/lib/python3.6/site-packages/ipykernel_launcher.py:89: DeprecationWarning: on_trait_change is deprecated in traitlets 4.1: use observe instead

Fill in boxes to select time period of interest

Press the ‘Check dates’ button to check that the time period is valid

[12]:
#Set a standard time range for people to change as desired
start_of_epoch = '2013-05-01'
end_of_epoch =  '2013-10-31'

#create a date picker widget for the start date
from_date_picker = widgets.Text(value=start_of_epoch,
                                placeholder='update this field YYYY-MM-DD',
    description='start date',
    disabled=False)

#define a function to check the start date
def handle_from_date(sender):
    '''takes the start of epoch date from the date picker and checks it is in the
    correct format'''
    start_of_epoch=from_date_picker.value
    validate(start_of_epoch)

#when you press enter this runs the event handler (handle_from_date)
from_date_picker.on_submit(handle_from_date)

#create a date picker widget for the end date
to_date_picker = widgets.Text(value = end_of_epoch,
                              placeholder='update this field YYYY-MM-DD',
    description='end date',
    disabled=False)

def handle_to_date(sender):
    '''takes the end of epoch date from the date picker and checks it is in the
    correct format'''
    end_of_epoch = to_date_picker.value
    validate(end_of_epoch)

#when you press enter this runs the event handler (handle_to_date)
to_date_picker.on_submit(handle_to_date)

#make a validation button to validate the dates
#first write the function the button runs
def validate(date_text):
    '''this function checks to see whether your text is year-month-day formatted'''
    try:
        datetime.datetime.strptime(date_text, '%Y-%m-%d')
        return True
    except ValueError:
        return False

#make the button run the function to check the dates
button_epoch =widgets.Button(
    description='Check dates',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Click me',
)

#This widget prints whether or not your start date is valid
starttextwidget = widgets.Label(
value='Is my start date valid?'
)

#This widget prints whether or not your end date is valid
endtextwidget = widgets.Label(
value='Is my end date valid?'
)

#this is the event handler for when you click on the button
def on_button_end_epoch_clicked(button):
    '''check whether the dates are in the format we need'''
    if not validate(to_date_picker.value):
        endtextwidget.value = ('Incorrect data format, should be YYYY-MM-DD')
    else:
        validate(to_date_picker.value)
        endtextwidget.value = ('End date is valid')


    if not validate(from_date_picker.value):
        starttextwidget.value = ('Incorrect data format, should be YYYY-MM-DD')
    else:
        validate(from_date_picker.value)
        starttextwidget.value = ('Start date is valid')
        start_of_epoch = from_date_picker.value


button_epoch.on_click(on_button_end_epoch_clicked)

#display widgets next to each other
end_hbox = widgets.HBox([to_date_picker, endtextwidget])
start_hbox = widgets.HBox([from_date_picker, starttextwidget])
display(start_hbox,end_hbox)
display(button_epoch)

Load the data

Push the ‘Load data’ button. This may take a few minutes, please be patient.

[13]:
loadButton = widgets.Button(description="Load")

def on_load_button_clicked(b):
    start_of_epoch = from_date_picker.value
    end_of_epoch = to_date_picker.value
    global GEOM
    global SHAPE_NAME
    global NBART
    if GEOM is None:
        GEOM, SHAPE_NAME = DEADataHandling.open_polygon_from_shapefile(shape_file_text.value)

    query = {
        'time': (start_of_epoch, end_of_epoch), 'geopolygon': GEOM,
    }

    sensor3_nbart, crs, affine=DEADataHandling.load_nbarx(dc,sensor='ls8', query=query)

    NBART=sensor3_nbart.sortby('time')

    print("Data loaded")

loadButton.on_click(on_load_button_clicked)
display(loadButton)
Loading ls8_nbart_albers
Loaded ls8_nbart_albers
Generating mask ls8_pq_albers
Data loaded
[19]:
#use_map_geom_button = widgets.Button(description="Use Drawn Geom")

calculate_button = widgets.ToggleButtons(
    options=['Geomedian', 'Tasselled Cap'],
    description='Click selection to calculate stats:',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltips=['Calculate Geomedian', 'Calculate Tasselled Cap Stats'],
#     icons=['check'] * 3
)

@interact
def clicked(a=calculate_button):
    if NBART is not None:
        if a == 'Geomedian':
            #geomedian transform
            print('Calculating Geomedian')
            global NBART_GM
            NBART_GM=GeoMedian().compute(NBART)
            print('Geomedian Calculated')
        else:
            global NBART_TC
            print('Calculating Tasselled Cap Stats')
            NBART_TC=TCWStats().compute(NBART)
            print('Tasselled Cap Stats Calculated')

Plot the stat composite in RGB:

[18]:
plotButton = widgets.RadioButtons(description="Plot Which Stat?",
                                 options=['Geomedian','TCWetness Summary','TCI RGB Summary','TCI means RGB',
                                         'TCI std dev RGB'])
@interact
def on_plot_button_click(b=plotButton):
    global NBART
    global NBART_GM
    global NBART_TC
    if NBART is not None:
        if b == 'Geomedian':
            if NBART_GM is not None:
                DEAPlotting.three_band_image(NBART_GM, title='Geomedian composite', contrast_enhance=True)
            else:
                print('Please calculate Geomedian first')
        else:
            if NBART_TC is None:
                print('Please calculate Tasselled Cap Stats first')
            else:
                if b == 'TCWetness Summary':
                    plt.figure(figsize=(10,10))
                    ax = plt.gca()
                    im = ax.imshow(NBART_TC.pct_exceedance_wetness, cmap='gist_earth_r',vmin=0,vmax=1)
                    # create an axes on the right side of ax. The width of cax will be 5%
                    # of ax and the padding between cax and ax will be fixed at 0.05 inch.
                    divider = make_axes_locatable(ax)
                    cax = divider.append_axes("right", size="3%", pad=0.2)
                    plt.colorbar(im, cax=cax)
                    plt.show()
                if b == 'TCI RGB Summary':
                    DEAPlotting.three_band_image(NBART_TC, bands=['pct_exceedance_brightness',
                                                                  'pct_exceedance_greenness',
                                                                  'pct_exceedance_wetness'],
                                                 title=str(b),reflect_stand=1) #0.1 for high contrast)
                if b == 'TCI means RGB':
                    DEAPlotting.three_band_image(NBART_TC, bands=['mean_brightness',
                                                                  'mean_greenness',
                                                                  'mean_wetness'],
                                                 title=str(b),contrast_enhance=True)
                if b == 'TCI std dev RGB':
                    DEAPlotting.three_band_image(NBART_TC, bands=['std_brightness',
                                                                  'std_greenness',
                                                                  'std_wetness'],
                                                 title=str(b),contrast_enhance=True)
[17]:
savechoiceButton = widgets.RadioButtons(description="Save Which Stat?",
                                 options=['Geomedian','Tasselled Cap'])
@interact
def on_save_button_click(c=savechoiceButton):
    global NBART
    global NBART_GM
    global NBART_TC
    global OUTPUT_STAT
    global OUTPUT_STAT_NAME
    if NBART is not None:
        OUTPUT_STAT_NAME = c
        if c == 'Geomedian':
            if NBART_GM is not None:
                OUTPUT_STAT=NBART_GM
            else:
                print('Please calculate Geomedian first')
        else:
            if NBART_TC is None:
                print('Please calculate Tasselled Cap Stats first')
            else:
                OUTPUT_STAT=NBART_TC

Save as netcdf:

[14]:
saveFileB = SaveFileButton()
def handle_file_selector(sender):
    filename =saveFileB.file
    global OUTPUT_STAT
    global OUTPUT_STAT_NAME
    if OUTPUT_STAT is None:
        print('Calculate output stats first')
    else:
        DEADataHandling.write_your_netcdf(OUTPUT_STAT,OUTPUT_STAT_NAME,filename=filename, crs=OUTPUT_STAT.crs)
        print("saved netcdf ...check you added a .nc to the end of your filename..")
saveFileB.on_click(handle_file_selector)
display(saveFileB)

Save as GeoTiff:

[21]:
saveFileB2 = SaveFileButton()
def handle_file_selector(sender):
    filename =saveFileB2.file
    global OUTPUT_STAT
    global OUTPUT_STAT_NAME
    if OUTPUT_STAT is None:
        print('Calculate output stats first')
    write_geotiff(filename, OUTPUT_STAT)
    print("saved geotiff, make sure you added a .tif to your filename")
saveFileB2.on_click(handle_file_selector)
display(saveFileB2)

References

  1. Roberts, D., Dunn,B., Mueller, N. 2018, Open Data Cube Products Using High-Dimensional Statistics of Time Series, in press.

      1. Crist, A tm tasseled cap equivalent transformation for reflectance factor data, Remote Sensing of Environment, vol. 17, no. 3, pp. 301-306, 1985.

[ ]: