Landsat Surface Reflectance CDR preprocessing with teamlucc

Alex Zvoleff bio photo By Alex Zvoleff


This post outlines how to use the teamlucc package to preprocess imagery from the Landsat Surface Reflectance Climate Data Record (CDR) archive. The teamlucc package supports a number of common preprocessing steps, including:

  • Conversion of CDR files to any GDAL supported file format
  • Cropping Landsat tiles to a given area of interest (AOI)
  • Mosaicking and cropping of DEM tiles (such as ASTER or SRTM) to a given AOI or Landsat path/row
  • Topographic correction of CDR scenes

Getting started

First load the teamlucc package, and the rgdal package:

## Loading required package: Rcpp
## Loading required package: raster
## Loading required package: sp
## Warning: replacing previous import by 'raster::buffer' when loading
## 'teamlucc'
## Warning: replacing previous import by 'raster::interpolate' when loading
## 'teamlucc'
## Warning: replacing previous import by 'raster::rotated' when loading
## 'teamlucc'
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
## Path to GDAL shared files: C:/Users/azvoleff/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/azvoleff/R/win-library/3.1/rgdal/proj

DEM setup

Before CDR surface reflectance images can be topographically corrected with teamlucc, a digital elevation model (DEM) needs to be obtained that is in the same resolution and coordinate system as the CDR imagery, and slope and aspect need to be calculated. The auto_setup_dem function in the teamlucc package facilitates this task. There are many options you can supply to auto_setup_dem - see ?auto_setup_dem for more information on these options. If you do not want to topographically correct your imagery, you can skip this step and move ahead to the next section.

auto_setup_dem assembles a DEM to cover a given area of interest (AOI), and can mosaic multiple DEM files if needed to cover an AOI. To use auto_setup_dem, the user must first define the area of interest with an AOI polygon. If you have a shapefile of an area of interest, load it into R using the readOGR command. The readOGR command needs the folder the shapefile is in (in this example the current working directory, specified by .) as the first parameter, and the filename of the shapefile (without the “.shp” extension) as the second parameter (PA_VB) in this example). This example uses a shapefile of the boundary of the Braulio Carrillo National Park, in which the Volcán Barva TEAM site is located:

aoi <- readOGR('.', 'PA_VB')
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "PA_VB"
## with 5 features and 8 fields
## Feature type: wkbPolygon with 2 dimensions


As seen in the above plot, there are a number of adjoining polygons in this shapefile. The functions in teamlucc expect the AOI to be of length one. So calculate the convex hull of the AOI using the functions in the rgeos package:

## rgeos version: 0.3-8, (SVN revision 460)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
aoi <- gConvexHull(aoi)

The AOI is now a single polygon:



If you do not have an AOI, but know the Landsat path and row you want to work with, an alternative is to define the AOI based on the path and row, using the wrspathrow R package, and supplying the desired path and row numbers. Here 127 is the WRS-2 path number, and 47 is the WRS-2 row number for the path/row at the center of the above AOI:

aoi_from_pathrow <- pathrow_poly(127, 47)

In addition to the AOI, auto_setup_dem needs to know the location and spatial extents of the DEM files you have available on your machine. This list can be assembled automatically using the get_extent_polys function in teamlucc.
See ?get_extent_polys for more information. The below code will fail on your machine because you will not have the proper DEMs for this example. Download the proper DEMs for the area in which you are working and place them in a folder on your machine if you want to test this function.

dem_files <- dir('H:/Data/CGIAR_SRTM/Tiles', pattern='.tif$', full.names=TRUE)
dems <- lapply(dem_files, raster)
dem_extents <- get_extent_polys(dems)

For flexibility, auto_setup_dem can optionally crop each output DEM to the spatial extent of the supplied AOI. If crop_to_aoi=TRUE, then auto_setup_dem will crop the DEMs to the spatial extent of the supplied AOI. If crop_to_aoi=FALSE, then auto_setup_dem will crop the DEMs to the extent of the Landsat path/rows needed to cover the AOI.

Lastly, auto_setup_dem needs to know where to save its output. For this example, save the output to the current working directory .. Now that all of the essential inputs are defined, auto_setup_dem can be called:

auto_setup_dem(aoi, '.', dem_extents, crop_to_aoi=TRUE)
## [1] "2014-11-13 16:49:49: started \"Processing DEMS for 1 path/rows\""
## [1] "2014-11-13 16:49:50: started \"Processing 1 of 1: 015-053\""
## [1] "2014-11-13 16:50:00: finished \"Processing 1 of 1: 015-053\" (10.047s elapsed)"
## [1] "2014-11-13 16:50:00: finished \"Processing DEMS for 1 path/rows\" (10.285s elapsed)"

The result will be a mosaicked DEM, in the current working directory, that can be used for topographically correcting imagery with the auto_preprocess_landsat function.


Images are delivered from the CDR archive in either ENVI format, GeoTIFF format, or Hierarchical Data Format (HDF). The teamlucc package will, by default, convert HDF or ENVI format images to a GeoTIFF format, as these image files can be easily read in R and in most commonly used remote sensing and GIS software packages. This example assumes you want GeoTIFF output - see the help files for teamlucc for other output options.

First you will need to acquire a time series of CDR imagery for your site. The espa_download function can facilitate downloading files from an ESPA order. Also see the post on Filtering and downloading Landsat scenes for more details on how teamlucc can help with image acquisition.

Once you have downloaded your imagery from ESPA, I recommend you put all of the zip files from your download in a single folder. You can then use the teamlucc espa_extract function to automate extracting your CDR image files, including placing the extracted files in consistently named output folders.
If, for example, your files are located in espa_downloads and you want to extract them to espa_extracts, run the following (this block, and the next block of code, will both fail on your computer since you do not have the required imagery - this is only an example, download your own imagery to follow along):

download_folder <- 'espa_downloads'
extract_folder <- 'espa_extracts'
espa_extract(download_folder, extract_folder)
## 1 of 2. Extracting LT50150532000044-SC20140514195145.tar.gz to espa_extracts/015-053_2000-044_LT5
## 2 of 2. Extracting LT50150532001014-SC20140514195632.tar.gz to espa_extracts/015-053_2001-014_LT5

Now that the CDR format image files are extracted, you are ready to run auto_preprocess_landsat. As with auto_setup_dem, there are many options you can supply to auto_preprocess_landsat - see ?auto_preprocess_landsat. The below is a simple example of how to call auto_preprocess_landsat.

The image_dirs line below is just a fancy way of finding all the folders located in extract_folder that contain CDR imagery. You could just as easily specify these folders individually as as a list of strings, like: image_dirs <- c('C:/folder1', 'C:/folder2') if you had two CDR Landsat scenes located in C:/folder1 and C:/folder2, respectively.

The prefix parameter specifies a string that will be used in naming files output by auto_preprocess_landsat. For prefix I suggest you use a short (2 or 3 character) site name or site code that is meaningful to you.

There are two other options we provide below to auto_preprocess_landsat.
tc=TRUE tells auto_preprocess_landsat to perform topographic correction.
Because of this, we also need to specify dem_path (where the DEM files preprocessed by auto_setup_dem are located), so that the DEM files for this scene can be found. Here we set dem_path='.' as the DEM is in our current working directory. We supply an AOI (same AOI as above) to crop the output images. verbose=TRUE indicates that we want detailed progress messages to print while the script is running. The output of auto_preprocess_landsat is a data.frame with a list of the preprocessed files and their file formats.

image_dirs <- dir('espa_extracts',
filelist <- auto_preprocess_landsat(image_dirs, prefix='VB', tc=TRUE, 
                                    dem_path='.', aoi=aoi, verbose=TRUE)
## Warning: executing %dopar% sequentially: no parallel backend registered
## [1] "2014-11-13 16:50:52: started \"Preprocessing 015-053_2000-044_L5TSR\""
## [1] "2014-11-13 16:50:52: started \"cropping and reprojecting\""
## Warning in build_mask_vrt(file_base, mask_vrt_file, file_format): Using
## "fmask_band" instead of newer "cfmask_band" band name
## [1] "2014-11-13 16:51:03: finished \"cropping and reprojecting\" (10.787s elapsed)"
## [1] "2014-11-13 16:51:03: started \"topocorr\""
## [1] "2014-11-13 16:53:25: finished \"topocorr\" (142.071s (~2.37 minutes) elapsed)"
## [1] "2014-11-13 16:53:25: started \"writing data\""
## [1] "2014-11-13 16:53:29: finished \"writing data\" (4.052s elapsed)"
## [1] "2014-11-13 16:53:29: finished \"Preprocessing 015-053_2000-044_L5TSR\" (156.919s (~2.62 minutes) elapsed)"
## [1] "2014-11-13 16:53:30: started \"Preprocessing 015-053_2001-014_L5TSR\""
## [1] "2014-11-13 16:53:30: started \"cropping and reprojecting\""
## Warning in build_mask_vrt(file_base, mask_vrt_file, file_format): Using
## "fmask_band" instead of newer "cfmask_band" band name
## [1] "2014-11-13 16:53:44: finished \"cropping and reprojecting\" (14.161s elapsed)"
## [1] "2014-11-13 16:53:44: started \"topocorr\""
## [1] "2014-11-13 16:55:55: finished \"topocorr\" (131.256s (~2.19 minutes) elapsed)"
## [1] "2014-11-13 16:55:55: started \"writing data\""
## [1] "2014-11-13 16:55:59: finished \"writing data\" (4.236s elapsed)"
## [1] "2014-11-13 16:55:59: finished \"Preprocessing 015-053_2001-014_L5TSR\" (149.662s (~2.49 minutes) elapsed)"

The result is two cropped, reprojected, and topographically corrected Landsat images covering the specified AOI. One image from 2000:

ls_2000 <-  brick('espa_extracts/015-053_2000-044_LT5/VB_015-053_2000-044_L5TSR_tc.tif')
ls_2000 <- linear_stretch(ls_2000, pct=2, max_val=255)
browse_image(ls_2000, r=4, g=3, b=2)


And one image from 2001:

ls_2001 <- brick('espa_extracts/015-053_2001-014_LT5/VB_015-053_2001-014_L5TSR_tc.tif')
ls_2001 <- linear_stretch(ls_2001, pct=2, max_val=255)
browse_image(ls_2001, r=4, g=3, b=2)


There is a fair amount of cloud cover in the 2000 image. See the post on cloud removal for one means of addressing this issue.