R Script & Challenge Code: NEON data lessons often contain challenges that reinforce
learned skills. If available, the code for challenge solutions is found in the
downloadable R script of the entire lesson, available in the footer of each lesson page.

Recommended Skills

About Hyperspectral Data

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.

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.

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. You can bring
this GeoTIFF into any GIS program!

#write out final raster
#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: Other band combinations

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)

Plot 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
map(database="state", region="california")
points(spInfo$LL_lat~spInfo$LL_lon,pch = 15)
#add title to map.
title(main="NEON San Joaquin (SJER) Field Site - Central 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.

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!