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

LiDAR (9)
R programming (70)
Remote Sensing (12)
Data Visualization (4)
Hyperspectral Remote Sensing (7)
Hierarchical Data Formats (HDF5) (24)
Spatial Data & GIS (18)
Time Series (15)
Phenology (7)
Raster Data (8)
Vector Data (6)
Metadata (1)
Git & GitHub (6)
(1) (1) (1)

Tutorial by R Package

dplyr (8)
ggplot2 (17)
h5py (1)
lubridate (time series) (6)
maps (1)
maptools (3)
plyr (2)
raster (32)
rasterVis (raster time series) (3)
rgdal (GIS) (23)
rgeos (5)
rhdf5 (21)
sp (7)
scales (4)
gridExtra (4)
ggtheme (0)
grid (2)
reshape2 (3)
plotly (6)

View ALL Tutorial Series




Twitter Youtube Github


Blog.Roll

R Bloggers
R Skill Level: Intermediate

Goals / Objectives

After completing this activity, you will be able to:
  1. Subset an HDF5 file.

Getting Started

In this tutorial, we will extract a single-pixel’s worth of reflectance values from an HDF5 file and plot a spectral profile for that pixel.

# load packages
library(rhdf5)
library(raster)
library(plyr)
library(rgeos)
library(rgdal)
library(ggplot2)

# be sure to set your working directory
# setwd("~/Documents/data/NEONDI-2016") # Mac
# setwd("~/data/NEONDI-2016")  # Windows

Import H5 Functions

We have built a suite of functions that allow us to quickly open and read NEON hyperspectral imagery data from an hdf5 file.

# install devtools (only if you have not previously intalled it)
# install.packages("devtools")
# call devtools library
#library(devtools)

## install from github
# install_github("lwasser/neon-aop-package/neonAOP")
## call library
library(neonAOP)



# your file will be in your working directory! This one happens to be in a diff dir
# than our data
# source("/Users/lwasser/Documents/GitHub/neon-aop-package/neonAOP/R/aop-data.R")

Open H5 File

Next, let’s define a few variables that we will need to access the H5 file.

# Define the file name to be opened
f <- "NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# Look at the HDF5 file structure
h5ls(f, all=T)

# define the CRS in EPGS format for the file
epsg <- 32611

Read Wavelength Values

Next, let’s read in the wavelength center associated with each band in the HDF5 file.

# read in the wavelength information from the HDF5 file
wavelengths<- h5read(f,"wavelength")
# convert wavelength to nanometers (nm)
# NOTE: this is optional!
wavelengths <- wavelengths*1000

Get Subset

Next, we might want to extract just a subset of our data to pull a spectral signature. For example, a plot boundary.

# get of H5 file map tie point
h5.ext <- create_extent(f)

# turn the H5 extent into a polygon to check overlap
h5.ext.poly <- as(extent(h5.ext), "SpatialPolygons")

# open file clipping extent
clip.extent <- readOGR("NEONdata/D17-California/TEAK/vector_data", "TEAK_plot")

## OGR data source with driver: ESRI Shapefile 
## Source: "NEONdata/D17-California/TEAK/vector_data", layer: "TEAK_plot"
## with 1 features
## It has 1 fields

# assign crs to h5 extent
# paste0("+init=epsg:", epsg) -- it is better to use the proj string here

crs(h5.ext.poly) <- CRS("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")


# ensure the two extents overlap
gIntersects(h5.ext.poly, clip.extent)

## [1] TRUE

# if they overlap, then calculate the extent
# this doesn't currently account for pixel size at all 
# and these values need to be ROUNDED

yscale <- 1
xscale <- 1

# define index extent
# xmin.index, xmax.index, ymin.index,ymax.index
# all units will be rounded which means the pixel must occupy a majority (.5 or greater)
# within the clipping extent
index.bounds <- calculate_index_extent(extent(clip.extent),
																			 h5.ext)

# open a band that is subsetted using the clipping extent
b58_clipped <- neonAOP::open_band(fileName=f,
								bandNum=58,
								epsg=32611,
								subsetData = TRUE,
								dims=index.bounds)

# plot clipped bands
plot(b58_clipped,
     main="Band 58 Clipped")

Run Subset over Many Bands

# within the clipping extent
index.bounds <- calculate_index_extent(extent(clip.extent),
																			 h5.ext)
# create  alist of the bands
bands <- list(19,34,58)
# within the clipping extent
index.bounds <- calculate_index_extent(extent(clip.extent),
              h5.ext)
# clip out raster
rgbRast.clip <- neonAOP::create_stack(file=f,
            bands=bands,
            epsg=epsg,
            subset=TRUE,
            dims=index.bounds)

plotRGB(rgbRast.clip,
        stretch="lin")

Unclipped Data

And here are the same data not subsetted but with the subset region overlayed on top.

rgbRast <- create_stack(file=f,
                         bands=bands,
                         epsg=epsg,
                         subset=FALSE)

plotRGB(rgbRast,
        stretch="lin")

plot(clip.extent,
     add=T,
     border="yellow",
     lwd=3)

View Spectra Within Clipped Raster

# array containing the index dimensions to slice
H5close()

subset.h5 <- h5read(f, "Reflectance",
                    index=list(index.bounds[1]:index.bounds[2],
                    					 index.bounds[3]:index.bounds[4],
                    					 1:426)) # the column, row

final.spectra <- data.frame(apply(subset.h5,
              MARGIN = c(3), # take the mean value for each z value
              mean)) # grab the mean value in the z dimension
final.spectra$wavelength <- wavelengths

# scale the data
names(final.spectra) <- c("reflectance","wavelength")
final.spectra$reflectance <- final.spectra$reflectance / 10000


# plot the data
qplot(x=final.spectra$wavelength,
      y=final.spectra$reflectance,
      xlab="Wavelength (nm)",
      ylab="Reflectance",
      main="Mean Spectral Signature for Clip Region")


Get Lesson Code

(some browsers may require you to right click.)