Working with specific remote sensing data

Landsat

Usage

Pygaarst offers an intuitive and friendly way to access Landsat scene data and metadata as distributed as at-sensor scaled radiance data files (“Level1”) in GeoTIFF format, one file per imaging band, as zipped tar archives.

The work is done by two classes: Landsatscene, which takes the name of the unzipped data directory, and Landsatband, which can be instantiated from a Landsatscene object:

>>> from pygaarst import raster
>>> basedir = "path/to/data"
>>> landsatscene = "LE70690142004201EDC01"
>>> sc = raster.Landsatscene(os.path.join(basedir, landsatscene))
>>> type(sc)
<class 'pygaarst.raster.Landsatscene'>
>>> b2 = sc.band2
>>> type(b2.data)
<type 'numpy.ndarray'>
>>> b2.data.shape
(8021, 8501)

The Landsatscene knows how to calculate often used radiometric quantities such as vegetation indices, or access the Landsat metadata. Both old and new metadata format is supported:

>>> ndvi = sc.NDVI
>>> type(ndvi)
<type 'numpy.ndarray'>
>>> sc.meta['IMAGE_ATTRIBUTES']['SUN_ELEVATION']
44.26873508
>>> sc.meta['PRODUCT_METADATA']['DATE_ACQUIRED']
datetime.date(2004, 7, 19)

Landsatscene also supports a filename infix in case all scene files have been pre-processed. A typical example would be if the GeoTIFF band files have been sub-setted and saved under a new file name that adds a label at the end of the filename, before the extension. For example:

>>> b2.data.shape
(8021, 8501)
>>> sc.infix = '_CLIP'
>>> b3 = sc.band3
>>> b3.data.shape
(347, 373)
>>> b2.filepath
'path/to/data/LE70690142004201EDC01_B2.TIF'
>>> b3.filepath
'path/to/data/LE70690142004201EDC01_B3_CLIP.TIF'

The Landsatband object knows its own geographic attributes and can translate the raw data into radiance and reflectance. For thermal IR bands, calculation of radiant (brightness) temperatures is supported:

>>> b3.radiance
array([[ 23.59606299,  24.83937008,  23.59606299, ...,  29.81259843,
         29.81259843,  29.81259843],
       [ 25.46102362,  24.21771654,  23.59606299, ...,  30.43425197,
         30.43425197,  29.19094488],
       ...,
       [ 41.0023622 ,  39.75905512,  40.38070866, ...,  90.11299213,
         91.35629921,  92.5996063 ],
       [ 40.38070866,  38.51574803,  41.0023622 , ...,  90.11299213,
         90.73464567,  90.73464567]])
>>> b3.Lat
array([[ 65.23978665,  65.23979005,  65.23979346, ...,  65.24086228,
         65.24086467,  65.24086706],
       [ 65.2400558 ,  65.2400592 ,  65.24006261, ...,  65.24113144,
         65.24113383,  65.24113622],
       ...,
       [ 65.33291157,  65.332915  ,  65.33291841, ...,  65.3339918 ,
         65.3339942 ,  65.3339966 ],
       [ 65.33318072,  65.33318414,  65.33318756, ...,  65.33426096,
         65.33426336,  65.33426576]])
>>> b6 = sc.band6L
>>> b6.tKelvin
array([[ 298.51857637,  298.51857637,  299.01776653, ...,  301.97175896,
         301.48420756,  301.48420756],
       [ 297.51409689,  298.01736166,  298.51857637, ...,  302.45745107,
         301.97175896,  301.48420756],
       ...,
       [ 294.96609241,  294.96609241,  294.96609241, ...,  290.23752626,
         290.77248753,  290.23752626],
       [ 294.96609241,  294.96609241,  294.96609241, ...,  290.77248753,
         290.77248753,  290.77248753]])

All parameters are, where provided, taken from the scene-specific metadata. Where not, they are taken from the published literature.

The Landsatscene class

class pygaarst.raster.Landsatscene(dirname)

A container object for TM/ETM+ L5/7 and OLI/TIRS L8 scenes

Arguments:
dirname (str): the full or relative file path that contains all files
NBR

The Normalized Burn Index

NDVI

The Normalized Difference Vegetation Index

TIRband

Label of a suitable thermal infrared band for the scene. Used in loops over Landsat scenes from multiple platform/sensor combinations. Will return B6 for L4/5, B6H for L7, B10 for L8.

ltkcloud

Cloud masking and landcover classification with the Luo–Trishchenko–Khlopenkov algorithm (doi:10.1109/LGRS.2010.2095409).

Returns:
A numpy array of the same shape as the input data The classes signify: - 1.0 - bare ground - 2.0 - ice/snow - 3.0 - water - 4.0 - cloud - 5.0 - vegetated soil
naivecloud

Heuristic cloud masking via thermal IR threshold

The Landsatband class

class pygaarst.raster.Landsatband(filepath, band=None, scene=None)

Represents a band of a Landsat scene.

Implemented: TM/ETM+ L5/7 and OLI/TIRS L8, old and new metadata format

newmetaformat

Checks if band uses old or new metadata format

radiance

Radiance in W/um/m^2/sr derived from DN and metadata, as numpy array

reflectance

Reflectance (0 .. 1) derived from DN and metadata, as numpy array

tKelvin

Radiant (brightness) temperature at the sensor in K, implemented for Landsat thermal infrared bands.

EO-1 Hyperion, ALI and general USGS Level 1 data

Hyperion and ALI

Hyperion imaging spectroscopy data and ALI multi-spectral imagery, processed to a Level 1 product, are distributed by the USGS in a format very similar to Landsat data. Pygaarst provides classes similar to the Landsat-specific classes for these sensors.

class pygaarst.raster.Hyperionscene(dirname)

A container object for EO-1 Hyperion scenes. Input: directory name, which is expected to contain all scene files.

get_datacube(outfn, bandlist, islice=None, jslice=None, set_fh=False)

Creates a rasterhelpers.Datacube object from a bandlist and the radiance data from the whole Hyperion scene.

Arguments:
outfn (str): file path of the HDF5 file that stores the cube bandlist (sequence of str): a list or array of band names islice (sequence of int): list or array of row indices jslice (sequence of int): list or array of column indices set_fh (bool): should an open filehandle be set as an argument?
spectrum(i_idx, j_idx, bands='calibrated', bdsel=[])

Calculates the radiance spectrum for one pixel.

Arguments:

i_idx (int): first coordinate index of the pixel j_idx (int): second coordinate index of the pixel bands (str): indicates the bands that are used

‘calibrated’ (default): only use calibrated bands ‘high’: use uncalibrated bands 225-242 ‘low’: use uncalibrated bands 1-7 ‘all’: use all available bands ‘selected’: use bdsel attribute or argument

bdsel: sequence data type containing band indices to select

class pygaarst.raster.Hyperionband(filepath, band=None, scene=None)

Represents a band of an EO-1 Hyperion scene.

radiance

Radiance in W / um / m^2 / sr derived from digital number and metadata, as numpy array

class pygaarst.raster.ALIscene(dirname)

A container object for EO-1 ALI scenes. Input: directory name, which is expected to contain all scene files.

class pygaarst.raster.ALIband(filepath, band=None, scene=None)

Represents a band of an EO-1 ALI scene.

radiance

Radiance in W / um / m^2 / sr derived from digital number and metadata, as numpy array

The USGSL1 parent classes

All of them inherit from parent classes called USGSL1scene and USGSL1band respectively.

class pygaarst.raster.USGSL1scene(dirname)

A container object for multi- and hyperspectral satellite imagery scenes as provided as Level 1 (at-sensor calibrated scaled radiance data) by various USGS data portals: Landsat 4/5 TM, Landsat 7 ETM+, Landsat 7 OLI/TIRS, EO1 ALI and EO1 Hyperion

Arguments:
dirname (str): name of directory that contains all scene files.
get_normdiff(label1, label2)

Calculate a generic normalized difference index

Arguments:
label1, label2 (str): valid band labels, usually of the form ‘bandN’
class pygaarst.raster.USGSL1band(filepath, band=None, scene=None)

Represents a band of a USGSL1scene.

This class is intended to be used via chlid classes: Landsatband, Hyperionband, ALIband

sensor

Returns sensor name

spacecraft

Returns spacecraft name (L4, L5, L7, L8)

The metatadata parser

All these USGS-provided satellite remote sensing data use essentially the same metadata format. It is provided in a text file that can be recognized by the letters MTL in the filename. The data structure is nested naturally maps onto a dictionary of dictionaries, with unlimited nesting depth. Unfortunately it appears to be quite badly documented, or in any event I have been unable to locate a data description or specification.

For objects of type pygaarst.raster.USGSL1scene or its children (pygaarst.raster.Landsatscene, pygaarst.raster.ALIscene or pygaarst.raster.Hyperionscene), the parsed metadata dictionary is provided in the meta attribute:

>>> print sc.meta
{'IMAGE_ATTRIBUTES': {'GEOMETRIC_RMSE_MODEL': 2.868, 'CLOUD_COVER': 31.0, 'GEOMETRIC_RMSE_MODEL_X': 1.73, 'SUN_AZIMUTH': 162.9104406, 'SUN_ELEVATION': 44.26873508, 'GROUND_CONTROL_POINTS_MODEL': 210, 'IMAGE_QUALITY': 9, 'GEOMETRIC_RMSE_MODEL_Y': 2.287}, 'RADIOMETRIC_RESCALING': {'RADIANCE_MULT_BAND_7': 0.044, 'RADIANCE_MULT_BAND_5': 0.126, 'RADIANCE_MULT_BAND_4': 0.969, ...

The parser is implemented via the module pygaarst.mtlutils. An arbitrary MTL file can be parsed by calling pygaarst.mtlutils.parsemeta():

pygaarst.mtlutils

Helpers for parsing MTL metadata files used by USGS Landsat and EO-1 (Hyperion and ALI) Level 1 GeoTIFF scene data archives

Created by Chris Waigl on 2014-04-20.

[2014-04-20] Refactoring original landsatutils.py, as MTL file format is also used by Hyperion and ALI.

exception pygaarst.mtlutils.MTLParseError

Custom exception: parse errors in Landsat or EO-1 MTL metadata files

pygaarst.mtlutils.parsemeta(metadataloc)

Parses the metadata.

Arguments:
metadataloc: a filename or a directory.

Returns metadata dictionary

MODIS swath data

[This is a little harder...]

Suomi/NPP VIIRS

class pygaarst.raster.VIIRSHDF5(filepath, geofilepath=None, variable=None)

A class providing access to a VIIRS SDS HDF5 file or dataset Parameters: filepath: full or relative path to the data file geofilepath (optional): override georeference array file from metadata; full or relative path to georeference file variable (optional): name of a variable to access

geodata

Object representing the georeference data, in its entirety

lats

Latitudes as provided by georeference array

lons

Longitudes as provided by georeference array

MODAPS data download API client

This is a REST-full web service client that implements the NASA LAADSWEB data API (http://ladsweb.nascom.nasa.gov/data/api.html). Usage:

>>> from pygaarst import modapsclient as m
>>> a = m.ModapsClient()
>>> retvar = a.[methodname](args)

Module pygaarst.modapsclient

A client to do talk to MODAPS web services. See http://ladsweb.nascom.nasa.gov/data/web_services.html

Created by Chris Waigl on 2013-10-22.

class pygaarst.modapsclient.ModapsClient

Implements a client for MODAPS web service retrieval of satellite data, without post-processing

See http://ladsweb.nascom.nasa.gov/data/quickstart.html

getAllOrders(email)

All orders for an email address

getBands(product)

Available bands for a product

getBrowse(fileId)

fileIds is a single file-ID

getCollections(product)

Available collections for a product

getDataLayers(product)

Available data layers for a product

getDateCoverage(collection, product)

Available dates for a collection/product combination

TODO: add some result postprocessing - not a good format

getFileOnlineStatuses(fileIds)

fileIds is a comma-separated list of file-IDs

getFileProperties(fileIds)

fileIds is a comma-separated list of file-IDs

getFileUrls(fileIds)

fileIds is a comma-separated list of file-IDs

getMaxSearchResults()

Max number of search results that can be returned

getOrderStatus(OrderID)

Order status for an order ID. TODO: implement

getOrderUrl(OrderID)

Order URL(?) for order ID. TODO: implement

getPostProcessingTypes(products)

Products: comma-concatenated string of valid product labels

listCollections()

Lists all collections. Deprecated: use getCollections

listMapProjections()

Lists all available map projections

listProductGroups(instrument)

Lists all available product groups

listProducts()

Lists all available products

listProductsByInstrument(instrument, group=None)

Lists all available products for an instrument

listReprojectionParameters(reprojectionName)

Lists reprojection parameters for a reprojection

listSatelliteInstruments()

Lists all available satellite instruments

orderFiles(FileIDs)

Submits an order. TODO: implement

searchForFiles(products, startTime, endTime, north, south, east, west, coordsOrTiles=u'coords', dayNightBoth=u'DNB', collection=5)

Submits a search based on product, geography and time

searchForFilesByName(collection, pattern)

Submits a search based on a file name pattern