Demo: Landsat background and data download

Demo: Landsat background and data download#

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean, Eric Gagliano, Quinn Brencher

Introduction#

In this demo, we’re going to download the Landsat 8 data that we’ll need for this module. While the data is being downloaded, we’ll briefly go over the Landsat program and the Landsat 8 data we’ll be working with.

The Landsat 8 mission#

The joint NASA and USGS Landsat program has been acquiring multi-spectral satellite imagery since 1972.

landsat timeline

Timeline of the Landsat program, beginning with the launch of Landsat 1 in 1972. Landsat Next, consisting of a trio of satellite observatories, is expected to launch in late 2030. As the tenth Landsat mission, it will continue the legacy of the Landsat program.

Orbit and data collection#

Landsat 8 orbits the the Earth in a sun-synchronous, near-polar orbit, at an altitude of 705 km (438 mi), inclined at 98.2 degrees, and completes one Earth orbit every 99 minutes. The satellite has a 16-day repeat cycle with an equatorial crossing time: 10:00 a.m. +/- 15 minutes.

Landsat 8

Landsat 8 aquires about 740 scenes a day on the Worldwide Reference System-2 (WRS-2) path/row system, with a swath overlap (or sidelap) varying from 7 percent at the equator to a maximum of approximately 85 percent at extreme latitudes. A Landsat 8 scene size is 185 km x 180 km (114 mi x 112 mi).

Instruments and bands#

  • Landsat-8 has two instruments…

    • Operational Land Imager (OLI)

    • Thermal Infrared Sensor (TIRS) LS bands chart The OLI collects data for two new bands, a coastal/aerosol band (band 1) and a cirrus band (band 9), as well as the heritage Landsat multispectral bands. Additionally, the bandwidth has been refined for six of the heritage bands. The Thermal Instrument (TIRS) carries two additional thermal infrared bands. Note: atmospheric transmission values for this graphic were calculated using MODTRAN for a summertime mid-latitude hazy atmosphere (circa 5 km visibility). Graphic via USGS

Image resolution - Ground Sample Distance (GSD)#

  • Panchromatic (PAN) band (band number 8) has 15 m ground sample distance (GSD)

  • Multispectral (MS) bands are 30 m GSD

  • Thermal IR are 100 m GSD, but are often oversampled to match MS bands

Landsat-8 spatial resolution

Dynamic range#

  • LS8 OLI provides 12-bit dynamic range, which improves characterization and signal to noise ratio. 12 bits means 2^12 (or 4096) unique combinations to represent brightness in the image.

  • We don’t have a convenient mechanism to store 12-bit data, so the LS8 images are stored as 16-bit unsigned integers. The initial values (spanning 0-4095) are scaled across 55000 of the total 65535 brightness levels in the 16-bit images.

Data products#

The standard data products are “Level 1” images, which are radiometrically corrected and orthorectified (terrain corrected) in the approprate UTM projection

For more sophisticated analysis, you typically want to use the higher-level “Level 2”, calibrated/corrected surface reflectance products, often considered “Analysis Ready Data (ARD)”

You’ll often see Digital Number (DN), top of atmosphere (TOA) and surface reflectance (SR)

  • “Top-of-atmosphere reflectance (or TOA reflectance) is the reflectance measured by a space-based sensor flying higher than the earth’s atmosphere. These reflectance values will include contributions from clouds and atmospheric aerosols and gases.”

  • “Surface reflectance (SR) is the amount of light reflected by the surface of the Earth. It is a ratio of surface radiance to surface irradiance, and as such is unitless, with values between 0 and 1.”

    • “Surface reflectance improves comparison between multiple images over the same region by accounting for atmospheric effects such as aerosol scattering and thin clouds, which can help in the detection and characterization of Earth surface change. “

    More information here.

Data availability#

USGS/NASA hosts the official Landsat products on EarthExplorer. This option is great for one-off interactive data searches, but can be clunky and requires a lot of manual steps

Commercial cloud providers now mirror the entire USGS archive, and provide a much more efficient API (application programming interface) to access the data. This is especially important when you need to access 100s-1000s of images.

Finding imagery#

The process of identifying and downloading data has evolved considerably throughout the Landsat missions, and modern approaches use on-demand access to cloud-hosted archives, often without local downloads. We will opt for a programmatic approach, using pystac/pystac-client to query a SpatioTemporal Asset Catalog (STAC)

Other approaches include…

  • Interactive, manual approach. EarthExplorer for visual queries, and either download directly, or save image ID and use another approach.

  • Other approaches

For now, we’ll use the following code to download some pre-identified images for WA state to the hub. Note that we could do all of this on the fly, but for experimentation and development, having a local copy of sample images will speed up reading.

Download#

import pystac
from pystac_client import Client
import planetary_computer
import odc.stac
import matplotlib.pyplot as plt
import os
from pathlib import Path
import urllib
import rioxarray

imgdir = f'{Path.home()}/gda_demo_data/LS8_data'

if not os.path.exists(imgdir):
    os.makedirs(imgdir)

#Collection and bands (specific bands for L2 surface reflectance/temperature products)
collection = 'landsat-c2-l2'
asset_id_list = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'ST_B10', 'reduced_resolution_browse']
asset_id_list = ['coastal', 'blue', 'green', 'red', 'nir08', 'swir16', 'lwir11', 'rendered_preview']

#Bounding box
bbox = [-121.9,46.7,-121.6,47.0]

#Date range
dt = '2018-08-17/2018-12-26'

catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",modifier=planetary_computer.sign_inplace)

results = catalog.search(
    collections=[collection],
    bbox=bbox,
    datetime=dt,
    query={
        "eo:cloud_cover": {"lt": 50},
        "landsat:wrs_path": {"eq": "046"}, # Path 46, Row 27 is the tile that contains Seattle and Mt. Rainier!
        "landsat:wrs_row": {"eq": "027"},
    },
)

# Check how many items were returned
items = list(results.items())
print(f"Returned {len(items)} Items")

Returned 7 Items

# let's just take the first and last set of images so we get an image in august and an image from december...
items = [items[0],items[6]]
items

[<Item id=LC08_L2SP_046027_20181224_02_T1>,
 <Item id=LC08_L2SP_046027_20180818_02_T1>]

items[0]

for item in items:
    for asset_id in asset_id_list:
        asset = item.assets[asset_id]
        print(asset.href)

        if asset_id == 'rendered_preview':
            out_fn = item.id + "_rendered_preview.png"
        else:
            out_fn = asset.href.split('?')[0].split('/')[-1]
            
        out_fp = os.path.join(imgdir, out_fn)

        # Check to see if file already exists
        if not os.path.exists(out_fp):
            print("Saving:", out_fp)
            # Download the file
            urllib.request.urlretrieve(asset.href, out_fp)

https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2018/046/027/LC08_L2SP_046027_20181224_20200829_02_T1/LC08_L2SP_046027_20181224_20200829_02_T1_SR_B1.TIF?st=2026-01-25T22%3A17%3A34Z&se=2026-01-26T23%3A02%3A34Z&sp=rl&sv=2025-07-05&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2026-01-26T06%3A46%3A59Z&ske=2026-02-02T06%3A46%3A59Z&sks=b&skv=2025-07-05&sig=5DWlK65HZBjyTlT6xdtasG5EwDRZzeM%2Bw04gH%2Bspbz4%3D
Saving: /home/jovyan/gda_demo_data/LS8_data/LC08_L2SP_046027_20181224_20200829_02_T1_SR_B1.TIF
https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2018/046/027/LC08_L2SP_046027_20181224_20200829_02_T1/LC08_L2SP_046027_20181224_20200829_02_T1_SR_B2.TIF?st=2026-01-25T22%3A17%3A34Z&se=2026-01-26T23%3A02%3A34Z&sp=rl&sv=2025-07-05&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2026-01-26T06%3A46%3A59Z&ske=2026-02-02T06%3A46%3A59Z&sks=b&skv=2025-07-05&sig=5DWlK65HZBjyTlT6xdtasG5EwDRZzeM%2Bw04gH%2Bspbz4%3D
Saving: /home/jovyan/gda_demo_data/LS8_data/LC08_L2SP_046027_20181224_20200829_02_T1_SR_B2.TIF

Inspect the downloaded data#

Do a quick ls -lh on the local image directory#

  • Note relative file sizes for the different bands of each image

    • Revisit the chart above - does this make sense for resolution of these bands?

!ls -lh $imgdir

Use the command line gdalinfo to get some basic information about one of the files#

  • Anything unexpected?

  • How many bands does each file contain? Why do you think this is?

    • Is the rendered preview different? Why?

  • What are “overviews”??

august_band1_filename = "LC08_L2SP_046027_20180818_20200831_02_T1_SR_B1.TIF"
august_band10_filename = "LC08_L2SP_046027_20180818_20200831_02_T1_ST_B10.TIF"
august_preview_filename = "LC08_L2SP_046027_20180818_02_T1_rendered_preview.png"
december_preview_filename = "LC08_L2SP_046027_20181224_02_T1_rendered_preview.png"

!gdalinfo $imgdir/$august_band1_filename
#!gdalinfo $imgdir/$august_band10_filename
#!gdalinfo $imgdir/$august_preview_filename

Let’s check out the rendered previews#

  • Do the dimensions match the output of gdalinfo? Do they make sense?

  • What do you make of the data type?

  • Do you think this image is georeferenced? (hint: try commenting out the two axis('off') in the plot)

  • Any differences between the two images?

august_preview_da = rioxarray.open_rasterio(f'{imgdir}/{august_preview_filename}')
december_preview_da = rioxarray.open_rasterio(f'{imgdir}/{december_preview_filename}')
august_preview_da

f,ax=plt.subplots(1,2,figsize=(12,6))

august_preview_da.plot.imshow(ax=ax[0])
ax[0].set_title('August 2018')
ax[0].invert_yaxis()
ax[0].axis('off')
ax[0].set_aspect('equal')

december_preview_da.plot.imshow(ax=ax[1])
ax[1].set_title('December 2018')
ax[1].invert_yaxis()
ax[1].axis('off')
ax[1].set_aspect('equal')

f.tight_layout()