Visualizing Landsat Data With R
EarthExplorer
The US Geological Service's EarthExplorer is a portal for downloading a variety of remotely sensed data.
You will need to create a (free) user account and password in order to download data.
After selecting a point on the map, you will be given an option of data sets to choose from.
For recent data, you should choose the Landsat Archive, Collection 1 Level-1 and then Landsat 8 (L8 OLI/TIRS). Collections utilize a naming and quality ranking scheme that was implemented in 2017. Data organized using the older system is referred to as pre-collection.
Because Landsat only passes over any given location on the surface of the earth every 16 days, if that happens to be a cloudy day at your point of interest on the day of interest, you will be out of luck. Therefore finding a clear scene will often be a challenge, especially in cloudy northern climates.
Once you have found a scene that looks promising, you will be given a list of different types of downloads:
- LandsatLook Images with Geographic Reference should be a choice if you need color visible light imagery
- Level-1 GeoTIFF Data Product is the appropriate product if you are using other bands, such as the near infrared band needed to calculate NDVI
Memory and Bandwidth Considerations
Note that rasters are often large files that take a long time to load and consume considerable disk space, memory, and connection bandwith.
If you are storing your data on a network drive, you should consider downloading and performing your work on your local hard drive, and then saving your work at completion to the network drive. Having to constantly transfer your data across a network can consume a considerable amount of unnecessary time.
You will also likely find performance better if connected to the internet through a wired connection rather than Wi-Fi.
Rasters in R
There are three diffent Raster* structures that you choose from depending on how many bands you need for a scene, and how many files they come from:
- RasterLayer = one file, one band
- RasterBrick = one file, multiple bands. Commonly used with RGB images
- RasterStack = multiple files, multiple bands (Landsat GeoTIFFs)
For example, the LandsatLook RGB images are distributed as a single GeoTIFF file (*T1.tif) containing three bands (red, green, blue). Use the brick() function to load the image into memory and the plotRGB function to display it as a single image. The exact name of the file will depend on what location and date you download.
library(raster)
image = brick("LC08_L1TP_043027_20150330_20170228_01_T1.tif")
plotRGB(image)
In contrast, a thermal infrared (surface temperature) raster (with a TIR at the end of the name) included in the .zip archive with the LandsatLook image is a single-band file that can be loaded with the raster() function and displayed with the regular plot() (not plotRGB()) function.
The Kelvin temperature values calculated from the sensor are fitted to a value range of 0-255 for smaller storage, so the values in this raster have only relative value for plotting and have no absolute value significance.
thermal = raster("LC08_L1TP_043027_20150330_20170228_01_T1_TIR.tif")
plot(thermal)
Single band images are commonly viewed in false-color, which assigns different colors to values (as with a choropleth) that makes areas of clustered high or low values more-easily visible.
Note that, unlike the plot() function with points, lines or polygons where the col= parameter is used to pass individual colors for individual features, by default the col= parameter when plotting a raster is treated as a set of colors (ranked from low to high) where plot() assigns specific ranges of values for to each palette color. This default behavior can be altered with the breaks parameter as described later.
thermal = raster('LC08_L1TP_043027_20150330_20170228_01_T1_TIR.tif')
palette = colorRampPalette(c('blue', 'white', 'red'))(256)
plot(thermal, col=palette)
A LandsatLook archive also contains a quality file (QB.tif at the end of the name) with assessment of the quality of each pixel. The range of values is from zero to 255, with higher values usually indicating cloud cover that obscures the surface of the earth.
# Quality image
quality = raster("LC08_L1TP_043027_20150330_20170228_01_T1_QB.tif")
palette = colorRampPalette(c("blue", "white", "red"))(256)
plot(quality, col=palette)
Multispectral Data and NDVI
Scenes from the Level-1 GeoTIFF Data Product are distributed as .tar.gz files, which are Unix gzipped tape archives that can be opened with a program like 7-Zip. Each archive contains individual GeoTIFF files for each of the 11 bands captured by the Landsat sensor.
It is possible to construct an RGB RasterStack using the red, green and blue bands (numbers 4, 3 and 2 respectively), although the color-enhanced LandsatLook rasters are generally preferred if you need visible-light color imagery.
blue = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B2.TIF")
green = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B3.TIF")
red = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B4.TIF")
rgb = stack(red, green, blue)
plotRGB(rgb, stretch='hist')
The Landsat multispectral data is of considerable value for the remote sensing of the environment.
One commonly used technique is the calculation of Normalized Difference Vegetation Index (NDVI), which can be used to detect the presence of and changes in photosynthetic vegetation.
NDVI = (infrared - red) / (infrared + red)
NDVI is based on the observation that green plants reflect infrared light (to avoid overheating), but absorb visible light (to power photosynthesis). This gives green plants a light reflection signature that can distinguish them from objects that absorb all spectra of light or different spectra of light from plants.
Rasters with the same extents and pixel dimensions can be used in mathematical operations like vectors or matrices. This is called raster algebra.
# Load visible (red) and infrared bands
red = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B4.TIF")
infrared = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B5.TIF")
# Perform raster algebra to calculate a raster of NDVI values
ndvi = (infrared - red) / (infrared + red)
# Plot the NDVI as a false-color image
palette = colorRampPalette(c("blue", "white", "red"))(256)
plot(ndvi, col=palette)
The range of NDVI is -1 to +1, although NDVI values will generally be in a small area in the middle of the full potential range of -1 to +1.
hist(ndvi, maxpixels = ncell(ndvi))
In such cases, visualization may be improved by passing a breaks parameter to plot() that assigns colors to specific ranges for plotting:
breaks = seq(-0.4, 0.4, 0.1)
palette = colorRampPalette(c("blue", "white", "red"))(8)
plot(ndvi, breaks=breaks, col=palette)
Plotting Over a Base Map For Context
Since the ultimate end of such analysis is usually to know what is going on at a specific location or set of locations, it can be helpful to plot NDVI over a base map that gives geographic context to the data:
library(OpenStreetMap)
# Find extent of raster in lat/long
wgs84 = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
outline = extent(projectExtent(ndvi, wgs84))
upperLeft = c(extent(outline)@ymax, extent(outline)@xmin)
lowerRight = c(extent(outline)@ymin, extent(outline)@xmax)
# Plot a monochrome base map
basemap = openmap(upperLeft, lowerRight, type='stamen-toner')
plot(basemap)
# Calculate NDVI
red = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B4.TIF")
infrared = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B5.TIF")
ndvi = (infrared - red) / (infrared + red)
# Reproject the NDVI raster to web mercator projection
ndvi = projectRaster(ndvi, projectExtent(ndvi, osm()))
# Plot the NDVI with a semi-transparent (alpha) palette
breaks = seq(-0.5, 0.5, 0.1)
palette = colorRampPalette(c("#0000FFC0", "#FFFFFFC0", "#FF0000C0"), alpha=T)(10)
plot(ndvi, breaks = breaks, col=palette, add=T)
Cropping To Smaller Areas
Landsat scenes cover a fairly large area of around 10,000 square miles, with each pixel representing about 10,000 square feet.
If you are interested in a smaller, more-specific area, the raster crop() function can be used to isolate a portion of the raster.
In the example below, the NDVI raster is cropped to the urban area around downtown Spokane. The low NDVI of water and paved urban areas clearly contrasts against the high NDVI of parks and greenspaces.
# Create a base map of Spokane and plot it
upperLeft = c(47.676, -117.453)
lowerRight = c(47.646, -117.381)
basemap = openmap(upperLeft, lowerRight, type='stamen-toner')
plot(basemap)
# Reproject the NDVI raster to web mercator projection
ndvi = projectRaster(ndvi, projectExtent(ndvi, osm()))
# Crop the raster to the region displayed on the base map
spokane = crop(ndvi, raster(basemap))
# Plot the cropped raster
breaks = seq(-0.5, 0.5, 0.1)
palette = colorRampPalette(c("#0000FFC0", "#FFFFFFC0", "#FF0000C0"), alpha=T)(10)
plot(spokane, breaks=breaks, col=palette, add=T)
Plotting NDVI Changes
The same type of raster algebra used to calculate NDVI can be used to calculate changes in NDVI between scenes.
When looking for a scene for for comparison a different date, it is important that you use an image from the same path and row. Although Landsat scenes overlap, using the same path and row will maximize the amount of common area and, therefore the amount of area covered in the analysis.
The raster showing the difference between the summer and winter scenes shows the clearly higher NDVI in the later summer scene.
# Create a base map of Spokane
upperLeft = c(47.676, -117.453)
lowerRight = c(47.646, -117.381)
basemap = openmap(upperLeft, lowerRight, type='stamen-toner')
plot(basemap)
# Calculate NDVI for the winter scene
red = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B4.TIF")
infrared = raster("LC08_L1TP_043027_20150330_20170228_01_T1_B5.TIF")
winter = (infrared - red) / (infrared + red)
# Reproject, crop and plot
winter = projectRaster(winter, projectExtent(winter, osm()))
winter = crop(winter, raster(basemap))
# Color palette for plotting
breaks = seq(-0.5, 0.5, 0.1)
palette = colorRampPalette(c("#0000FFC0", "#FFFFFFC0", "#FF0000C0"), alpha=T)(10)
plot(winter, breaks=breaks, col=palette, add=T)
# Calculate NDVI for summer scene
red = raster("LC08_L1TP_043027_20160823_20170223_01_T1_B4.TIF")
infrared = raster("LC08_L1TP_043027_20160823_20170223_01_T1_B5.TIF")
summer = (infrared - red) / (infrared + red)
# Reproject, crop and plot
summer = projectRaster(summer, projectExtent(summer, osm()))
summer = crop(summer, raster(basemap))
dev.new()
plot(basemap)
plot(summer, breaks=breaks, col=palette, add=T)
# Calculate the difference and plot it
change = summer - winter
dev.new()
breaks = seq(-0.2, 0.2, 0.05)
palette = colorRampPalette(c("#0000FFC0", "#FFFFFFC0", "#FF0000C0"), alpha=T)(8)
plot(basemap)
plot(change, breaks = breaks, col=palette, add=T)
# Put an annotation on an area of significant change
coords = data.frame(x = -117.419378, y = 47.662274)
point = SpatialPoints(coords, proj4string=wgs84)
point = spTransform(point, osm())
points(point, pch = 19, col='black')
text(point, pos=4, font=2, labels="Riverfront Park")
