Chapter Nine - Landsat Data in Python
In this chapter, you will learn how to work with Landsat multi-band raster data stored in
.tif format in Python using rasterio.
After completing this chapter, you will be able to:
glob()to create a subsetted list of file names within a specified directory on your computer.
- Create a raster stack from a list of
.tiffiles in Python.
- Crop rasters to a desired extent in Python.
- Plot various band combinations using a numpy array in Python.
What You Need
You will need a computer with internet access to complete this lesson and the Cold Springs Fire data.
Or use earthpy
In the previous lessons in this textbook, you learned how to import a multi-band image into Python using rasterio. You then plotted the data as a composite, RGB (and CIR) image using
imshow() and calculated NDVI.
In that case, all bands of the data were stored in a single
.tif file. However, sometimes data are downloaded in individual bands rather than a single file.
In this chapter, you will learn how to work with Landsat data in Python. Each band in a landsat scene is often stored in an individual
.tif file. Thus you will need to grab the bands that you want to work with and then bring them into a
About Landsat Data
At over 40 years, the Landsat series of satellites provides the longest temporal record of moderate resolution multispectral data of the Earth’s surface on a global basis. The Landsat record has remained remarkably unbroken, proving a unique resource to assist a broad range of specialists in managing the world’s food, water, forests, and other natural resources for a growing world population. It is a record unmatched in quality, detail, coverage, and value. Source: USGS
Landsat data are spectral and collected using a platform mounted on a satellite in space that orbits the earth. The spectral bands and associated spatial resolution of the first 9 bands in the Landsat 8 sensor are listed below.
Landsat 8 Bands
|Band||Wavelength range (nanometers)||Spatial Resolution (m)||Spectral Width (nm)||Units||Data Type||Fill Value (no data)||Range||Valid Range||Scale Factor|
|Band 1 - Coastal aerosol||430 - 450||30||2.0||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 2 - Blue||450 - 510||30||6.0||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 3 - Green||530 - 590||30||6.0||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 4 - Red||640 - 670||30||0.03||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 5 - Near Infrared (NIR)||850 - 880||30||3.0||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 6 - SWIR 1||1570 - 1650||30||8.0||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
|Band 7 - SWIR 2||2110 - 2290||30||18||Reflectance||16-bit signed integer (int16)||-9999||-2000 to 16000||0 to 10000||0.0001|
Review the Landsat 8 Surface Reflectance Product Guide for more details.
There are additional collected bands that are not distributed within the Landsat 8 Surface Reflectance Product such as the panchromatic band, which provides a finer resolution, gray scale image of the landscape, and the cirrus cloud band, which is used in the quality assessment process:
|Band||Wavelength range (nanometers)||Spatial Resolution (m)||Spectral Width (nm)|
|Band 8 - Panchromatic||500 - 680||15||18|
|Band 9 - Cirrus||1360 - 1380||30||2.0|
Understand Landsat Data
When working with landsat, it is important to understand both the metadata and the file naming convention. The metadata tell you how the data were processed, where the data are from and how they are structured.
The file names, tell you what sensor collected the data, the date the data were collected, and more.
Landsat File Naming Convention
Landsat and many other satellite remote sensing data is named in a way that tells you a about:
- When the data were collected and processed
- What sensor was used to collect the data
- What satellite was used to collect the data.
Here you will learn a few key components of the landsat 8 collection file name. The first scene that you work with below is named:
First, we have LC08
- L: Landsat Sensor
- C: OLI / TIRS combined platform
08: Landsat 8 (not 7)
- 034032: The next 6 digits represent the path and row of the scene. This identifies the spatial coverage of the scene
Finally, you have a date. In your case as follows:
- 20160723: representing the year, month and day that the data were collected.
The second part of the file name above tells you more about when the data were last processed. You can read more about this naming convention using the link below.
As you work with these data, it is good to double check that you are working with the sensor (Landsat 8) and the time period that you intend. Having this information in the file name makes it easier to keep track of this as you process your data.
import os from glob import glob # File manipulation import matplotlib.pyplot as plt import numpy as np import geopandas as gpd import rasterio as rio import earthpy as et import earthpy.spatial as es import earthpy.plot as ep # Download data and set working directory data = et.data.get_data('cold-springs-fire') os.chdir(os.path.join(et.io.HOME, 'earth-analytics'))
Downloading from https://ndownloader.figshare.com/files/10960109 Extracted output to /root/earth-analytics/data/cold-springs-fire/.
# Get list of all pre-cropped data and sort the data path = os.path.join("data", "cold-springs-fire", "landsat_collect", "LC080340322016072301T1-SC20180214145802", "crop") all_landsat_post_bands = glob(path + "/*band*.tif") all_landsat_post_bands.sort()
# Create an output array of all the landsat data stacked landsat_post_fire_path = os.path.join("data", "cold-springs-fire", "outputs", "landsat_post_fire.tif") # This will create a new stacked raster with all bands land_stack, land_meta = es.stack(all_landsat_post_bands, landsat_post_fire_path)
Open The New Raster Stack
Once you have stacked your data, you can import it and work with it as you need to!
with rio.open(landsat_post_fire_path) as src: landsat_post_fire = src.read()
Plot RGB image
Just like you did with NAIP data, you can plot 3 band color composite images with Landsat too. Below you will plot an RGB image using landsat. Refer to the landsat bands in the table at the top of this page to figure out the red, green and blue bands. Or read the ESRI landsat 8 band combinations post.
To make plotting less code intensive we have created a plot_rgb function that allows you to quickly make a 3 band raster plot. To use it simply provide
- the numpy array containing the bands that you wish to plot.
IMPORTANT: this array should be in rasterio band order (bands first).
- The bands that you wish to plot in the array.
Optionally you can chose to provide a title for the plot, and the figure size if you’d like.
ep.plot_rgb(landsat_post_fire, rgb=[3, 2, 1], title="RGB Composite Image\n Post Fire Landsat Data") plt.show()
Notice that the image above looks dark. You can stretch the image as you did with the NAIP data, too. Below you use the stretch argument built into the earthpy
plot_rgb() function. The
str_clip argument allows you to specify how much of the tails of the data that you want to clip off. The larger the number, the most the data will be stretched or brightened.
ep.plot_rgb(landsat_post_fire, rgb=[3, 2, 1], title="Landsat RGB Image\n Linear Stretch Applied", stretch=True, str_clip=1) plt.show()
# Adjust the amount of linear stretch to futher brighten the image ep.plot_rgb(landsat_post_fire, rgb=[3, 2, 1], title="Landsat RGB Image\n Linear Stretch Applied", stretch=True, str_clip=4) plt.show()
Raster Pixel Histograms
You can create a histogram to view the distribution of pixel values in the rgb bands plotted above.
# Plot all band histograms using earthpy band_titles = ["Band 1", "Blue", "Green", "Red", "NIR", "Band 6", "Band7"] ep.hist(landsat_post_fire, title=band_titles) plt.show()
Now you’ve created a red, green blue color composite image. Remember red green and blue are colors that your eye can see.
Next, create a color infrared image (CIR) using landsat bands: 4,3,2.
ep.plot_rgb(landsat_post_fire, rgb=[4, 3, 2], title="CIR Landsat Image Pre-Cold Springs Fire", figsize=(10, 10)) plt.show()
Data Tip: Landsat 8 Pre Collections Data
If you are working with Landsat data downloaded pre USGS collections, your data may be formatted and named slightly differently than the example shown on this page. Below is an explanation of the legacy Landsat 8 naming convention.
|Sensor||Sensor||Satellite||WRS path||WRS row|
|Landsat||OLI & TIRS||Landsat 8||path = 034||row = 032||year = 2016||Month = 7, day = 7||Processing Date: 2017-02-21 (feb 21 2017)||Archive (second version): 01|
- L: Landsat
- X: Sensor
- C = OLI & TIRS O = OLI only T = IRS only E = ETM+ T = TM M = MSS
- S Satelite
- YYYY = Year
- DDD = Julian DAY of the year
- GSI - Ground station ID
- VV = Archive Version
We won’t spend a lot of time on Julian days. The julian day used to be used in Landsat pre collections file naming. However recently they have switched to a normal year-month-date format
See this link that provide tables to help you convert julian days to actual date.