Landsat-8 on AWS¶
How can I efficiently find, visualize, and analyze a large set of Landsat-8 imagery on AWS?
We will examine raster images from the Landsat-8 instrument. The Landsat program is the longest-running civilian satellite imagery program, with the first satellite launched in 1972 by the US Geological Survey. Landsat 8 is the latest satellite in this program, and was launched in 2013. Landsat observations are processed into “scenes”, each of which is approximately 183 km x 170 km, with a spatial resolution of 30 meters and a temporal resolution of 16 days. The duration of the landsat program makes it an attractive source of medium-scale imagery for land surface change analyses.
This notebook is a simplified update of a blog post written in October 2018, highlighting the integration of several modern Python libraries: satsearch, intake-stac, geopandas, xarray, dask, holoviz
Finding data on the Cloud¶
Finding geospatial data on the Cloud has been difficult until recent years. Earth scientists are accustomed to going to specific portals to find data, for example NASA’s Earthdata Search, ESA’s Copernicus Hub (https://scihub.copernicus.eu), USGS National Map (https://www.usgs.gov/core-science-systems/national-geospatial-program/national-map). AWS has had a registry of open datasets stored on AWS for many years now https://aws.amazon.com/opendata/public-datasets/. Earth-science specific data is also listed here - https://aws.amazon.com/earth/. But what if you want to search for Landsat8 data over an area of interest? Browsing through lists of files is cumbersome.
An effort is underway to make geospatial data on the Cloud more easy to discover by having a bare-bones web-friendly and search-friendly metadata catalog standard: The SpatioTemporal Asset Catalog (STAC). Once the standard is set, many tools can co-exist to build upon it. For example https://www.element84.com/earth-search/ allows us to programmatically and visually search for data on AWS! Here we will use the satsearch Python library for searching Landsat8 on AWS. Note you could also search for Sentinel2
[1]:
# Suppress library deprecation warnings
import warnings
warnings.filterwarnings('ignore')
[2]:
import satsearch
print(satsearch.__version__)
print(satsearch.config.API_URL)
from satsearch import Search
0.2.3
https://earth-search-legacy.aws.element84.com
[3]:
# NOTE this STAC API endpoint does not currently search the entire catalog
bbox = (-124.71, 45.47, -116.78, 48.93) #(west, south, east, north)
timeRange = '2019-01-01/2020-10-01'
# STAC metadata properties
properties = ['eo:row=027',
'eo:column=047',
'landsat:tier=T1']
results = Search.search(collection='landsat-8-l1',
bbox=bbox,
datetime=timeRange,
property=properties,
sort=['<datetime'],
)
print('%s items' % results.found())
items = results.items()
items.save('subset.geojson')
19 items
[4]:
# Remember that it is easy to load geojson with geopandas!
import geopandas as gpd
gf = gpd.read_file('subset.geojson')
gf.head()
[4]:
| id | collection | eo:gsd | eo:platform | eo:instrument | eo:off_nadir | eo:bands | datetime | eo:sun_azimuth | eo:sun_elevation | eo:cloud_cover | eo:row | eo:column | landsat:product_id | landsat:scene_id | landsat:processing_level | landsat:tier | landsat:revision | eo:epsg | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | LC80470272019096 | landsat-8-l1 | 15 | landsat-8 | OLI_TIRS | 0 | [ { "name": "B1", "common_name": "coastal", "g... | 2019-04-06T19:01:20 | 152.713795 | 46.106183 | 79 | 027 | 047 | LC08_L1TP_047027_20190406_20190422_01_T1 | LC80470272019096LGN00 | L1TP | T1 | 00 | 32610 | POLYGON ((-124.30921 48.51312, -121.80425 48.0... |
| 1 | LC80470272019128 | landsat-8-l1 | 15 | landsat-8 | OLI_TIRS | 0 | [ { "full_width_half_max": 0.02, "center_wavel... | 2019-05-08T19:01:19 | 149.213545 | 56.592420 | 19 | 027 | 047 | LC08_L1TP_047027_20190508_20190521_01_T1 | LC80470272019128LGN00 | L1TP | T1 | 00 | 32610 | POLYGON ((-124.27756 48.51300, -121.77234 48.0... |
| 2 | LC80470272019160 | landsat-8-l1 | 15 | landsat-8 | OLI_TIRS | 0 | [ { "full_width_half_max": 0.02, "center_wavel... | 2019-06-09T19:01:37 | 143.652070 | 61.668260 | 63 | 027 | 047 | LC08_L1TP_047027_20190609_20190619_01_T1 | LC80470272019160LGN00 | L1TP | T1 | 00 | 32610 | POLYGON ((-124.28944 48.51325, -121.78398 48.0... |
| 3 | LC80470272019176 | landsat-8-l1 | 15 | landsat-8 | OLI_TIRS | 0 | [ { "full_width_half_max": 0.02, "center_wavel... | 2019-06-25T19:01:42 | 141.834366 | 61.713605 | 52 | 027 | 047 | LC08_L1TP_047027_20190625_20190705_01_T1 | LC80470272019176LGN00 | L1TP | T1 | 00 | 32610 | POLYGON ((-124.28314 48.51335, -121.77799 48.0... |
| 4 | LC80470272019208 | landsat-8-l1 | 15 | landsat-8 | OLI_TIRS | 0 | [ { "full_width_half_max": 0.02, "center_wavel... | 2019-07-27T19:01:51 | 143.987513 | 57.532989 | 76 | 027 | 047 | LC08_L1TP_047027_20190727_20190801_01_T1 | LC80470272019208LGN00 | L1TP | T1 | 00 | 32610 | POLYGON ((-124.26661 48.51255, -121.76179 48.0... |
[5]:
# Tidy display of band information from the 'eo:bands column'
import ast
import pandas as pd
band_info = pd.DataFrame(ast.literal_eval(gf.iloc[0]['eo:bands']))
band_info
[5]:
| name | common_name | gsd | center_wavelength | full_width_half_max | |
|---|---|---|---|---|---|
| 0 | B1 | coastal | 30 | 0.44 | 0.02 |
| 1 | B2 | blue | 30 | 0.48 | 0.06 |
| 2 | B3 | green | 30 | 0.56 | 0.06 |
| 3 | B4 | red | 30 | 0.65 | 0.04 |
| 4 | B5 | nir | 30 | 0.86 | 0.03 |
| 5 | B6 | swir16 | 30 | 1.60 | 0.08 |
| 6 | B7 | swir22 | 30 | 2.20 | 0.20 |
| 7 | B8 | pan | 15 | 0.59 | 0.18 |
| 8 | B9 | cirrus | 30 | 1.37 | 0.02 |
| 9 | B10 | lwir11 | 100 | 10.90 | 0.80 |
| 10 | B11 | lwir12 | 100 | 12.00 | 1.00 |
[6]:
# Plot search AOI and frames on a map using Holoviz Libraries (more on these later)
import geoviews as gv
import hvplot.pandas
cols = gf.loc[:,('id','eo:row','eo:column','geometry')]
footprints = cols.hvplot(geo=True, line_color='k', hover_cols=['eo:row','eo:column'], alpha=0.1, title='Landsat 8 T1')
tiles = gv.tile_sources.CartoEco.options(width=700, height=500)
labels = gv.tile_sources.StamenLabels.options(level='annotation')
tiles * footprints * labels