Earth analytics
earth-analytics Home

Lesson 6. Landsat Remote Sensing tif Files in R

Leah Wasser

Learning Objectives

After completing this tutorial, you will be able to:

  • Use list.files() to create a subsetted list of file names within a specified directory on your computer.
  • Create a raster stack from a list of .tif files in R.
  • Plot various band combinations using a rasterstack in R with plotRGB().

What You Need

You will need a computer with internet access to complete this lesson and the data for week 7 of the course.

Download Week 7 - 9 Data (~500 MB)

In the previous lesson, you learned how to import a multi-band image into R using the stack() function. You then plotted the data as a composite, RGB (and CIR) image using plotRGB(). However, sometimes data are downloaded in individual bands rather than a composite raster stack.

In this lesson you will learn how to work with Landsat data in R. In this case, your data are downloaded in .tif format with each .tif file representing a single band rather than a stack of bands.

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 40 year timeline source: USGS.
The 40 year history of Landsat missions. Source: USGS - USGS Landsat

Landsat data is a multispectral dataset collected from space. The multispectral bands and associated spatial resolution of the first 9 bands in the Landsat 8 sensor are listed below.

Landsat 8 Bands

BandWavelength range (nanometers)Spatial Resolution (m)Spectral Width (nm)
Band 1 - Coastal aerosol430 - 450302.0
Band 2 - Blue450 - 510306.0
Band 3 - Green530 - 590306.0
Band 4 - Red640 - 670300.03
Band 5 - Near Infrared (NIR)850 - 880303.0
Band 6 - SWIR 11570 - 1650308.0
Band 7 - SWIR 22110 - 22903018
Band 8 - Panchromatic500 - 6801518
Band 9 - Cirrus1360 - 1380302.0

Get to Know Landsat 8 Filenames

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 file naming convention
Landsat file names Source: USGS Landsat - <a href="https://landsat.usgs.gov/what-are-naming-conventions-landsat-scene-identifiers

Let’s have a look at one of the files and use the image above to guide us through understanding the file name.

File: LC80340322016205LGN00_sr_band1_crop.tif

SensorSensorSatelliteWRS pathWRS row    
LC80340322016205LGN00
LandsatOLI & TIRSLandsat 8path = 034row = 032year = 2016Julian day= 205Ground station: LGNArchive (first version): 00
  • L: Landsat
  • X: Sensor
    • C = OLI & TIRS O = OLI only T = IRS only E = ETM+ T = TM M = MSS
  • S Satelite
  • PPP
  • RRR
  • YYYY = Year
  • DDD = Julian DAY of the year
  • GSI - Ground station ID
  • VV = Archive Version

More here breaking down the file name.

Julian Day

In class, we won’t spend a lot of time on Julian days. For the purpose of working with Landsat and MODIS data, what you need to know is that the calendar year Julian day represents the numeric day of the year. So Jan 1 = day 1. Feb 1 = day 32. And so on.

There are several links at the bottom of this page that provide tables that help you convert Julian days to actual date.

Landsat tif Files in R

Next, you can open the Landsat data in R.

# load spatial packages
library(raster)
library(rgdal)
library(rgeos)

# turn off factors
options(stringsAsFactors = FALSE)

If you look at the Landsat directory for the week_07 data, you will see that each of the individual bands is stored individually as a GeoTIFF rather than being stored as a stacked or layered, multi-band raster.

Why would they store the data this way?

Conventionally Landsat was stored in a file format called HDF - hierarchical data format. However that format, while extremely efficient, is a bit more challenging to work with. In recent years USGS has started to make each band of a Landsat scene available as a .tif file. This makes it a bit easier to use across many different programs and platforms.

You have already been working with the geotiff file format in this class! You will thus use many of the same functions you used previously, to work with Landsat.

Get List of Files

To begin, let’s explore your file directory in R, You can use list.files() to grab a list of all files within any directory on your computer.

# get list of all tifs
list.files("data/week-07/landsat/LC80340322016205-SC20170127160728/crop")
##  [1] "LC80340322016205LGN00_bqa_crop.tif"        
##  [2] "LC80340322016205LGN00_cfmask_conf_crop.tif"
##  [3] "LC80340322016205LGN00_cfmask_crop.tif"     
##  [4] "LC80340322016205LGN00_sr_band1_crop.tif"   
##  [5] "LC80340322016205LGN00_sr_band2_crop.tif"   
##  [6] "LC80340322016205LGN00_sr_band3_crop.tif"   
##  [7] "LC80340322016205LGN00_sr_band4_crop.tif"   
##  [8] "LC80340322016205LGN00_sr_band5_crop.tif"   
##  [9] "LC80340322016205LGN00_sr_band6_crop.tif"   
## [10] "LC80340322016205LGN00_sr_band7_crop.tif"   
## [11] "LC80340322016205LGN00_sr_cloud_crop.tif"   
## [12] "LC80340322016205LGN00_sr_ipflag_crop.tif"

You can also use list.files() with the pattern argument. This allows you to specify a particular pattern that further subsets your data. In this case, you just want to look at a list of files with the extension: .tif. Note that it is important that the file ends with .tif. So you use the dollar sign at the end of your pattern to tell R to only grab files that end with .tif.

pattern = ".tif$"

# but really you just want the tif files
all_landsat_bands <- list.files("data/week-07/landsat/LC80340322016205-SC20170127160728/crop",
                      pattern = ".tif$",
                      full.names = TRUE) # make sure you have the full path to the file
all_landsat_bands
##  [1] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_bqa_crop.tif"        
##  [2] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_cfmask_conf_crop.tif"
##  [3] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_cfmask_crop.tif"     
##  [4] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band1_crop.tif"   
##  [5] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band2_crop.tif"   
##  [6] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band3_crop.tif"   
##  [7] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band4_crop.tif"   
##  [8] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band5_crop.tif"   
##  [9] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band6_crop.tif"   
## [10] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band7_crop.tif"   
## [11] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_cloud_crop.tif"   
## [12] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_ipflag_crop.tif"

Above, you use the $ after .tif to tell R to look for files that end with .tif. This is a good start but there is one more condition that we’d like to meet. We only want the .tif files that are spectral bands. Notice that some of your files have text that includes “mask”, flags, etc. Those are all additional layers that you don’t need right now. You just need the spectral data saved in bands 1_7.

Mini Introduction to Regular Expressions

Thus, you want to grab all bands that both end with .tif AND contain the text “band” in them. To do this you use the function glob2rx() which allows you to specify both conditions using what is called regular expressions.

A regular expression is a sequence of characters that defines a search pattern. Here you tell R to select all files that have the word band in the filename. You use a * sign before and after band because you don’t know exactly what text will occur before or after band. You use .tif$ to tell R that each file needs to end with .tif.

all_landsat_bands <- list.files("data/week-07/landsat/LC80340322016205-SC20170127160728/crop",
           pattern = glob2rx("*band*.tif$"),
           full.names = TRUE) # use the dollar sign at the end to get all files that END WITH
all_landsat_bands
## [1] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band1_crop.tif"
## [2] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band2_crop.tif"
## [3] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band3_crop.tif"
## [4] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band4_crop.tif"
## [5] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band5_crop.tif"
## [6] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band6_crop.tif"
## [7] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band7_crop.tif"

Open the .tif Files in R

Now you have a list of all of the Landsat bands in your folder. You could chose to open each file individually using the raster() function.

# get first file
all_landsat_bands[2]
## [1] "data/week-07/landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band2_crop.tif"
landsat_band2 <- raster(all_landsat_bands[2])
plot(landsat_band2,
     main = "Landsat cropped band 2\nCold Springs fire scar",
     col = gray(0:100 / 100))

Landsat band 2 plot of Cold Springs fire scar.

However, that is not a very efficient approach. It’s more efficiently to open all of the layers together as a stack. Then you can access each of the bands and plot / use them as you want. You can do that using the stack() function.

# stack the data
landsat_stack_csf <- stack(all_landsat_bands)
# then turn it into a brick
landsat_csf_br <- brick(landsat_stack_csf)
# view stack attributes
landsat_csf_br
## class      : RasterBrick 
## dimensions : 177, 246, 43542, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 455655, 463035, 4423155, 4428465  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : LC8034032//band1_crop, LC8034032//band2_crop, LC8034032//band3_crop, LC8034032//band4_crop, LC8034032//band5_crop, LC8034032//band6_crop, LC8034032//band7_crop 
## min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values :                  3488,                  3843,                  4746,                  5152,                  5674,                  4346,                  3767

Let’s plot each individual band in your brick.

plot(landsat_csf_br,
     col = gray(20:100 / 100))

Plot all individual Landsat bands.

# get list of each layer name
names(landsat_csf_br)
## [1] "LC80340322016205LGN00_sr_band1_crop"
## [2] "LC80340322016205LGN00_sr_band2_crop"
## [3] "LC80340322016205LGN00_sr_band3_crop"
## [4] "LC80340322016205LGN00_sr_band4_crop"
## [5] "LC80340322016205LGN00_sr_band5_crop"
## [6] "LC80340322016205LGN00_sr_band6_crop"
## [7] "LC80340322016205LGN00_sr_band7_crop"

# remove the filename from each band name for pretty plotting
names(landsat_csf_br) <- gsub(pattern = "LC80340322016205LGN00_sr_", replacement = "", names(landsat_csf_br))
plot(landsat_csf_br,
     col = gray(20:100 / 100))

Plot all individual Landsat bands with clean names.

Plot RGB Composite Band Images with Landsat in R

Next, let’s 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.

par(col.axis = "white", col.lab = "white", tck = 0)
plotRGB(landsat_csf_br,
     r = 4, g = 3, b = 2,
     stretch = "lin",
     axes = TRUE,
     main = "RGB composite image\n Landsat Bands 4, 3, 2")
box(col = "white")

Plot Landsat RGB composite.

Now we’ve created a red, green blue color composite image. Remember this is what your eye would see. What happens if you plot the near infrared band instead of red? Try the following combination:

par(col.axis = "white", col.lab = "white", tck = 0)
plotRGB(landsat_csf_br,
     r = 5, g = 4, b = 3,
     stretch = "lin",
     axes = TRUE,
     main = "Color infrared composite image\n Landsat Bands 5, 4, 3\n Post Fire Data")
box(col = "white")

Plot Landsat CIR composite.

Optional Challenge

Using the ESRI Landsat 8 band combinations post as a guide. Plot the following Landsat band combinations:

  • False color
  • Color infrared
  • Agriculture
  • Healthy vegetation

Be sure to add a title to each of your plots that specifies the band combination.

Additional Resources

Faster Raster Calculations

Leave a Comment

You May Also Enjoy

Plot Data With Matplotlib

1 minute read

This Week

Calculate Seasonal Summary Values from Climate Data Variables Stored in NetCDF 4 Format: Work With MACA v2 Climate Data in Python

26 minute read

Learn how to calculate seasonal summary values for MACA 2 climate data using xarray and region mask in open source Python.

Calculate Summary Values Using Spatial Areas of Interest (AOIs) including Shapefiles for Climate Data Variables Stored in NetCDF 4 Format: Work With MACA v2 Climate Data in Python

32 minute read

Climate datasets stored in netcdf 4 format often cover the entire globe or an entire country. Learn how to subset climate data spatially and by time slices u...

How to Open and Process NetCDF 4 Data Format in Open Source Python

26 minute read

Historic and projected climate data are most often stored in netcdf 4 format. Learn how to open and process MACA version 2 climate data for the Continental U...