select-landsat8

#!/usr/bin/env python

import os
import json
import gzip
import requests

from datetime import datetime
from dateutil.parser import parse

import click
import fiona as fio
from shapely.geometry import shape


class SceneTree(object):

    def __init__(self, scene_list):

        self.data = {
            '%03d' % p: {
                '%03d' % r: []
                for r in range(1, 249)
            } for p in range(1, 234)
        }

        for s in scene_list:
            path, row = s[3:6], s[6:9]
            self.data[path][row].append(s)

    def __contains__(self, scene):
        path, row = scene[3:6], scene[6:9]
        return True if scene in self.data[path][row] else False


def get_scene_list(cloud_cover):
    """
    Get the scene list hosted on landsat-pds. This will be used
    to ensure scenes are not ingested more than once.
    """

    scene_list_path = '/tmp/scene_list.gz'

    if not os.path.exists(scene_list_path):
        SCENE_LIST_URL = 'https://landsat-pds.s3.amazonaws.com/scene_list.gz'

        with open(scene_list_path, 'wb') as f:
            r = requests.get(SCENE_LIST_URL, stream=True)
            for block in r.iter_content(1024):
                f.write(block)

    with gzip.open(scene_list_path, 'rb') as f:

        # Pop the first line containing the header
        f.readline()

        scene_list = [
            s.decode('utf-8').split(',')[0]
            for s in f.readlines()
            if float(s.decode('utf-8').split(',')[2]) < cloud_cover
        ]

    return scene_list


def parse_date_input(ctx, param, value):
    return None if value is None else parse(value)


def get_pathrows(wrspath, region):
    """
    Select pathrows that intersect a region.
    """

    with fio.open(wrspath) as src:

        if region is None:
            for idx, item in src.items():
                path = str(item['properties']['PATH']).zfill(3)
                row = str(item['properties']['ROW']).zfill(3)
                yield (path, row)
        else:

            for idx, item in src.items():
                s = shape(item['geometry'])
                if s.intersects(region):
                    path = str(item['properties']['PATH']).zfill(3)
                    row = str(item['properties']['ROW']).zfill(3)
                    yield (path, row)


def get_acquisition_date(scene):
    return datetime.strptime(scene[9:16], '%Y%j')


@click.command('select-landsat8')
@click.option('--geojson', type=click.Path(exists=True), required=False, help='GeoJSON geometry for the desired region.')
@click.option('--acquired-after', required=False, callback=parse_date_input, help='A scene is acquired after this date.')
@click.option('--acquired-before', required=False, callback=parse_date_input, help='A scene is acquired before this date.')
@click.option('--cloud-cover', type=click.FLOAT, default=100.0)
def select_landsat8(geojson, acquired_after, acquired_before, cloud_cover):
    """
    Select Landsat 8 scenes that intersect with a specified geometry and date range.
    """

    # Get all the scenes available on landsat-pds
    scene_list = get_scene_list(cloud_cover)
    scene_tree = SceneTree(scene_list)

    # Assume the WRS2 descending shapefile is local
    wrspath = os.path.join('/', 'Users', os.environ['USER'], 'Data', 'wrs2_descending', 'wrs2_descending.shp')

    region = None

    if geojson is not None:

        # Open and parse the geometry
        with open(geojson) as f:
            data = json.load(f)
            region = shape(data)
    
    for path, row in get_pathrows(wrspath, region):

        for scene in scene_tree.data[path][row]:
            if (get_acquisition_date(scene) >= acquired_after) and (get_acquisition_date(scene) < acquired_before):
                click.echo(scene)


if __name__ == '__main__':
    select_landsat8()

./select-landsat8 --cloud-cover 50 --acquired-after 2016-12-21 --acquired-before 2017-03-20 > scenes.txt