NEON EDUCATION bio photo

NEON EDUCATION

Devoted to open data and open source in science and education.

View All Tutorials

This tutorial is a part of a series!

Click below to view all lessons in the series!

Tags

R programming (52)
Hierarchical Data Formats (HDF5) (15)
Spatial Data & GIS (22)
LiDAR (10)
Raster Data (14)
Remote Sensing (25)
Data Visualization (4)
Hyperspectral Remote Sensing (18)
Time Series (15)
Phenology (7)
Vector Data (6)
Metadata (1)
Git & GitHub (7)
(1) (1) (14)

Tutorial by R Package

dplyr (7)
ggplot2 (16)
h5py (2)
lubridate (time series) (6)
maps (1)
maptools (1)
plyr (2)
raster (26)
rasterVis (raster time series) (3)
rgdal (GIS) (24)
rgeos (2)
rhdf5 (11)
sp (5)
scales (4)
gridExtra (4)
ggtheme (0)
grid (2)
reshape2 (3)
plotly (5)

View ALL Tutorial Series




Twitter Youtube Github


Blog.Roll

R Bloggers

Overview

R Skill Level: Intermediate

Goals / Objectives

After completing this activity, you will:
  1. Extract a "slice" of data from a hyperspectral data cube.
  2. Create a rasterstack in R which can then be used to create RGB images from bands in a hyperspectral data cube.
  3. Plot data spatially on a map.
  4. Create basic vegetation indices like NDVI using raster based calculations in R.

R Libraries to Install

  • rhdf5: source("http://bioconductor.org/biocLite.R") ; biocLite("rhdf5")
  • raster: install.packages('raster')
  • rgdal: install.packages('rgdal')
  • maps: install.packages('maps')

Data to Download

Download NEON Teaching Data Subset: Meteorological Data for Harvard Forest
The data used in this lesson were collected at the National Ecological Observatory Network's Harvard Forest field site. These data are proxy data for what will be available for 30 years on the NEON data portal for the Harvard Forest and other field sites located across the United States. </p>

Pre-reqs

We highly recommend you work through the - Introduction to Working with hyperspectral data in R activity before moving on to this activity.

About

We often want to generate a 3 band image from multi or hyperspectral data. The most commonly recognized band combination is RGB which stands for Red, Green and Blue. RGB images are just like the images that your camera takes. But there are other band combinations that are useful too. For example, near infrared images emphasize vegetation and help us classify or identify where vegetation is located on the ground.

SJER image using 3 different band combinations. Left: typical red, green and blue (bands 58,34,19), middle: color infrared: near infrared, green and blue (bands 90, 34, 19).

Data Tip - Band Combinations: the Biodiversity Informatics group created a great interactive tool that lets you explore band combinations. Check it out:Learn more about band combinations using a great online tool!

About This Activity

In this activity, we will learn how to create multi (3) band images. We will also learn how to perform some basic raster calculations (known as raster math in the GIS world).

Creating a Raster Stack in R

In the previous activity, we exported a subset of the NEON Reflectance data from a HDF5 file. In this activity, we will create a full color image using 3 (red, green and blue - RGB) bands. We will follow many of the steps we followed in the intro to working with hyperspectral data activity. These steps included loading required packages, reading in our file and viewing the file structure.

#Load required packages
  library(raster)
	library(rhdf5)


#Read in H5 file
f <- 'NEON-DS-Imaging-Spectrometer-Data.h5'
#View HDF5 file structure 
h5ls(f,all=T)

##   group        name         ltype corder_valid corder cset       otype
## 0     / Reflectance H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 1     /        fwhm H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 2     /    map info H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 3     / spatialInfo H5L_TYPE_HARD        FALSE      0    0   H5I_GROUP
## 4     /  wavelength H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
##   num_attrs  dclass          dtype  stype rank             dim
## 0         6 INTEGER  H5T_STD_I16LE SIMPLE    3 477 x 502 x 426
## 1         2   FLOAT H5T_IEEE_F32LE SIMPLE    2         426 x 1
## 2         1  STRING     HST_STRING SIMPLE    1               1
## 3        11                                  0                
## 4         2   FLOAT H5T_IEEE_F32LE SIMPLE    2         426 x 1
##            maxdim
## 0 477 x 502 x 426
## 1         426 x 1
## 2               1
## 3                
## 4         426 x 1

To spatially locate our raster data, we need a few key attributes:

  • The coordinate reference system
  • The lower left hand corner X, Y location of the raster
  • The dimensions (number of pixels in the x and y directions) of the raster to use to calculate the extent.

We’ll begin by grabbing these key attributes from the H5 file.

#r get spatial info and map info using the h5readAttributes function
spInfo <- h5readAttributes(f,"spatialInfo")
#define coordinate reference system
myCrs <- spInfo$projdef
#define the resolution
res <- spInfo$xscale

#Populate the raster image extent value. 
mapInfo<-h5read(f,"map info")
#the map info string contains the lower left hand coordinates of our raster
#let's grab those next
# split out the individual components of the mapinfo string
mapInfo<-unlist(strsplit(mapInfo, ","))

#grab the utm coordinates of the lower left corner
xMin<-as.numeric(mapInfo[4])
yMax<-as.numeric(mapInfo[5]) 

#r get attributes for the Reflectance dataset
reflInfo <- h5readAttributes(f,"Reflectance")

#create objects represents the dimensions of the Reflectance dataset
#note that there are several ways to access the size of the raster contained
#within the H5 file
nRows <- reflInfo$row_col_band[1]
nCols <- reflInfo$row_col_band[2]
nBands <- reflInfo$row_col_band[3]

#grab the no data value
myNoDataValue <- reflInfo$`data ignore value`
myNoDataValue

## [1] "15000"

Next, we’ll write a function that will perform the processing that we did step by step in the intro to working with hyperspectral data activity. This will allow us to process multiple bands in bulk.

The function band2Rast slices a band of data from the HDF5 file, and extracts the reflectance. It them converts the data to a matrix, converts it to a raster and returns a spatially corrected raster for the specified band.

The function requires the following variables:

  • file: the file
  • band: the band number we wish to extract
  • noDataValue: the noDataValue for the raster
  • xMin, yMin: the X,Y coordinate left hand corner locations for the raster.
  • res: the resolution of the raster
  • crs: the Coordinate Reference System for the raster

The function output is a spatially referenced, r raster object.

# file: the hdf file
# band: the band you want to process
# returns: a matrix containing the reflectance data for the specific band

band2Raster <- function(file, band, noDataValue, xMin, yMin, res, crs){
    #first read in the raster
    out<- h5read(f,"Reflectance",index=list(1:nCols,1:nRows,band))
	  #Convert from array to matrix
	  out <- (out[,,1])
	  #transpose data to fix flipped row and column order 
    #depending upon how your data are formated you might not have to perform this
    #step.
	  out <-t(out)
    #assign data ignore values to NA
    #note, you might chose to assign values of 15000 to NA
    out[out == myNoDataValue] <- NA
	  
    #turn the out object into a raster
    outr <- raster(out,crs=myCrs)
 
    # define the extents for the raster
    #note that you need to multiple the size of the raster by the resolution 
    #(the size of each pixel) in order for this to work properly
    xMax <- xMin + (outr@ncols * res)
    yMin <- yMax - (outr@nrows * res)
 
    #create extents class
    rasExt  <- extent(xMin,xMax,yMin,yMax)
   
    #assign the extents to the raster
    extent(outr) <- rasExt
   
    #return the raster object
    return(outr)
}

Now that the function is created, we can create our list of rasters. The list specifies which bands (or dimensions in our hyperspectral dataset) we want to include in our raster stack. Let’s start with a typical RGB (red, green, blue) combination. We will use bands 58, 34, and 19.

Data Tip - wavelengths and bands: Remember that you can look at the wavelengths dataset to determine the center wavelength value for each band.

#create a list of the bands we want in our stack
rgb <- list(58,34,19)
#lapply tells R to apply the function to each element in the list
rgb_rast <- lapply(rgb,band2Raster, file = f, 
                   noDataValue=myNoDataValue, 
                   xMin=xMin, yMin=yMin, res=1,
                   crs=myCrs)

#check out the properties or rgb_rast
#note that it displays properties of 3 rasters.

rgb_rast

## [[1]]
## class       : RasterLayer 
## dimensions  : 502, 477, 239454  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 256521, 256998, 4112069, 4112571  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11N +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 123, 15453  (min, max)
## 
## 
## [[2]]
## class       : RasterLayer 
## dimensions  : 502, 477, 239454  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 256521, 256998, 4112069, 4112571  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11N +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 116, 15677  (min, max)
## 
## 
## [[3]]
## class       : RasterLayer 
## dimensions  : 502, 477, 239454  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 256521, 256998, 4112069, 4112571  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11N +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 84, 13805  (min, max)

#finally, create a raster stack from our list of rasters
hsiStack <- stack(rgb_rast)

More about Lapply here.

NOTE: We are using the raster stack object in R to store several rasters that are of the same CRS and extent.

Next, add the names of the bands to the raster so we can easily keep track of the bands in the list.

#Add the band numbers as names to each raster in the raster list

#Create a list of band names
bandNames <- paste("Band_",unlist(rgb),sep="")

names(hsiStack) <- bandNames
#check properties of the raster list - note the band names
hsiStack

## class       : RasterStack 
## dimensions  : 502, 477, 239454, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 256521, 256998, 4112069, 4112571  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11N +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## names       : Band_58, Band_34, Band_19 
## min values  :     123,     116,      84 
## max values  :   15453,   15677,   13805

#scale the data as specified in the reflInfo$Scale Factor

hsiStack <- hsiStack/reflInfo$`Scale Factor`

##  Plot one raster in the stack to make sure things look OK.
plot(hsiStack$Band_58, main="Band 58")

We can play with the color ramps too if we want:

#change the colors of our raster 
myCol=terrain.colors(25)
image(hsiStack$Band_58, main="Band 58", col=myCol)

#adjust the zlims or the stretch of the image
myCol=terrain.colors(25)
image(hsiStack$Band_58, main="Band 58", col=myCol, zlim = c(0,.5))

#try a different color palette
myCol=topo.colors(15, alpha = 1)
image(hsiStack$Band_58, main="Band 58", col=myCol, zlim=c(0,.5))

The plotRGB function allows you to combine three bands to create an image. More on plotRGB here.

# create a 3 band RGB image
plotRGB(hsiStack,
        r=1,g=2,b=3, scale=300, 
        stretch = "Lin")

A note about image stretching: Notice that the scale is set to 300 on the RGB image that we plotted above. We can adjust this number and notice that the image gets darker - or lighter.

Once you’ve created your raster, you can export it as a geotiff.

#write out final raster	
#note - you can bring this tiff into any GIS program!
#note: if you set overwrite to TRUE, then you will overwite or lose the older
#version of the tif file! keep this in mind.
writeRaster(hsiStack, file="rgbImage.tif", format="GTiff", overwrite=TRUE)

Data Tip - False color and near infrared images: Use the band combinations listed at the top of this page to modify the raster list. What type of image do you get when you change the band values?

Challenge

Use different band combinations to create other “RGB” images. Suggested band combinations are below:

  • Color Infrared / False Color: rgb: (90,34,19)
  • SWIR, NIR,Red Band – rgb (152,90,58)
  • False Color: – rgb (363,246,55)

More on Band Combinations: http://gdsc.nlr.nl/gdsc/en/information/earth_observation/band_combinations

Plotting spectral data on a map.

We can plot the location of our image on a map of the US. For this we’ll use the lower left coordinates of the raster, extracted from the SPINFO group. Note that these coordinates are in latitude and longitude (geographic coordinates) rather than UTM coordinates.

#Create a Map showing the location of our dataset in R
library(maps)
map(database="state",region="california")
points(spInfo$LL_lat~spInfo$LL_lon,pch = 15)
#add title to map.
title(main="NEON San Joaquin Field Site - Southern California")

Raster Math - Creating NDVI and other Vegetation Indices in R

In this last part, we will calculate some vegetation indices using raster math in R! We will start by creating NDVI or Normalized Difference Vegetation Index.

Data Tip - About NDVI: NDVI is a ratio between the near infrared (NIR) portion of the electromagnetic spectrum and the red portion of the spectrum. Please keep in mind that there are different ways to aggregate bands when using hyperspectral data. This example is using individual bands to perform the NDVI calculation. Using individual bands is not necessarily the best way to calculate NDVI from hyperspectral data!

#Calculate NDVI
#select bands to use in calculation (red, NIR)
ndvi_bands <- c(58,90)


#create raster list and then a stack using those two bands
ndvi_rast <- lapply(ndvi_bands,band2Raster, file = f, noDataValue=15000, 
                    xMin=xMin, yMin=yMin,
                    crs=myCRS,res=1)
ndvi_stack <- stack(ndvi_rast)

#make the names pretty
bandNDVINames <- paste("Band_",unlist(ndvi_bands),sep="")
names(ndvi_stack) <- bandNDVINames

ndvi_stack

## class       : RasterStack 
## dimensions  : 502, 477, 239454, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 256521, 256998, 4112069, 4112571  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11N +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## names       : Band_58, Band_90 
## min values  :     123,     315 
## max values  :   15453,   15293

#calculate NDVI
NDVI <- function(x) {
	  (x[,2]-x[,1])/(x[,2]+x[,1])
}
ndvi_calc <- calc(ndvi_stack,NDVI)
plot(ndvi_calc, main="NDVI for the NEON SJER Field Site")

#play with breaks and colors to create a meaningful map

#add a color map with 5 colors
myCol <- terrain.colors(3)
#add breaks to the colormap (6 breaks = 5 segments)
brk <- c(0, .4, .7, .9)

#plot the image using breaks
plot(ndvi_calc, main="NDVI for the NEON SJER Field Site", col=myCol, breaks=brk)

LEFT: NDVI for the NEON SJER field site, created in R. RIGHT: EVI for the NEON SJER field site created in R.

Challenge

Try any of the following:

  1. Calculate EVI using the following formula : EVI<- 2.5 * ((b4-b3) / (b4 + 6 * b3- 7.5*b1 + 1))
  2. Calculate NDNI using the following equation: log(1/p1510)-log(1/p1680)/ log(1/p1510)+log(1/p1680)
  3. Explore the bands in the hyperspectral data. What happens if you average reflectance values across multiple red and NIR bands and then calculate NDVI?

Get Lesson Code

(some browsers may require you to right click.)