This script allows extracting Landsat data from the AGDC for a selected location

agdc_extractor.py

#!/g/data/v10/public/modules/dea-env/20181015/bin/python

import datacube as dc
from datacube.helpers import ga_pq_fuser
from datacube.storage import masking
import pyproj
import numpy as np
import xarray as xr
import datetime
import warnings
import argparse
import os

class BurnCube(dc.Datacube):
    def __init__(self):
        super(BurnCube, self).__init__(app='TreeMapping.getLandsatStack')

        self.band_names = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2']
        self.dataset = None
        self.geomed = None
        self.dists = None
        self.outlrs = None
        

    def to_netcdf(self, path):
        """Saves input data to a netCDF4 on disk
            Args:
                path string: path to a file on disk
        """
        self.dataset.to_netcdf(path)

    def open_dataset(self, path):
        """Loads input data from a netCDF4 on disk
            Args:
                path string: path to a file on disk
        """
        self.dataset = xr.open_dataset(path)

    def _load_pq(self, x, y, res, period, n_landsat):
        query = {
            'time': period,
            'x': x,
            'y': y,
            'crs': 'EPSG:3577',
            'measurements': ['pixelquality'],
            'resolution': res,
        }

        pq_stack = []
        for n in n_landsat:
            pq_stack.append(self.load(product='ls{}_pq_albers'.format(n),
                                      group_by='solar_day', fuse_func=ga_pq_fuser,
                                      resampling='nearest', **query))

        pq_stack = xr.concat(pq_stack, dim='time').sortby('time')

        # Land/sea mask isn't used at the moment. Possible alternatives are WOFS and ITEM.
        #pq_stack['land'] = masking.make_mask(pq_stack.pixelquality, land_sea='land')
        # masking cloud, saturation and invalid data (contiguous)
        pq_stack['good_pixel'] = masking.make_mask(pq_stack.pixelquality, cloud_acca='no_cloud',
                                                   cloud_fmask='no_cloud', cloud_shadow_acca='no_cloud_shadow',
                                                   cloud_shadow_fmask='no_cloud_shadow',
                                                   blue_saturated=False,
                                                   green_saturated=False,
                                                   red_saturated=False,
                                                   nir_saturated=False,
                                                   swir1_saturated=False,
                                                   swir2_saturated=False,
                                                   contiguous=True)

        return pq_stack

    def _load_nbart(self, x, y, res, period, n_landsat):
        query = {
            'time': period,
            'x': x,
            'y': y,
            'crs': 'EPSG:3577',
            'measurements': self.band_names,
            'resolution': res,
        }

        nbart_stack = []
        for n in n_landsat:
            print(n)
            dss = self.find_datasets(product='ls{}_nbart_albers'.format(n), **query)
            nbart_stack.append(self.load(product='ls{}_nbart_albers'.format(n),
                                         group_by='solar_day', datasets=dss, resampling='bilinear',
                                         **query))

        print(nbart_stack)
        nbart_stack = xr.concat(nbart_stack, dim='time').sortby('time')

        return nbart_stack

    def load_cube(self, x, y, res, period, n_landsat):
        """Loads the Landsat data for the selected region and sensors
            Note:
                This method loads the data into the self.dataset variable.
            Args:
                x list float: horizontal min max range
                y list float: vertical min max range
                res float: pixel resolution for the input data
                period list datetime: temporal range for input data
                n_landsat int: number of the Landsat mission
        """
        
        nbart_stack = self._load_nbart(x, y, res, period, n_landsat)
        pq_stack = self._load_pq(x, y, res, period, n_landsat)
        pq_stack, nbart_stack = xr.align(pq_stack, nbart_stack, join='inner')
        mask = np.nanmean(pq_stack.good_pixel.values.reshape(pq_stack.good_pixel.shape[0], -1), axis=1) > .2
    
        self.dataset = nbart_stack.sel(time=mask).where(pq_stack.good_pixel.sel(time=mask), 0, drop=False) # keep data as integer
        del self.dataset['red'].attrs['spectral_definition']
        del self.dataset['red'].attrs['crs']
        del self.dataset['green'].attrs['spectral_definition']
        del self.dataset['green'].attrs['crs']
        del self.dataset['blue'].attrs['spectral_definition']
        del self.dataset['blue'].attrs['crs']
        del self.dataset['nir'].attrs['spectral_definition']
        del self.dataset['nir'].attrs['crs']
        del self.dataset['swir1'].attrs['spectral_definition']
        del self.dataset['swir1'].attrs['crs']
        del self.dataset['swir2'].attrs['spectral_definition']
        del self.dataset['swir2'].attrs['crs']
        del self.dataset.attrs['crs']
        del self.dataset['time'].attrs['units']

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="""Modis Vegetation Analysis argument parser""")
    parser.add_argument('-lat', '--latitude', type=float, required=True, help="Latitude of the selected point.")
    parser.add_argument('-lon', '--longitude', type=float, required=True, help="Longitude of the selected point.")
    parser.add_argument('-n', '--landsat_n', type=int, required=False, default=8, help="Date with format YYYYMMDD or YYYY to update with latest data.")
    parser.add_argument('-dst', '--destination', required=True, type=str, help="Full path to destination.")
    args = parser.parse_args()

    # convert to projected centre coordinates
    wgs84 = pyproj.Proj(init='epsg:4326')
    gda94 = pyproj.Proj(init='epsg:3577')
    easting, northing = pyproj.transform(wgs84, gda94, args.longitude, args.latitude)

    # define projected region extent
    x = (easting+1500, easting-1500)   # 1,500 m is half of the required window size
    y = (northing+1500, northing-1500)

    bc = BurnCube()
    bc.load_cube(x, y, (25, 25), ('1989-01-01', '2019-01-01'), [args.landsat_n])

    bc.dataset.to_netcdf(args.destination)

To execute this on the VDI:

1.- Download and make this code executable
chmod u+x agdc_extractor.py

2.- Load the DEA module
module use /g/data/v10/public/modules/modulefiles
module load dea

3.- Run for the selected location. For example for Landsat 5 at a location in Canberra:
./agdc_extractor.py -lat -35.281 -lon 149.13 -n 5 -dst /g/data/fj4/scratch/Canberra_LS5.nc