lsru (Landsat Surface Reflectance Utils)¶
Query, order and download Landsat surface reflectance data programmatically
lsru allows interaction with Usgs and Espa APIs programmatically from python.
It has 3 main classes:
Usgsis the interface to the USGS json API. It is mostly used to query the Landsat catalog for available scenes intersecting with a given spatio-temporal window.Espais the interface to the ESPA API. Espa is a platform that proposes on demand pre-processing of Landsat data to surface reflectance. Orders can be placed directly from python using that class.Orderis the interface to each individual orders placed to the espa platform; it allows retrieving order status and downloading corresponding scenes.
lsru also contains various utilities to smoothen workflows for various use cases of the module.
User Guide¶
Setup¶
The package requires a configuration file in which usgs credentials are written.
By default the file is called ~/.lsru (this can be modified if you want to join
this configuration with the configuration of another project) and has the following structure.
[usgs]
username=your_usgs_username
password=your_very_secure_password
API reference¶
lsru¶
Interfaces to Usgs and Espa APIs
Usgs([version, conf]) |
Interface to the Usgs API |
Usgs.login() |
Login to the Usgs api |
Usgs.search(collection, bbox[, begin, end, …]) |
Perform a spatio temporal query on Landsat catalog |
Usgs.get_collection_name(num) |
Get Earth Explorer Landsat collection names |
Espa([conf]) |
Interface to the Espa API |
Espa.order(scene_list, products[, format, …]) |
Place a pre-procesing order to espa |
Espa.get_available_products(scene_list) |
Get the list of available products for each elements of a list of scene ids |
Order(orderid[, conf]) |
Class to deal with espa orders |
Order.download_all_complete(path[, unpack, …]) |
Download all completed scenes of the order to a folder |
Order.cancel() |
Cancel the order |
utils¶
utils.bounds(geom) |
Return a bounding box from a geometry |
utils.geom_from_metadata(meta) |
Return a geometry from a Landsat scene metadata as returned by USGS api |
utils.is_valid(id) |
Landsat scene id validity checker |
utils.url_retrieve(url, filename[, …]) |
Generic file download function |
utils.url_retrieve_and_unpack(url, path[, …]) |
Generic function to combine download and unpacking of tar archives |
Basic Usage¶
This example shows how to query the data from a bounding box, order surface reflectance processing for full scenes without reprojection and download.
Define variables and query available scenes to Usgs¶
from lsru import Usgs
import datetime
# Define query extent
bbox = (3.5, 43.4, 4, 44)
# Instantiate Usgs class and login
usgs = Usgs()
usgs.login()
# Query the Usgs api to find scene intersecting with the spatio-temporal window
scene_list = usgs.search(collection='LANDSAT_8_C1',
bbox=bbox,
begin=datetime.datetime(2013,1,1),
end=datetime.datetime(2016,1,1),
max_results=10,
max_cloud_cover=40)
# Extract Landsat scene ids for each hit from the metadata
scene_list = [x['displayId'] for x in scene_list]
Place a processing order to Espa¶
The scene list can be used to send a processing order to Espa via the Espa API.
from lsru import Espa
from pprint import print
# Instantiate Espa class
espa = Espa()
# Place order (full scenes, no reprojection, sr and pixel_qa)
order = espa.order(scene_list=scene_list, products=['sr', 'pixel_qa'])
print(order.orderid)
# espa-loic.dutrieux@wur.nl-10212018-102816-245'
Check current orders status¶
for order in espa.orders:
# Orders have their own class with attributes and methods
print('%s: %s' % (order.orderid, order.status))
# espa-loic.dutrieux@wur.nl-10222018-062836-330: ordered
# espa-loic.dutrieux@wur.nl-10212018-174321-508: complete
# espa-loic.dutrieux@wur.nl-10212018-174430-792: complete
# espa-loic.dutrieux@wur.nl-10212018-102816-245: complete
# espa-loic.dutrieux@wur.nl-10182018-100137-786: complete
Download completed orders¶
When Espa finishes pre-processing an order, its status
changes to complete, we can then download the processed scenes.
for order in espa.orders:
if order.is_complete:
order.download_all_complete('/media/landsat/download/dir')
Order data using a polygon¶
The following example details the steps to place a pre-processing order of scenes intersecting with a polygon.
Such task requires a little bit of boilerplate code since the Usgs API requires an extent and not a polygon, but thanks to shapely and some additional utilities provided by lsru it can be done with a few lines of code.
The steps are roughly to:
- Read the feature (e.g. using
fionaor directly from a geojson file withjson) - Compute the bounding box of the feature
- Query scenes intersecting with the bounding box to the Usgs API
- Filter out scenes that do not intersect with the geometry
- Place the pre-processing order to Espa
To run this script, we’ll need shapely, and some utils to manipulate geojson geometries provided by lsru. You may need fiona as well in case you need to read a feature from a geospatial vector file (shapefile, geopackage, etc)
from datetime import datetime
from shapely.geometry import shape
from lsru import Usgs, Espa, Order
from lsru.utils import bounds, geom_from_metadata
For the present example, we’ll use the following feature, which roughly corresponds to the contours of the state of Baja California Sur in Mexico.
feature = {'geometry': {'coordinates': [[[-114.192, 27.975],
[-114.456, 27.8],
[-115.093, 27.878],
[-114.543, 27.313],
[-113.928, 26.922],
[-113.423, 26.844],
[-112.72, 26.412],
[-112.148, 25.741],
[-112.324, 24.827],
[-111.687, 24.467],
[-110.984, 24.147],
[-110.72, 23.725],
[-110.281, 23.504],
[-110.061, 22.877],
[-109.446, 23.06],
[-109.424, 23.423],
[-110.259, 24.347],
[-110.347, 24.187],
[-110.61, 24.327],
[-110.588, 24.707],
[-111.445, 26.215],
[-112.654, 27.684],
[-112.961, 28.053],
[-114.192, 27.975]]],
'type': 'Polygon'},
'properties': {},
'type': 'Feature'}
# Generate the bounding box of the geometry used for the spatial query
bbox = bounds(feature['geometry'])
# Instantiate Usgs class and login
usgs = Usgs()
usgs.login()
collection = usgs.get_collection_name(8)
# Query for records of Landsat 8 intersecting with the spatio-temporal window
meta_list = usgs.search(collection=collection, bbox=bbox,
begin=datetime(2018, 1, 1), end=datetime(2018, 2, 28),
max_cloud_cover=30)
print(len(meta_list))
# Printed:
# 52
Because we queried the data using the extent, it is highly probable that some scenes do not intersect with the initial geometry but only with its extent. We can therefore filter the list of scene metadata by testing for intersection between the scene bounds and the geometry. This is done using shapely
region_geom = shape(feature['geometry'])
meta_list = [x for x in meta_list if
shape(geom_from_metadata(x)).intersects(region_geom)]
print(len(meta_list))
# printed:
# 27
The amount of element has reduced by half compared to the total API hits and we are now sure to have retained only scenes that actually intersect with the initial geometry.
We can now proceed to preparing the scene list for placing the order to espa
scene_list = [x['displayId'] for x in meta_list]
espa = Espa()
order = espa.order(scene_list,
products=['pixel_qa', 'sr_ndmi'])
We can then track the status of the order and eventually download it once processing is completed
print(order.status)
# printed:
# ordered