Downloading Landsat Satellite Images with Python | By Conor O’Sullivan | May 2023

Simplify Landsat scene downloads with the landsatxplore Python package

Landsat satellites are the most commonly used source of Earth observation data. They have been providing high-quality images of the planet’s surface for more than four decades. However, manually downloading these images can be tedious! Fortunately, with the landsatxplore package, you can easily download and process Landsat scenes with just a few lines of code.

We will explore the landsatxplore package and walk through the steps to download Landsat satellite images using Python. It includes:

• Setting up API access with a USGS account
• Search and filter Landsat scenes
• Download and run scenes in Python

Say goodbye to manual downloads and hello to an automated, efficient workflow!

Step 1: Register with USGS

To get started, you need to set up a USGS account. This is the same account you use to download scenes with EarthExplorer. keep yours username and password As we will need later.

Once you’re registered, you can use the USGS M2M API. However, it will take some work to set it up. Instead, we’ll use the launchxplore package. It abstracts most of the technical details for you.

Step 2: Install launchxplore

Follow the instructions on GitRepo.

Step 3: Test the API connection

Use the code below to verify that everything is set up correctly. You have to change username and password with the details you used to register your USGS account.

from landsatxplore.api import API
# Your USGS  credentials
username = "XXXXXXXXXXXX"
password = "XXXXXXXXXXXX"
# Initialize a new API instance
api = API(username, password)
# Perform a request
response = api.request(endpoint="dataset-catalogs")
print(response)

The result of the response should look like this:

‘EE’: ‘EarthExplorer’, ‘GV’: ‘GloVis’, ‘HDDS’: ‘HDDS Explorer’

This is a dataset that is available through an API. For this tutorial, we will only consider the EarthExplorer database.

Before we move on to downloading scenes with the API, we’ll do a manual search with EarthExplorer. This is so we can compare the results to what we see using Python. If you are not familiar with the EarthExplorer portal, this tutorial may help.

We are looking for scenes with the following criteria:

• They must contain a point at a given latitude and longitude. This point falls on Bull Island in Dublin.
• taken between 01/01/2020 that 31.12.2022
• Maximum cloud cover 50%
• Part of a Level 2 Landsat 8 or 9 collection

You can see the search results below. We note a few things to compare with our Python search:

• They are there 54 scenes that matches the search criteria.
• They are there 2 tiles which contains a point on Bull Island. They have path and row values (206, 023) and (205023).
• The ID of the first scene is LC08_L2SP_206023_20221118_20221128_02_T1. If you’re wondering what this ID means, see the Landsat naming convention.

Search for scenes

Let’s do an equivalent search using landsatxplore. We do this using the search function below:

• data set — Determines which satellite we want scenes for. The value we use is the dataset ID for Landsat 8 and 9. See GitRepo for IDs for Landsat 5 and 7.
• The latitude and latitude Give the same point on Bull Island. Except that we converted the coordinates to decimals.
• The start date, Completion date and max_cloud_cover Also same as before.

# Search for Landsat TM scenes
scenes = api.search(
dataset='landsat_ot_c2_l2',
latitude=53.36305556,
longitude=-6.15583333,
start_date='2020-01-01',
end_date='2022-12-31',
max_cloud_cover=50
)
# log out
api.logout()

The search result will return information in a JSON dictionary. We convert this to a Pandas DataFrame (line 4), where each row represents a unique scene. A lot of metadata is back! So we filter out what is necessary for this application (line 5). Finally, we order this purchase_date – The date the scene was captured by Landsat.

import pandas as pd
# Create a DataFrame from the scenes
df_scenes = pd.DataFrame(scenes)
df_scenes = df_scenes[['display_id','wrs_path', 'wrs_row','satellite','cloud_cover','acquisition_date']]
df_scenes.sort_values('acquisition_date', ascending=False, inplace=True)

You can see an image of this data set below. If we compare this to our search using EarthExplorer, we can be sure that the results are the same. This data set has 54 rows And there are two unique ones wrs_path and wrs_row couples – (206, 23) and (205,23). Პirveli display_id It’s also the same as what we saw earlier.

If we wanted, we could filter this database. We can use it satellite column to select only images from Landsat 8 or 9. Besides, cloud_cover The column gives the percentage of the image covered by clouds. When you are satisfied with the final list of scenes, you can proceed to download them.

Data is being downloaded

Below, we have the code used to download the Landsat scene. We use the EarthExplorer function (line 5). This is initialized the same way as before – using your USGS credentials. To download the scene, we need to use it display_id and specify the output directory (line 12). we use display_id of the first scene we mentioned above (line 8).

from landsatxplore.earthexplorer import EarthExplorer
import os
# Initialize the API
ee = EarthExplorer(username, password)
# Select the first scene
ID = 'LC08_L2SP_206023_20221118_20221128_02_T1'
# Download the scene 
try: 
ee.download(ID, output_dir='./data')
print(' succesful'.format(ID))
# Additional error handling
except:
if os.path.isfile('./data/.tar'.format(ID)):
print(' error but file exists'.format(ID))
else:
print(' error'.format(ID))
ee.logout()

You may notice additional error handling above. This is due to a problem with the package. In some cases, the scenes will download correctly, but the error will still occur. An additional bug double checks that the scene file exists.

Working with data

The scene will be downloaded as a tar file. will be the file name display_id followed .tar:

LC08_L2SP_206023_20221118_20221128_02_T1.tar

We can work with this data directly in Python. To begin, we need to open the tarfile (lines 4–6). The new folder name is set to the stage display_id (line 5).

import tarfile
# Extract files from tar archive
tar = tarfile.open('./data/.tar'.format(ID))
tar.extractall('./data/'.format(ID))
tar.close()

You can see the unzipped folder and all available files below. This includes all information available for Lansat-2 level science products. The applications for this data are endless! For example, we visualize this scene using bands of visible light. These are available in the files highlighted below.

We load the blue, green, and red bars (lines 6–8). We arrange these bands (line 11), align them (line 12), and glue them to enhance the contrast (line 15). Finally, we display this image using matplotlib (lines 18–20). You can see this image below.

import tifffile as tiff
import numpy as np
import matplotlib.pyplot as plt
# Load Blue (B2), Green (B3) and Red (B4) bands
B2 = tiff.imread('./data//_SR_B2.TIF'.format(ID, ID))
B3 = tiff.imread('./data//_SR_B3.TIF'.format(ID, ID))
B4 = tiff.imread('./data//_SR_B4.TIF'.format(ID, ID))
# Stack and scale bands
RGB = np.dstack((B4, B3, B2))
RGB = np.clip(RGB*0.0000275-0.2, 0, 1)
# Clip to enhance contrast
RGB = np.clip(RGB,0,0.2)/0.2
# Display RGB image
fig, ax = plt.subplots(figsize=(10, 10))
plt.imshow(RGB)
ax.set_axis_off()

Check out this article if you want more details on working with RGB channels of satellite images:

Visualizing RGB channels is just the beginning. After downloading the data, we can perform any remote sensing task – from calculating indices and training models. The best part is that we can do all this without leaving our Python notebook.