How to Download Landsat Imagery: EarthExplorer and Python

Landsat is the longest record of the Earth from space. Its satellites have imaged the planet since the 1970s, and the surface reflectance products of Collection 2 reach back to 1982. The cell is 30 m, coarser than Sentinel-2, but the long history makes Landsat the choice when you need to look back decades, or when you want a second source to confirm what Sentinel-2 shows. Landsat 8 and Landsat 9 together pass over a place about every 8 days.

This guide shows two routes to the same data. Route A uses USGS EarthExplorer, a website where you click. Route B uses a short Python script that reads only your area.

Decide what to download

The product. Collection 2 Level-2 is surface reflectance, corrected for the atmosphere, and it also has surface temperature. Choose Level-2 for vegetation, water and land cover. Level-1 is top-of-atmosphere data that still contains the haze.

The scene. Landsat scenes follow a fixed grid called the Worldwide Reference System 2. Each scene has a path and a row, about 185 km wide. The name of a scene carries all of this. For example LC09_L2SP_148039_20230224_20230307_02_T1 reads:

• LC09: the sensor and satellite, Landsat 9 OLI and TIRS. Landsat 8 is LC08.
• L2SP: Level-2 with both surface reflectance and surface temperature.
• 148039: path 148, row 039.
• 20230224: the date of the pass.
• 20230307: the date it was processed.
• 02: Collection 2.
• T1: Tier 1, the best quality class. Use Tier 1 for time series.

The cloud. A scene can have 30 percent cloud over its whole 185 km and be clear over your small area. Always preview.

Route A: USGS EarthExplorer

EarthExplorer is the USGS portal for Landsat and many other datasets. Searching is open, and downloading needs a free USGS account. The steps below are partly checked on the live site, as noted after them.

1. Open earthexplorer.usgs.gov. Searching needs no account. Downloading needs a free USGS account, and the Login link is at the top right.
2. The page opens on the Search Criteria tab, headed 1. Enter Search Criteria. There are two ways to set the area: the Geocoder, where you pick a geocoding method and type a place name, and GeoJSON/KML/Shapefile Upload. You can also click or draw on the map, and set a date range.
3. Open the Data Sets tab, headed 2. Select Your Data Set(s). The data sets are in a list of categories, among them Digital Elevation, Landsat and Radar, and a plus sign next to a category opens it. Open Landsat, then Landsat Collection 2 Level-2, and tick Landsat 8-9 Operational Land Imager and Thermal Infrared Sensor Collection 2 Level-2. The same group has Landsat 7 and Landsat 4-5 for earlier years. A small cloud icon beside a data set means that it supports a cloud cover filter.
4. Open the Additional Criteria tab. Set a maximum land cloud cover, such as 20 percent.
5. Click the Results tab. Each result has a footprint button and a browse image so that you can preview it. Check that cloud is not sitting on your area.
6. Click the Download Options icon of the scene you want. USGS lets you download single bands, so you do not have to take the whole bundle. For vegetation work you need the red band (SR_B4) and the near infrared band (SR_B5), plus the quality band QA_PIXEL. For other indices, take the bands in the table below.

The tab names, the headings, the Login link and the data set names in steps 1 to 4 were read on the live site. The Results tab and the download options, steps 5 and 6, follow the USGS description of the portal. The screens change from time to time, so check each label against what you see.

Each band arrives as a GeoTIFF with the scene name in its file name. The values are whole numbers and need scaling, as the next section explains.

The bands of Landsat 8 and 9 surface reflectance:

Route B: Python, no account

Route B asks a catalogue for the scenes, then reads only the pixels inside your area from a file that sits online. The catalogue is the Microsoft Planetary Computer, which is free to read. It follows a standard called STAC (SpatioTemporal Asset Catalog), the same one Sentinel-2 guides use. The files are Cloud Optimized GeoTIFFs, so a program can read a window without downloading the rest.

The files are protected by a short-lived access key called a SAS token. You do not need an account. The script asks the Planetary Computer for a token and attaches it to each file address. A token lasts for a limited time, so ask for a new one each session.

You need Python with requests, numpy, rasterio and shapely.

import requests

# Step 1: describe what you want. The box is west, south, east, north in longitude and latitude.
study_area = [75.40, 30.74, 75.55, 30.84]        # a box of about 14 km by 11 km near Jagraon, Punjab

search = {
    "collections": ["landsat-c2-l2"],             # Collection 2, Level-2
    "bbox": study_area,
    "datetime": "2023-02-01T00:00:00Z/2023-02-28T23:59:59Z",
    "query": {"eo:cloud_cover": {"lt": 20}, "platform": {"in": ["landsat-8", "landsat-9"]}},
    "limit": 20,
}

# Step 2: ask the catalogue. It answers with a list of scenes.
url = "https://planetarycomputer.microsoft.com/api/stac/v1/search"
scenes = requests.post(url, json=search, timeout=60).json()["features"]

# Step 3: the catalogue returns every scene that touches your box. Keep only the scenes that cover all of it.
from shapely.geometry import shape, box
area = box(*study_area)
covering = [s for s in scenes if shape(s["geometry"]).contains(area)]
print(len(scenes), "scenes touch the box,", len(covering), "cover all of it")

# Step 4: print the date, the cloud percentage, the path and row, and the name of every scene that covers it
for scene in covering:
    p = scene["properties"]
    print(p["datetime"][:10], round(p["eo:cloud_cover"], 1), p["landsat:wrs_path"] + "/" + p["landsat:wrs_row"], scene["id"])

3 scenes touch the box, 2 cover all of it
2023-02-24 4.6 148/039 LC09_L2SP_148039_20230224_02_T1
2023-02-16 18.1 148/039 LC08_L2SP_148039_20230216_02_T1

Step 3 matters. A scene is returned when its footprint touches your box, and touching is not the same as covering. A scene from the next row up can overlap your box by a thin strip, and if you chose it only because it has little cloud, you would read a window that is almost all empty. The shapely test keeps the scenes whose footprint contains the whole box. The first line shows how many were dropped.

Pick the clearest scene and look at what it offers.

# Step 1: take the scene with the lowest cloud cover among those that cover your box
scene = min(covering, key=lambda s: s["properties"]["eo:cloud_cover"])
print("chosen:", scene["id"], "| cloud:", round(scene["properties"]["eo:cloud_cover"], 1), "percent")

# Step 2: every file is an asset. List the ones that matter here.
assets = scene["assets"]
print([name for name in assets if name in ("red", "nir08", "qa_pixel", "swir16")])

# Step 3: the scale and offset needed for reflectance sit in the metadata of each band
print(assets["red"]["raster:bands"][0])

chosen: LC09_L2SP_148039_20230224_02_T1 | cloud: 4.6 percent
['red', 'nir08', 'swir16', 'qa_pixel']
{'scale': 2.75e-05, 'nodata': 0, 'offset': -0.2, 'data_type': 'uint16', 'spatial_resolution': 30}

The asset names are descriptive, red, green, blue, nir08, swir16 and swir22, so you do not need to remember that red is band 4 on Landsat 8 and 9. Their numbers differ on older Landsat satellites, and the names stay the same.

Get an access key and read your area

import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

# Step 1: ask for a token. It is a short piece of text that unlocks the files for a limited time.
token = requests.get("https://planetarycomputer.microsoft.com/api/sas/v1/token/landsat-c2-l2", timeout=60).json()["token"]

def read_window(asset, bbox):
    # Step 2: open the file online with the token added to its address. Nothing is downloaded yet.
    with rasterio.open(asset["href"] + "?" + token) as src:
        # Step 3: your box is in longitude and latitude, the file is in metres, so convert the box first
        left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *bbox)
        window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()

        # Step 4: read just that window. Only these pixels travel over the network.
        values = src.read(1, window=window)
        transform = src.window_transform(window)
        crs = src.crs
    return values, transform, crs

red_raw, transform, crs = read_window(assets["red"], study_area)
nir_raw, _, _ = read_window(assets["nir08"], study_area)
print("window:", red_raw.shape[1], "x", red_raw.shape[0], "cells | CRS:", crs.to_string(), "| data type:", red_raw.dtype)
print("cells with no data:", round(float((red_raw == 0).mean()) * 100, 1), "percent")
print("raw red values: from", int(red_raw[red_raw > 0].min()), "to", int(red_raw.max()))

window: 480 x 371 cells | CRS: EPSG:32643 | data type: uint16
cells with no data: 0.0 percent
raw red values: from 7898 to 17674

The raw values are whole numbers of several thousand. They become reflectance with a scale and an offset. Landsat Collection 2 uses a scale of 0.0000275 and an offset of -0.2. These are not the Sentinel-2 numbers, so never reuse one set for the other. Read them from the metadata.

# Step 1: whole numbers to reflectance. The value 0 means no data.
def to_reflectance(raw, asset):
    info = asset["raster:bands"][0]
    return np.where(raw == 0, np.nan, raw * info["scale"] + info["offset"]).astype("float32")

red = to_reflectance(red_raw, assets["red"])
nir = to_reflectance(nir_raw, assets["nir08"])
print("red reflectance: from", round(float(np.nanmin(red)), 3), "to", round(float(np.nanmax(red)), 3))

# Step 2: a quick sanity check. NDVI of winter crops in February should be high, and never above 1.
ndvi = (nir - red) / (nir + red)
print("NDVI: median", round(float(np.nanmedian(ndvi)), 2), "| highest", round(float(np.nanmax(ndvi)), 2))

red reflectance: from 0.017 to 0.286
NDVI: median 0.75 | highest 0.89

If the median NDVI of a green area is near 0.7 and no value is above 1, the scaling is right. If it is wildly different, check the scale and the offset first.

Mask cloud with the quality band

Landsat Level-2 comes with a quality band called qa_pixel. Each cell holds a number whose individual bits say what the cell is. Bit 3 is cloud, bit 4 is cloud shadow, and bit 2 is thin cirrus. A bit is one binary digit of the number, and the test below asks whether it is switched on.

# Step 1: read the quality band for the same window
qa, _, _ = read_window(assets["qa_pixel"], study_area)

# Step 2: test the bits. (qa >> 3) & 1 shifts bit 3 to the right edge and keeps only that bit.
cloud = ((qa >> 3) & 1) == 1
shadow = ((qa >> 4) & 1) == 1
cirrus = ((qa >> 2) & 1) == 1
bad = cloud | shadow | cirrus
print("cloud, shadow or cirrus:", round(float(bad.mean()) * 100, 1), "percent of the window")

# Step 3: hide those cells in the red and near infrared bands
red_clear = np.where(bad, np.nan, red)
nir_clear = np.where(bad, np.nan, nir)

cloud, shadow or cirrus: 0.0 percent of the window

Save what you read

# Step 1: write the two bands to one GeoTIFF with the same position and CRS as the window
with rasterio.open(
    "landsat_red_nir.tif", "w", driver="GTiff",
    height=red_clear.shape[0], width=red_clear.shape[1], count=2,
    dtype="float32", crs=crs, transform=transform, nodata=np.nan, compress="deflate",
) as dst:
    dst.write(red_clear, 1)
    dst.write(nir_clear, 2)
    dst.set_band_description(1, "red SR_B4")
    dst.set_band_description(2, "nir SR_B5")

# Step 2: read it back to check
with rasterio.open("landsat_red_nir.tif") as check:
    print("saved:", check.count, "bands |", check.width, "x", check.height, "cells | cell size", check.res[0], "m")

saved: 2 bands | 480 x 371 cells | cell size 30.0 m

Landsat or Sentinel-2?

The two can be used together. Landsat's 30 m cell is nine times the area of a Sentinel-2 cell, so a field of a hectare is only about eleven Landsat cells.

Common mistakes

• Downloading Level-1 because it is listed first. Choose Level-2 for surface reflectance.
• Reusing the Sentinel-2 scale and offset on Landsat, or the reverse.
• Forgetting that the band numbers change between Landsat generations. The asset names red and nir08 avoid this.
• Judging a scene by the cloud percentage of the whole 185 km and skipping the preview.
• Choosing a scene that only touches your box. Keep the ones that cover all of it, and check the share of cells with no data.
• Taking Tier 2 scenes into a time series. Their position can be off by more than a cell, so they may not line up with Tier 1 scenes.
• Using a token from yesterday. Ask for a fresh one each session.

Key takeaways