OpenHDF5-Cattau.Rmd

---
title: "openHDF5"
author: "Megan Cattau"
date: "June 20, 2016"
output: html_document
---

```{r load-libraries}
# load the required packages
library(raster)
library(rhdf5)
library(rgdal)
```

<br>
We're going to open the reflectance dataset from the Teakettle site and look at the structure
```{r select-file}
# if you type "../" then hit tab, the file drop-down options will open automatically
f<-"../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# view H5 structure to explore it
h5ls(f)
```

<br>
Import spatial information - view the attributes (metadata) for the map info dataset within the NEON H5 file.
```{r import-spatial-info}
# tell R where the data are located
# in HDF viewer, the "map info" group tells you the coordinate reference system and the coordinate loaction of the upper left corner of the file and spatial resolution

mapInfo<-h5read(f,
                "map info", # one can find the name of this data in HDF view
                read.attributes = TRUE)
mapInfo
```

<br>
We're going to import the reflectance metadata and define the scale factor and no data value
```{r get-reflectance-metadata}

# read in reflectance data attributes
reflInfo<-h5readAttributes(f,
                           name = "Reflectance")
reflInfo

# define scale factor and no data value as objects
scaleFactor<-reflInfo$`Scale Factor`
noDataValue<-as.numeric(reflInfo$`data ignore value`)
```

<br>
Import data dimensions
```{r import-dims}
# open file for viewing. This is a direct connection to the file, but it doesn't read it in. We don't want to read in a huge data cube into memory.
fid<-H5Fopen(f)

# Use this to open the reflectance dataset itself
did<-H5Dopen(fid, "Reflectance")
did
# this has columns first rather than rows

# grab the dimensions of the object
# you could also extract dimensions from H5ls, but this creates a text object and then you have to split it out
# H5ls useful if you want to start looping through multiple files
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims
# Columns are the FIRST dimension, then rows. wavelength is the THIRD dimension rather than the first. 

# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)
```

<br>
Read Reflectance Data

Once we know the dimensions of the data, we can more efficiently slice out chunks or subsets of our data out. The power of HDF5 is that it allows us to store large heterogeneous data. However, we can quickly and efficiently access those data through “slicing” or extracting subsets of the data.

```{r read-reflectance-data}
# Extract or "slice" data (one layer of datacube) for band 56 from the HDF5 file
b56<- h5read(f,"Reflectance", index=list(1:dims[1],1:dims[2],56))

# note the data come in as an array
class(b56)
```

<br>
Convert to Matrix and plot the data
```{r convert-to-matrix-and-view}
# Convert from array to matrix so we can plot and convert to a raster
b56 <- b56[,,1]

# plot the data
image(b56)

# stretch it
image(log(b56),
			main="Band 56 with log Transformation")
# view distribution of reflectance data

# force non scientific notation
options("scipen"=100, "digits"=4)

hist(b56,
     col="springgreen",
     main="Distribution of Reflectance Values \nBand 56")
```

<br>
Data Clean-up
Set the no data value to 15,000.
Apply the scale factor to the data (10,000).
```{r clean-up}
# extract no data value from the attributes
noDataVal <- as.integer(reflInfo$`data ignore value`)
# set all reflectance values = 15,000 to NA
b56[b56 == noDataVal] <- NA

# Extract the scale factor as an object
scaleFactor <- reflInfo$`Scale Factor`

# divide all values in our B56 object by the scale factor to get a range of
# reflectance values between 0-1 (the valid range)
b56 <- b56/scaleFactor

# view distribution of reflectance values
hist(b56,
     main="Distribution with NoData Value Considered\nData Scaled")
```

<br>
Unflip the data
```{r unflip-data}
b56<-t(b56)
image(log(b56), main="Band 56\nTransposed Values")
```

<br>
Create Spatial Extent
```{r create-spatial-extent}

mapInfo<-unlist(strsplit(mapInfo, ","))

# X,Y left corner coordinate - change to numeric
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])

# get x and y resolution
xres <- as.numeric(mapInfo[6])
yres <- as.numeric(mapInfo[7])

# finally calculate the xMax value and the yMin value from the dimensions
# we grabbed above. The xMax is the left corner + number of columns* resolution.
xMax <- xMin + (dims[1]*xres)
yMin <- yMax - (dims[2]*yres)

# Define the raster extent
rasExt <- extent(xMin, xMax,yMin,yMax)

# Create a raster and assign it it's spatial extent
b56r <- raster(b56,
               crs=CRS("+init=epsg:32611"))
# assign CRS
extent(b56r) <- rasExt

# view raster object attributes
b56r

# plot the new image
plot(b56r, main="Raster for Lower Teakettle \nBand 56")
```

<br>
Use NEON functions for a way easier way to do what we just did!
```{r use-NEON-functions}
# install.packages("devtools")
library(devtools)
# install_github("lwasser/neon-aop-package/neonAOP")

library(neonAOP)
# open band function 
b55<-open_band(f,
               bandNum=55,
               epsg=32611)
plot(b55)

# import several bands
bands<-c(58, 34, 19)

# create a raster stack
RGBStack<-create_stack(f, 
                       bands = bands,
                       epsg = 32611)
plot(RGBStack) # will plot 3 bands separately
plotRGB(RGBStack,
        stretch="lin")

# epsg codes http://spatialreference.org/
```

Testing second time

<br>