Tiled time series

README.md

A demo of how to calculate the NDVI and compute statistics across the z-axis (time) in a tiled/chunking workflow.

ndvi_time_series_test.py

#!/usr/bin/env python

from datetime import datetime
import numpy
import gdal
from EOtools.DatasetDrivers import stacked_dataset as sd
from EOtools import tiling as tl
from EOtools import band_math as bm
from EOtools import pq_utils as pqu
from EOtools.stats import temporal_stats as ts
from gaip import GriddedGeoBox

"""
The time taken to calculate the NDVI and then derive stats from the
NDVI time series is miniscule compared to the data read times.
"""

red_fname = '/g/data1/v10/testing_ground/jps547/time_series/NBAR_144_-035_R.vrt'
nir_fname = '/g/data1/v10/testing_ground/jps547/time_series/NBAR_144_-035_NIR.vrt'
pq_fname = '/g/data1/v10/testing_ground/jps547/time_series/PQA_144_-035_PQA.vrt'

out_fname = '/g/data1/v10/testing_ground/jps547/time_series/ndvi_time_series_stats'

# Initialise the time series datasets
red_ds = sd.StackedDataset(red_fname)
nir_ds = sd.StackedDataset(nir_fname)
pq_ds = sd.StackedDataset(pq_fname)

# Get a geobox
geobox = GriddedGeoBox.from_dataset(gdal.Open(red_fname))

# Get the tiles to iterate over
# Inbuilt method with default settings but can be re-initialised with custom
# settings
#tiles = red_ds.tiles
# Or access the function directly
tiles = tl.generate_tiles(red_ds.samples, red_ds.lines, red_ds.samples, 60)

# Initialise the output stats file
out_ds = tl.TiledOutput(out_fname, samples=red_ds.samples, lines=red_ds.lines,
                        bands=14, dtype=6, nodata=numpy.nan, geobox=geobox)

# The number of bands to process i.e. all bands
bands_list = range(1, red_ds.bands + 1)

# Define the expression we wish to calculate
exp = "(b1 - b2) / (b1 + b2)"

st = datetime.now()
for tile in tiles:
    # Read the spectral data and calculate the NDVI
    nir_subset = (nir_ds.read_tile(tile, raster_bands=bands_list)).astype('float32')
    red_subset = red_ds.read_tile(tile, raster_bands=bands_list)
    var_key = {'b1': nir_subset, 'b2': red_subset}
    ndvi_result = bm.band_math(exp, var_key)

    # Extract the PQ data and apply the default bit extraction
    pq_subset = pq_ds.read_tile(tile, raster_bands=bands_list)
    pq_mask = pqu.extract_pq_flags(pq_subset)

    # Set PQ'd pixels to a null -999
    ndvi_result[~pq_mask] = -999

    # Compute stats across the z-axis (Internally turn -999 to NaN)
    stats = ts.temporal_stats(ndvi_result, no_data=-999)

    # Write the stats tile (with 14 bands) to disk
    out_ds.write_tile(stats, tile)

    et = datetime.now()
    print "Time taken for tile {}: {}".format(tile, et-st)
    st = datetime.now()

# Close the output file
out_ds.close()