Lesson 3. Landsat remote sensing tif files in R

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 6 of the course.

Download Week 6/7 Data (~500 MB)

In the previous lesson, we learned how to import a multi-band image into R using the stack() function. We 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 we will learn how to work with Landsat data in R. In this case, our 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: <a href=”USGS

Landsat 40 year timeline source: USGS.
The 40 year history of landsat missions. Source: USGS - <a href="https://landsat.usgs.gov/landsat-missions-timeline

Landsat data is a spectral dataset, collected from space. The spectral 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

Understanding landsat

When working with landsat, it is important to understand both the metadata and the file naming convention. The metadata tell us about how the data were processed, where the data are from and how they are structured.

The file names, tell us what sensor collected the data, the date the data were collected, and more.

<a href=”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 underanding 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

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 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 <a href=”convert julian days to actual date.

Landsat tif files in R

Now that we understand how our file is named.

# load spatial packages
library(raster)
library(rgdal)
library(rgeos)
# turn off factors
options(stringsAsFactors = F)

If we look at the directory that contains our landsat data, we will see that each of the individual bands is stored individually as a geotiff rather than being stored as a stacked or layered 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.

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

Get list of files

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

# get list of all tifs
list.files("data/week_06/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"

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

pattern=".tif$"

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

Above, we 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 our files have text that includes “mask”, flags, etc. Those are all additional layers that we don’t need right now. We just need the spectral data saved in bands 1_7.

Thus, we want to grab all bands that both end with .tif AND contain the text “band” in them. To do this we use the function glob2rx() which allows us to specify both conditions. Here we tell R to select all files that have the word band in the filename. We use a * sign before and after band because we don’t know exactly what text will occur before or after band. We use .tif$ to tell R that each file needs to end with .tif.

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

Now we have a list of all of the landsat bands in our folder. We could chose to open each file individually using the raster() function.

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

Landsat band 2 plot

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

# stack the data
landsat_stack_csf <- stack(all_landsat_bands)
# view stack attributes
landsat_stack_csf
## class       : RasterStack 
## dimensions  : 177, 246, 43542, 7  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 455655, 463035, 4423155, 4428465  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## 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 our stack.

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

plot individual landsat bands

# get list of each layer name
names(landsat_stack_csf)
## [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_stack_csf) <- gsub(pattern = "LC80340322016205LGN00_sr_", replacement = "", names(landsat_stack_csf))
plot(landsat_stack_csf,
     col=gray(20:100 / 100))

plot individual landsat bands good names

Plot RGB image

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 <a href=”ESRI landsat 8 band combinations post.

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

plot rgb composite

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

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

plot rgb composite

optional challenge

Using the <a href=”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

Leave a Comment