Tool for browsing landsat data in IPython

landsat.py

""" Tool for searching Landsat archives """

from urllib import request
from datetime import datetime
from xml.etree import ElementTree

from karta import Polygon
from karta.crs import LonLatWGS84
from IPython.display import Image, display

# XML documentation at http://earthexplorer.usgs.gov/EE/metadata.xsd


Landsat8     = "LANDSAT_8"
Landsat7     = "LANDSAT_ETM"
Landsat7SLC  = "LANDSAT_ETM_SLC_OFF"
Landsat45TM  = "LANDSAT_TM"
Landsat45MSS = "LANDSAT_MSS2"
Landsat13MSS = "LANDSAT_MSS1"
LandsatAll   = "LANDSAT_COMBINED"

class LandsatScene(object):
    schema = "{http://upe.ldcm.usgs.gov/schema/metadata}"
    def __init__(self, xml):
        self.xml = xml
        self.sceneid = xml.find(self.schema + "sceneID").text
        self.browseurl = xml.find(self.schema + "browseURL").text
        self.cloudcover = int(xml.find(self.schema + "cloudCover").text)
        self.cloudcover2 = float(xml.find(self.schema + "cloudCoverFull").text)
        self.row = int(xml.find(self.schema + "row").text)
        self.path = int(xml.find(self.schema + "path").text)
        self.starttime = xml.find(self.schema + "sceneStartTime").text
        self.poly = Polygon([(float(xml.find(self.schema + "upperLeftCornerLongitude").text),
                              float(xml.find(self.schema + "upperLeftCornerLatitude").text)),
                             (float(xml.find(self.schema + "lowerLeftCornerLongitude").text),
                              float(xml.find(self.schema + "lowerLeftCornerLatitude").text)),
                             (float(xml.find(self.schema + "lowerRightCornerLongitude").text),
                              float(xml.find(self.schema + "lowerRightCornerLatitude").text)),
                             (float(xml.find(self.schema + "upperRightCornerLongitude").text),
                              float(xml.find(self.schema + "upperRightCornerLatitude").text))],
                             crs=LonLatWGS84)
        return
    
    def __repr__(self):
        return "<LandsatScene: %s>" % self.sceneid

    def __str__(self):
        return ("ID: {ident}\nWRS: {row}/{path}\nDate: {date}\n"
                "Clouds: {clouds}".format(ident=self.sceneid, row=self.row,
                path=self.path, date=self.starttime, clouds=self.cloudcover2))

    @property
    def year(self):
        return int(self.starttime.split(":")[0])
    
    @property
    def day(self):
        return int(self.starttime.split(":")[1])
    
    @property
    def hour(self):
         return int(self.starttime.split(":")[2])
    
    def browse(self, **kw):
        """ Returns an IPython-embeddable view of the scene's sample image """
        return Image(self.browseurl, **kw)

def _append_query_constraints(s, **kw):
    """
    dates :: [2]datetime    Start and end dates for scenes

    maxclouds :: int        Maximum cloud cover, with 1 -> 10%, 2 -> 20%, etc.

    sensor :: str           Sensor/platform
    """

    if "dates" in kw:
        dates = kw["dates"]
    else:
        dates = (datetime(1900, 1, 1), datetime(2100, 1, 1))
    s = "&".join([s, "start_date=%s" % dates[0].strftime("%Y-%m-%d"),
                     "end_date=%s" % dates[1].strftime("%Y-%m-%d")])

    if "maxclouds" in kw:
        s = "&".join([s, "cc=%s" % kw["maxclouds"]])

    if "sensor" in kw:
        s = "&".join([s, "sensor=%s" & kw["sensor"]])

    return s

def search_by_lonlat(bbox=(-179.9999, 180.0, -90.0, 90.0), **kw):
    """ Query the USGS for XML bulk metadata.
    
    Arguments:
    ----------
    bbox    :: [4]float          [west, east, south, north]

    dates :: [2]datetime    Start and end dates for scenes

    maxclouds :: int        Maximum cloud cover, with 1 -> 10%, 2 -> 20%, etc.

    sensor :: str           Sensor/platform

    Returns xml.etree.ElementTree
    """


    s = ("http://earthexplorer.usgs.gov/EE/InventoryStream/latlong?"
         "north={n}&south={s}&east={e}&west={w}".format(w=bbox[0], e=bbox[1], s=bbox[2], n=bbox[3]))
    s = _append_query_constraints(s, **kw)
    response = request.urlopen(s)
    xml = ElementTree.fromstring(response.readall().decode("utf-8"))
    return xml

def search_by_rowpath(startpath, endpath, startrow, endrow, **kw):
    """ Query the USGS for XML bulk metadata.
    
    Arguments:
    ----------
    startrow :: int
    endrow :: int
    startpath :: int
    endpath :: int

    dates :: [2]datetime    Start and end dates for scenes

    maxclouds :: int        Maximum cloud cover, with 1 -> 10%, 2 -> 20%, etc.

    sensor :: str           Sensor/platform

    Returns xml.etree.ElementTree
    """

    s = ("http://earthexplorer.usgs.gov/EE/InventoryStream/pathrow?"
         "start_path={sp}&end_path={ep}&start_row={sr}&end_row={er}".format(
            sp=startpath, ep=endpath, sr=startrow, er=endrow))
    s = _append_query_constraints(s, **kw)
    response = request.urlopen(s)
    xml = ElementTree.fromstring(response.readall().decode("utf-8"))
    return xml

def get_scenes(xml):
    """ Extract Landsat scenes from an XML bulk metadata object.
    
    Returns []LandsatScene
    """
    scenes = []
    for child in xml.getchildren():
        if child.tag.endswith("metaData"):
            scenes.append(LandsatScene(child))
    return scenes