Tasselled Cap Wetness Index Calculation

Authors:

Bex Dunn

Created:

Feb 20, 2019

Last edited:

Feb 20, 2019

Requirements:

You need to run the following commands from the command line prior to launching jupyter notebooks from the same terminal so that the required libraries and paths are set:

module use /g/data/v10/public/modules/modulefiles

module load dea

This notebook assumes you have cloned the dea-notebooks repository by following the instructions on DEA notebooks using command line git or DEA notebooks using Github. If you didn’t, you will need to download any scripts that fail in the first cell into your relative directory so that the notebook will have access to those scripts. "../" means the directory above the directory you are working in, and "./" means the current directory.

If you find an error or bug in this notebook, please either create an ‘Issue’ in the Github repository, or fix it yourself and create a ‘Pull’ request to contribute the updated notebook back into the repository (See the repository README for instructions on creating a Pull request).

Background: Data from the Landsat 5,7 and 8 satellite missions are accessible through Digital Earth Australia (DEA).

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.This notebook produces the wetness index.

What does this notebook do?: This notebook takes a supplied shapefile of a polygon and queries the datacube Landsat surface reflectance data. It calculates tasselled cap wetness. The results and input surface reflectance are output as geotiffs for each individual scene.

Tags: plot,:index:gridspec, Landsat,:index:Landsat5,:index:Landsat7,:index:Landsat8, polygondrill, DEAPlotting, shapefile, geopolygon, datacube.utils.geometry, fiona, rasterio, query,:index:Scripts,:index:tasseled_cap, DEADataHandling, DEAPlotting, masking,:index:make mask,:index:load_clearlandsatwrite_geotiff, multi_timesteps_geotiffs

import modules and scripts

[1]:
import datacube
import datetime
import fiona
import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio.mask
import rasterio.features
from shapely import geometry
import seaborn as sns
import sys
import xarray as xr

import matplotlib.dates as mdates
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt

from datacube.storage import masking
from datacube.utils import geometry
from datacube.helpers import ga_pq_fuser, write_geotiff

sys.path.append('../10_Scripts')
import DEADataHandling, DEAPlotting, TasseledCapTools

dc = datacube.Datacube(app='tcw')

%load_ext autoreload

%autoreload 2

Set up polygon

[2]:
#change the path here if you want a different polygon
poly_path = '../Supplementary_data/Files/lake_george.shp'
[3]:
#open the polygon
with fiona.open(poly_path) as shapes:
        crs = geometry.CRS(shapes.crs_wkt)
        first_geometry = next(iter(shapes))['geometry']
        geom = geometry.Geometry(first_geometry, crs=crs)
[4]:
#plot polygon to check it looks ok
plt.clf()
shape_plot = gpd.read_file(poly_path)
shape_plot.plot()
plt.show()
<Figure size 432x288 with 0 Axes>
../../_images/notebooks_04_Index_calculation_Tasselled_Cap_Wetness_Index_Calculation_8_1.png

set up DEA query

[5]:
query = {'geopolygon': geom,
         'time': ('2019-01-01', '2019-08-30')
         }

Set up datasets

set cloudmasking threshold and load landsat nbart data

[6]:
#set cloudmasking threshold and load landsat nbart data
landsat_masked_prop = 0.99 #0 #remove scenes with more than 99% bad values.

#use the loadclearlandsat function to remove scenes that aren't nearly all clear - just so it looks pretty in the notebook and so cloud-shadow doesn't look 'wet'
ls578_ds = DEADataHandling.load_clearlandsat(dc=dc, query=query, product='nbart', masked_prop=landsat_masked_prop)
Loading ls5
    Skipping ls5; no valid data for query
Loading ls7
    Ignoring SLC-off observations for ls7
    Skipping ls7; no valid data for query
Loading ls8
    Loading 3 filtered ls8 timesteps
Returning ls8 data
    Replacing invalid -999 values with NaN (data will be coerced to float64)

mask the data with our original polygon to remove extra data

[7]:
data = ls578_ds
mask = rasterio.features.geometry_mask([geom.to_crs(data.geobox.crs)for geoms in [geom]],
                                           out_shape=data.geobox.shape,
                                           transform=data.geobox.affine,
                                           all_touched=False,
                                           invert=False)
[8]:
#for some reason xarray is not playing nicely with our old masking function
mask_xr = xr.DataArray(mask, dims = ('y','x'))
ls578_ds = data.where(mask_xr==False)

Plot our data as a false color image

[9]:
# plt.clf()
DEAPlotting.rgb(ls578_ds, bands=['swir1', 'nir', 'green'], col='time',col_wrap=4)
# plt.show()
../../_images/notebooks_04_Index_calculation_Tasselled_Cap_Wetness_Index_Calculation_18_0.png

Run the tasselled cap transform for our AOI

[10]:
#transform the nbart into tci. Note I've set a lowish threshold on the TCI but we're only really dealing with the unthresholded results in this notebook.
tci = TasseledCapTools.thresholded_tasseled_cap(ls578_ds,wetness_threshold=-1200, drop=True , drop_tc_bands=False)

get the wetness results

[11]:
tcw = tci['wetness']

plot the tasseled cap wetness for our AOI

[12]:
# plt.clf()
tcw.isel(time=1).plot(cmap='gist_earth_r', vmin=-4000, vmax=0)#note -4000 is very low
# plt.show()
[12]:
<matplotlib.collections.QuadMesh at 0x7f0450ced2b0>
../../_images/notebooks_04_Index_calculation_Tasselled_Cap_Wetness_Index_Calculation_24_1.png

have a look at areas where it might be wet (~>-350)

[21]:
# plt.clf()
tcw.isel(time=1).plot(cmap='gist_earth_r', vmin=-350, vmax=0)#note -4000 is very low
# plt.show()
[21]:
<matplotlib.collections.QuadMesh at 0x7f0450dc6978>
../../_images/notebooks_04_Index_calculation_Tasselled_Cap_Wetness_Index_Calculation_26_1.png

Write geotiffs of wetness for each time we have imagery

[13]:
#get attrs off our original dataset
tcw_dataset=tcw.to_dataset()
tcw_dataset.attrs=ls578_ds.attrs

set up paths to write out outputs to

[14]:
#get polygon name from the polygon path
polyname = poly_path.split('/')[-1].split('.')[0]

change the path below to a path where you’d like things saved!

[15]:
savefilepath = '/g/data/r78/rjd547/tmp/'+polyname
[16]:
#set dataset equal to wetness dataset
ds = tcw_dataset
#write a geotiff to file for each timestep
if len(ds.time)==1:
    print('one timestep')
    #drop the time dimension for only one timestep
    #write the dataset without the data percentage to file
    ds1 = ds.squeeze()
    #ds1 = ds1.drop('data_perc')
    write_geotiff(savefilepath+'_TCW_'+'.tif', ds1)

elif len(ds.time)>1:
    print('multiple timesteps')
    #remove data percentage as it breaks the geotiff writer
    #ds = ds.drop('data_perc')
    for timestep in range(len(ds.time)):
        timestep_date =np.datetime_as_string(ds.time.isel(time=timestep))[0:10]
        try:
            write_geotiff(savefilepath+'_TCW_'+timestep_date+'.tif', ds.isel(time =timestep))
            #complain if the file already exists but don't fail
            print('wrote to GeoTIFF' )
        except RuntimeError as err:
            print("RuntimeError: {0}".format(err))
multiple timesteps
wrote to GeoTIFF
wrote to GeoTIFF
wrote to GeoTIFF

Write geotiffs of surface reflectance for each time we have imagery

[17]:
for timestep in range(len(ds.time)):
    print(np.datetime_as_string(ds.time.isel(time=timestep))[0:10])
2019-01-13
2019-01-29
2019-03-02
[18]:
#set dataset equal to landsat dataset
ds = ls578_ds
#write a geotiff to file for each timestep
if len(ds.time)==1:
    print('one timestep')
    #drop the time dimension for only one timestep
    #write the dataset without the data percentage to file
    ds1 = ds.squeeze()
    #ds1 = ds1.drop('data_perc')
    write_geotiff(savefilepath+'_LS_'+'.tif', ds1)

elif len(ds.time)>1:
    print('multiple timesteps')
    #remove data percentage as it breaks the geotiff writer
    #ds = ds.drop('data_perc')
    for timestep in range(len(ds.time)):
        timestep_date =np.datetime_as_string(ds.time.isel(time=timestep))[0:10]
        try:
            write_geotiff(savefilepath+'_LS_'+timestep_date+'.tif', ds.isel(time =timestep))
            #complain if the file already exists but don't fail
            print('wrote to GeoTIFF' )
        except RuntimeError as err:
            print("RuntimeError: {0}".format(err))
multiple timesteps
wrote to GeoTIFF
wrote to GeoTIFF
wrote to GeoTIFF

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.

[ ]: