Spectral Products: Vegetation, Water, Urbanization


Notebook Summary

This notebook can be used to calculate several common spectral products from cloud-filtered Landsat mosaic images. These spectral products include: fractional cover, NDBI, NDVI, NDWI, SAVI, EVI, WOFS (Water Observations from Space algorithm from Australia). The notebook also creates a threshold product for any spectral product and can output results in GeoTIFF format.


Index

Import Dependencies and Connect to the Data Cube [▴](#top)


In [1]:
# Enable importing of utilities.
import sys
sys.path.append('..')

# Supress Warning 
import warnings
warnings.filterwarnings('ignore')

import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt

# Load Data Cube Configuration
import datacube
import utils.data_cube_utilities.data_access_api as dc_api  
api = dc_api.DataAccessApi()
dc= api.dc

Choose Platforms and Products [▴](#top)

List available products for each platform


In [2]:
# Get available products
products_info = dc.list_products()

# List LANDSAT 7 products
print("LANDSAT 7 Products:")
products_info[["platform", "name"]][products_info.platform == "LANDSAT_7"]


LANDSAT 7 Products:
Out[2]:
platform name
id
2 LANDSAT_7 ls7_usgs_sr_scene

In [3]:
# List LANDSAT 8 products
print("LANDSAT 8 Products:")
products_info[["platform", "name"]][products_info.platform == "LANDSAT_8"]


LANDSAT 8 Products:
Out[3]:
platform name
id
1 LANDSAT_8 ls8_usgs_sr_scene

Choose products

<p style="color:red";>CHANGE INPUTS BELOW


In [4]:
# Select a Product and Platform

# product = "ls7_usgs_sr_scene"
# platform = "LANDSAT_7"

product = "ls8_usgs_sr_scene"
platform = "LANDSAT_8"

Get the Extents of the Cube [▴](#top)


In [5]:
from utils.data_cube_utilities.dc_load import get_product_extents
from utils.data_cube_utilities.dc_time import dt_to_str

full_lat, full_lon, min_max_dates = get_product_extents(api, platform, product)

# Print the extents of the combined data.
print("Latitude Extents:", full_lat)
print("Longitude Extents:", full_lon)
print("Time Extents:", list(map(dt_to_str, min_max_dates)))


Latitude Extents: (-12.63361111121218, 18.40166666681388)
Longitude Extents: (-25.47250000020378, 44.01000000035208)
Time Extents: ['2013-03-21', '2020-01-27']

Visualize the available area


In [6]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(full_lat, full_lon)


Out[6]:

Define the Extents of the Analysis [▴](#top)

<p style="color:red";>CHANGE INPUTS BELOW


In [9]:
# Select an analysis region (Lat-Lon) within the extents listed above. 
# Select a time period (Min-Max) within the extents listed above (Year-Month-Day)
# This region and time period will be used for the cloud assessment

# Mount Kenya
# latitude = (-0.50, 0.20) 
# longitude = (37.0, 37.7) 

# Mombasa, Kenya
# latitude = (-4.1, -3.9) 
# longitude = (39.5, 39.7) 

# Nairobi, Kenya
# latitude = (-1.4347, -1.151) 
# longitude = (36.6985, 36.9978) 

# Kigoma, Tanzania
latitude = (-4.95, -4.84) 
longitude = (29.59, 29.75) 

# Dar es Salaam, Tanzania
# latitude = (-7.0, -6.7)
# longitude = (39.1, 39.4)

# Freetown, Sierra Leone
# latitude = (8.3267, 8.5123)
# longitude = (-13.3039, -13.139)

# Northern portion of Lake Naivasha, Kenya
# latitude = (-0.7355, -0.6901)
# longitude = (36.3248, 36.3735)

# Phu Thuy, Vietnam
# latitude = (10.8945, 10.9679)
# longitude = (108.0616, 108.1460)

# Time Period
time_extents = ('2018-01-01', '2018-12-31')

Visualize the selected area


In [10]:
display_map(latitude,longitude)


Out[10]:

Load, Clean, and Composite Data from the Data Cube [▴](#top)


In [11]:
landsat_dataset = dc.load(latitude = latitude,
                          longitude = longitude,
                          platform = platform,
                          time = time_extents,
                          product = product,
                          measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'])

In [12]:
# Displays the first few values of each data array to check the content
# Latitude and Longitude numbers = number of pixels in each dimension
# Time = number of time slices in the dataset

landsat_dataset


Out[12]:
<xarray.Dataset>
Dimensions:    (latitude: 397, longitude: 577, time: 68)
Coordinates:
  * time       (time) datetime64[ns] 2018-01-09T08:08:48.138332 ... 2018-12-27T08:08:58.772967
  * latitude   (latitude) float64 -4.84 -4.84 -4.84 ... -4.949 -4.95 -4.95
  * longitude  (longitude) float64 29.59 29.59 29.59 29.59 ... 29.75 29.75 29.75
Data variables:
    red        (time, latitude, longitude) int16 1975 1930 1880 ... 512 458 549
    green      (time, latitude, longitude) int16 2102 2027 1948 ... 772 789 805
    blue       (time, latitude, longitude) int16 1986 1959 1943 ... 461 457 502
    nir        (time, latitude, longitude) int16 2107 2082 2058 ... 3857 3432
    swir1      (time, latitude, longitude) int16 2150 2040 1897 ... 2145 2285
    swir2      (time, latitude, longitude) int16 1958 1875 1762 ... 1007 1168
    pixel_qa   (time, latitude, longitude) uint16 480 480 480 ... 322 322 322
Attributes:
    crs:      EPSG:4326

Create 2 masks - one that only retains land, and one that retains land and water.


In [13]:
from utils.data_cube_utilities.clean_mask import landsat_qa_clean_mask

cloud_mask = landsat_qa_clean_mask(landsat_dataset, platform=platform)
land_mask = landsat_qa_clean_mask(landsat_dataset, platform=platform, cover_types=['clear'])

# Land and Water Dataset = Land and Water pixels with NO Clouds and NO Cloud Shadows
land_and_water_dataset = landsat_dataset.where(cloud_mask)

# Land Dataset = Land ONLY pixels with NO Clouds, NO Cloud Shadows and NO Water pixels
land_dataset = landsat_dataset.where(land_mask)

Create composite "images" with mosaic functions.

<p style="color:red";>CHANGE INPUTS BELOW


In [14]:
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic, create_max_ndvi_mosaic, create_hdmedians_multiple_band_mosaic

# Select a compositing method to create your cloud-filtered mosaic
# Remove the comments from the pair of lines under one of the mosaic types
# Options are: Median, Geomedian, or Max_NDVI 

# This is the MEDIAN mosaic
land_and_water_composite = create_median_mosaic(land_and_water_dataset, cloud_mask)
land_composite = create_median_mosaic(land_dataset, land_mask)

# This is the GEOMEDIAN mosaic
# land_and_water_composite = create_hdmedians_multiple_band_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_hdmedians_multiple_band_mosaic(land_dataset, land_mask)

# This is the MAX_NDVI mosaic
# land_and_water_composite = create_max_ndvi_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_max_ndvi_mosaic(land_dataset, land_mask)

cloud_mask_composite = cloud_mask.max('time')

Show the land and water composite.


In [15]:
from utils.data_cube_utilities.dc_rgb import rgb

# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green

rgb(land_and_water_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=5000, width=10)
plt.show()


Show the land-only composite (water will be removed ... black pixels).


In [16]:
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green

rgb(land_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=5000, width=10)
plt.show()


Fractional Cover [▴](#top)

Fractional Cover (FC) is used for landcover type estimation (vegetation, non-green vegetation, bare soil) of each pixel.
We use a model from CSIRO (Juan Gerschmann) and apply it to a median mosaic.

There are 3 values for every pixel: Bare Soil (bs), Photosynthetic Vegetation (pv), and Non-Photosynthetic Vegetation (npv)
The range of each value is [0,100] - the integer percent of confidence in the pixel being the corresponding class.
Plot a False Color RGB result where RGB = bs/pv/npv


In [18]:
from utils.data_cube_utilities.dc_fractional_coverage_classifier import frac_coverage_classify 
frac_classes = frac_coverage_classify(land_composite)

In [19]:
# Plot of Fractional Cover
# RED = Bare Soil or Urban Areas
# BLUE = Non-Green Vegetation
# GREEN = Green Vegetation
# BLACK = Water

rgb(frac_classes, bands = ['bs', 'pv', 'npv'], min_possible=0, max_possible=100) 
plt.show()


Land Spectral Indices [▴](#top)


In [20]:
def NDBI(dataset):
    """
    The Normalized Difference Built-up Index (NDBI) monitors urbanization.
    Greater values generally indicate more urbanization.
    The range is [-1,1].
    """
    return (dataset.swir1 - dataset.nir)/(dataset.swir1 + dataset.nir)

In [21]:
def NDVI(dataset):
    """
    The Normalized Difference Vegetation Index (NDVI) monitors photosynthetic vegetation.
    Greater values generally indicate more photosynthetic vegetation.
    The range is [-1,1].
    """
    return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)

In [22]:
def SAVI(ds, L=0.5, normalize=True):
    """
    The Soil-Adjusted Vegetation Index (SAVI) monitors photosyntehtic vegetation, 
    correcting for soil differences.
    Greater values generally indicate more photosynthetic vegetation.
    For Landsat data, returned values should be in the range [-1,1] if `normalize == True`.
    If `normalize == False`, returned values should be in the range [-1-L,1+L].
    """
    savi = (ds.nir - ds.red) / (ds.nir + ds.red + L) * (1 + L)
    if normalize:
        savi.values = np.interp(savi.values, (-1-L, 1+L), (-1, 1))
    return savi

In [23]:
def EVI(ds, G=2.5, C1=6, C2=7.5, L=1, normalize=True):
    """
    The 3-band Enhanced Vegetation Index (EVI). EVI is superior to NDVI in accuracy 
    because it is less dependent on the solar incidence angle, atmospheric conditions 
    (e.g. particles and clouds), shadows, and soil appearance.
    Greater values generally indicate more photosynthetic vegetation.
    The range is [-1,1] if `normalize`, and [-1,2.5] otherwise..
    """
    evi = G * (ds.nir - ds.red) / (ds.nir + C1 * ds.red - C2 * ds.blue + L)
    # Clamp values to the range [-1,2.5].
    evi.values[evi.values < -1] = -1
    evi.values[2.5 < evi.values] = 2.5
    if normalize:
        # Scale values in the  range [0,2.5] to the range [0,1].
        pos_vals_mask = 0 < evi.values
        evi.values[pos_vals_mask] = np.interp(evi.values[pos_vals_mask], (0, 2.5), (0, 1))
    return evi

In [24]:
def NDWI(data, normalize=False, band_pair=0):
    """
    Computes the Normalized Difference Water Index for an `xarray.Dataset`.
    The range is [-1,1] if `not normalize`, and [0,1] otherwise.
    """
    return (data.nir - data.swir1) / (data.nir + data.swir1)

In [25]:
ndbi = NDBI(land_composite)  # Normalized Difference Build Up (Urbanization) Index
ndvi = NDVI(land_composite)  # Normalized Difference Vegetation Index
ndwi = NDWI(land_composite) # Normalized Difference Water Index
ndbi2 = NDBI(land_and_water_composite)  # Normalized Difference Build Up (Urbanization) Index
ndvi2 = NDVI(land_and_water_composite)  # Normalized Difference Vegetation Index
ndwi2 = NDWI(land_and_water_composite) # Normalized Difference Water Index
savi = SAVI(land_composite)  # Soil Adjusted Vegetation Index 
evi = EVI(land_composite) # Enhanced Vegetation Index

In [26]:
ds_ndvi = ndvi2.to_dataset(name = "NDVI")
ds_ndwi = ndwi2.to_dataset(name=  "NDWI")
ds_ndbi = ndbi2.to_dataset(name = "NDBI")
normalization_dataset = ds_ndvi.merge(ds_ndwi).merge(ds_ndbi)

Show NDBI, NDVI, NDWI False Color Composite.


In [27]:
# Plot of RGB = NDBI-NDVI-NDWI
# RED = Bare Soil or Urban Areas
# GREEN = Vegetation
# BLUE = Water

rgb(normalization_dataset, bands = ['NDBI','NDVI','NDWI'], min_possible=0, max_possible=1.0) 
plt.show()


Show the NDBI Composite.


In [28]:
%pylab inline
pylab.rcParams['figure.figsize'] = (10, 10)


Populating the interactive namespace from numpy and matplotlib

In [29]:
# May try cmap=GnBu or cmap=Greys, or cmap=Greens
(ndbi).plot(figsize=(12,10),cmap = "Reds", vmin=-1.0, vmax=1.0)
plt.show()


Show the NDVI Composite.


In [30]:
(ndvi).plot(figsize=(12,10), cmap = "Greens", vmin=-1.0, vmax=1.0)
plt.show()


Show the NDVI Composite, Including Water (usually very low NDVI for water).


In [31]:
# Create a custom colour map for NDVI
# Water (blue) = NDVI -1.0 to 0.0
# Urban or Bare Soil (brown) = NDVI 0.0 to 0.1
# Low Vegetation (tan) = NDVI 0.1 to 0.4
# Croplands (light green) = NDVI 0.4 to 0.6
# Dense Vegetation / Forests (dark green) = NDVI 0.6 to 1.0

ndvi_cmap = mpl.colors.ListedColormap(['blue', '#a52a2a','#ffffcc' ,  '#2eb82e',  '#006600'])
ndvi_bounds = [-1, 0, 0.1,  0.4,  0.6, 1]
ndvi_norm = mpl.colors.BoundaryNorm(ndvi_bounds, ndvi_cmap.N)

In [32]:
plt.figure(figsize=(10,10))
plt.imshow(ndvi2, cmap = ndvi_cmap, norm = ndvi_norm)
plt.show()


Show the EVI Composite.


In [33]:
(evi).plot(figsize=(10,10), cmap = "Greens", vmin=-1.0, vmax=2.5)
plt.show()


Show the SAVI Composite.


In [34]:
(savi).plot(figsize=(10,10), cmap = "Greens", vmin=-1.0, vmax=1.0)
plt.show()


Water Products [▴](#top)

Normalized Difference Water Index (NDWI)


In [35]:
def NDWI(dataset):
    """
    The Normalized Difference Water Index (NDWI) monitors water.
    Greater values generally indicate higher probability of water.
    The range is [-1,1].
    """
    return (dataset.green - dataset.nir)/(dataset.green + dataset.nir)

In [36]:
(ndwi2).plot(figsize=(10,10), cmap = "Blues", vmin=-1.0, vmax=1.0)
plt.show()


Water Observation from Space (WOFS)

WOFS is a water classifier developed by Geoscience Australia


In [37]:
from utils.data_cube_utilities.dc_water_classifier import wofs_classify

In [38]:
water_classification = wofs_classify(land_and_water_composite, mosaic = True)

In [39]:
# Plot of WOFS product
# BLUE = 1.0 = Water
# WHITE = 0.0 = Non-Water

water_classification.wofs.plot(cmap='Blues')
plt.show()


Create a Threshold Plot [▴](#top)

First we will define a minimum threshold and a maximum threshold. Then you will create a plot that colors the region between the threshold a single color (e.g. red) and the region outside the threshold will be BLACK or WHITE. Also, we will calculate the % of pixels and the number of pixels in the threshold range.


In [40]:
def threshold_plot(da, min_threshold, max_threshold, mask = None, width = 10, *args, **kwargs): 
    color_in    = np.array([255,0,0])
    color_out   = np.array([0,0,0])
    color_cloud = np.array([255,255,255])
    
    array = np.zeros((*da.values.shape, 3)).astype(np.int16)
    
    inside  = np.logical_and(da.values > min_threshold, da.values < max_threshold)
    outside = np.invert(inside)
    masked  = np.zeros(da.values.shape).astype(bool) if mask is None else mask
    
    array[inside] =  color_in
    array[outside] = color_out
    array[masked] =  color_cloud

    def figure_ratio(ds, fixed_width = 10):
        width = fixed_width
        height = len(ds.latitude) * (fixed_width / len(ds.longitude))
        return (width, height)


    fig, ax = plt.subplots(figsize = figure_ratio(da,fixed_width = width))
    
    x_tick_inds = np.linspace(0, len(da.longitude)-1, 10)
    x_tick_vals = np.interp(x_tick_inds, np.arange(len(da.longitude)), da.longitude.values)
    ax.set_xticks(x_tick_inds)
    ax.set_xticklabels(["{0:.3f}".format(lon) for lon in x_tick_vals], rotation=30)
    y_tick_inds = np.linspace(0, len(da.latitude)-1, 10)
    y_tick_vals = np.interp(y_tick_inds, np.arange(len(da.latitude)), da.latitude.values)
    ax.set_yticks(y_tick_inds)
    ax.set_yticklabels(["{0:.3f}".format(lat) for lat in y_tick_vals])
    
    plt.title("Threshold: {} < x < {}".format(min_threshold, max_threshold))
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    
    plt.imshow(array, *args, **kwargs)
    plt.show()

<p style="color:red";>CHANGE INPUTS BELOW


In [41]:
# Select a threshold range for your spectral variable and generate a plot.
# Remove comments from the set of 3 lines for your desired variable.
# Variable choices are NDBI, NDVI, EVI, SAVI, FC-Bare Soil, FC-Photosynthetic Vegetation.

# NDBI (Buildup Index) = -1.0 to 1.0 (full range)
# NDBI 0.0 to 0.2 is typical for urban areas
# -----------------------
# minimum_threshold = 0.0
# maximum_threshold = 0.3
# threshold_plot(ndbi, minimum_threshold, maximum_threshold, width = 8)

# NDVI (Vegetation Index) = -1.0 to 1.0
# NDVI < 0.0 = non-vegetation (bare soil)
# NDVI 0.2 to 0.6 = grasslands
# NDVI 0.6 to 0.9 = dense vegetation / trees
# -----------------------
minimum_threshold = 0.2
maximum_threshold = 0.6
threshold_plot(ndvi, minimum_threshold, maximum_threshold, width = 8)

# EVI (Vegetation Index) = -1.0 to 2.5
# EVI 2.0 to 2.5 is typical for dense vegetation
# -----------------------
# minimum_threshold = 2.0
# maximum_threshold = 2.5
# threshold_plot(evi, minimum_threshold, maximum_threshold, width = 8)

# SAVI (Vegetation Index) = -1.0 to 1.0
# -----------------------
# minimum_threshold = 0.6
# maximum_threshold = 0.9
# threshold_plot(savi, minimum_threshold, maximum_threshold, width = 8)

# Fractional Cover (pv,npv,bs) = 0 to 100
# Bare Soil (bs) >40 = urbanization / bare soil
# ----------------------
# minimum_threshold = 40.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.bs, minimum_threshold, maximum_threshold, width = 8)

# Fractional Cover (pv,npv,bs) = 0 to 100
# Photosynthetic Vegetation (pv) >80 = dense green vegetation
# ----------------------
# minimum_threshold = 80.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.pv, minimum_threshold, maximum_threshold, width = 8)



In [42]:
def threshold_count(da, min_threshold, max_threshold, mask = None):
    def count_not_nans(arr):
        return np.count_nonzero(~np.isnan(arr))
    
    in_threshold = np.logical_and( da.values > min_threshold, da.values < max_threshold)
    
    total_non_cloudy = count_not_nans(da.values) if mask is None else np.sum(mask.values)
    
    return dict(total = np.size(da.values),
                total_non_cloudy = total_non_cloudy,
                inside = np.nansum(in_threshold),
                outside = total_non_cloudy - np.nansum(in_threshold)
               )    
    
def threshold_percentage(da, min_threshold, max_threshold, mask = None):
    counts = threshold_count(da, min_threshold, max_threshold, mask = mask)
    return dict(percent_inside_threshold = (counts["inside"]   / counts["total"]) * 100.0,
                percent_outside_threshold = (counts["outside"] / counts["total"]) * 100.0,
                percent_clouds = ( 100.0-counts["total_non_cloudy"] / counts["total"] * 100.0))

<p style="color:red";>CHANGE INPUTS BELOW


In [43]:
# Select a threshold statistical function that matches your land spectral variable
# COUNT = number of pixels in each category
# PERCENTAGE = percent of pixels in each category
# ---------------------------------

# NDBI Threshold
threshold_count(ndbi,minimum_threshold,maximum_threshold, cloud_mask_composite)
# threshold_percentage(ndbi,minimum_threshold,maximum_threshold)

# NDVI Threshold
# threshold_count(ndvi,minimum_threshold,maximum_threshold)
# threshold_percentage(ndvi,minimum_threshold,maximum_threshold)

# EVI Threshold
# threshold_count(evi,minimum_threshold,maximum_threshold)
# threshold_percentage(evi,minimum_threshold,maximum_threshold)

# SAVI Threshold
# threshold_count(savi,minimum_threshold,maximum_threshold)
# threshold_percentage(savi,minimum_threshold,maximum_threshold)

# Fractional Cover - Bare Soil
# threshold_count(frac_classes.bs, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.bs, minimum_threshold, maximum_threshold)

# Fractional Cover - Photosynthetic Vegetation
# threshold_count(frac_classes.pv, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.pv, minimum_threshold, maximum_threshold)


Out[43]:
{'total': 229069, 'total_non_cloudy': 229069, 'inside': 139, 'outside': 228930}

In [44]:
threshold_percentage(ndbi,minimum_threshold,maximum_threshold)


Out[44]:
{'percent_inside_threshold': 0.06068040634044764,
 'percent_outside_threshold': 79.68472381684121,
 'percent_clouds': 20.254595776818334}

Water Counts


In [45]:
# NDWI value for water = 0.0 to 1.0

ndwi_minimum_threshold = 0.0
ndwi_maximum_threshold = 1.0

<p style="color:red";>CHANGE INPUTS BELOW


In [46]:
# Select a threshold "count" or "percentage" statistical result
# COUNT = number of pixels in each category
# PERCENTAGE = percent of pixels in each category

threshold_count(ndwi,ndwi_minimum_threshold,ndwi_maximum_threshold, cloud_mask_composite)
# threshold_percentage(ndwi,ndwi_minimum_threshold,ndwi_maximum_threshold)


Out[46]:
{'total': 229069,
 'total_non_cloudy': 229069,
 'inside': 123950,
 'outside': 105119}

In [47]:
threshold_percentage(ndwi,ndwi_minimum_threshold,ndwi_maximum_threshold, cloud_mask_composite)


Out[47]:
{'percent_inside_threshold': 54.110333567615,
 'percent_outside_threshold': 45.889666432385,
 'percent_clouds': 0.0}

In [48]:
# WOFS value for water = 0.0 or 1.0
# The threshold uses a range of 0.9 to 1.1 to identify water (1.0)

wofs_minimum_threshold = 0.0
wofs_maximum_threshold = 1.0

<p style="color:red";>CHANGE INPUTS BELOW


In [49]:
# Select a threshold "count" or "percentage" statistical result
# COUNT = number of pixels in each category
# PERCENTAGE = percent of pixels in each category

threshold_count(water_classification.wofs,wofs_minimum_threshold,wofs_maximum_threshold, cloud_mask_composite)
# threshold_percentage(water_classification.wofs,wofs_minimum_threshold,wofs_maximum_threshold)


Out[49]:
{'total': 229069, 'total_non_cloudy': 229069, 'inside': 0, 'outside': 229069}

In [50]:
threshold_percentage(water_classification.wofs,wofs_minimum_threshold,wofs_maximum_threshold, cloud_mask_composite)


Out[50]:
{'percent_inside_threshold': 0.0,
 'percent_outside_threshold': 100.0,
 'percent_clouds': 0.0}

Create GeoTIFF Output Products [▴](#top)

<p style="color:red";>CHANGE INPUTS BELOW


In [51]:
from utils.data_cube_utilities.dc_utilities import write_geotiff_from_xr
import os

# Remove the comment to create a GeoTIFF output product
# Change the name of the output file, or it will be overwritten for each run 
# Change the desired bands at the end of the function

output_dir = "output/geotiffs"
if not os.path.exists(output_dir):
    os.mkdir(output_dir)

# Fractional Coverage
# write_geotiff_from_xr(output_dir + "/frac_classes.tif", frac_classes, bands=['bs'])

# NDVI
# write_geotiff_from_xr(output_dir + "/ndvi_land.tif", ndvi)

# SAVI
# write_geotiff_from_xr(output_dir + "/savi.tif", savi)

# EVI
# write_geotiff_from_xr(output_dir + "/evi.tif", evi)

# WOFS
write_geotiff_from_xr(output_dir + "/wofs.tif", water_classification.wofs)

In [52]:
!ls -lah output/geotiffs/*.tif


-rw-r--r-- 1 jovyan users 1.8M Apr 24 01:07 output/geotiffs/wofs.tif