updates to use full reflectance tile, download from neonUtilities, te…

…rra package, etc.

tutorials/R/AOP/Hyperspectral/Select-Pixels-Compare-Spectral-Signatures/Select-Pixels-Compare-Spectral-Signatures.R

@@ -1,70 +1,86 @@
-## ----load-libraries, message=FALSE, warning=FALSE--------------------
+## ----load-libraries, message=FALSE, warning=FALSE-------------------------------------------------------------------------------------------
 
-# Load required packages
+# load required packages
 library(rhdf5)
 library(reshape2)
-library(raster)
+library(terra)
 library(plyr)
 library(ggplot2)
+library(grDevices)
 
-# set working directory to ensure R can find the file we wish to import and where
-# we want to save our files. Be sure to move the download into your working directory!
-wd <- "~/Documents/data/" #This will depend on your local environment
+# set working directory, you can change this if desired
+wd <- "~/data/" 
 setwd(wd)
 
+
+## ----download-h5, eval=FALSE----------------------------------------------------------------------------------------------------------------
+## byTileAOP(dpID = 'DP3.30006.001',
+##           site = 'SJER',
+##           year = '2021',
+##           easting = 257500,
+##           northing = 4112500,
+##           savepath = wd)
+
+
+## ----read-h5--------------------------------------------------------------------------------------------------------------------------------
 # define filepath to the hyperspectral dataset
-fhs <- paste0(wd,"NEON_hyperspectral_tutorial_example_subset.h5")
+h5_file <- paste0(wd,"DP3.30006.001/neon-aop-products/2021/FullSite/D17/2021_SJER_5/L3/Spectrometer/Reflectance/NEON_D17_SJER_DP3_257000_4112000_reflectance.h5")
 
 # read in the wavelength information from the HDF5 file
-wavelengths <- h5read(fhs,"/SJER/Reflectance/Metadata/Spectral_Data/Wavelength")
+wavelengths <- h5read(h5_file,"/SJER/Reflectance/Metadata/Spectral_Data/Wavelength")
 
 # grab scale factor from the Reflectance attributes
-reflInfo <- h5readAttributes(fhs,"/SJER/Reflectance/Reflectance_Data" )
+reflInfo <- h5readAttributes(h5_file,"/SJER/Reflectance/Reflectance_Data" )
 
 scaleFact <- reflInfo$Scale_Factor
 
 
+## ----read-in-RGB-and-plot, fig.cap="RGB image of a portion of the SJER field site using 3 bands from the raster stack. Brightness values have been stretched using the stretch argument to produce a natural looking image."----
 
-## ----read-in-RGB-and-plot, fig.cap="RGB image of a portion of the SJER field site using 3 bands from the raster stack. Brightness values have been stretched using the stretch argument to produce a natural looking image. At the top right of the image, there is dark, brakish water. Scattered throughout the image, there are several trees. At the center of the image, there is a baseball field, with low grass. At the bottom left of the image, there is a parking lot and some buildings with highly reflective surfaces, and adjacent to it is a section of a gravel lot."----
+# read in RGB image as a 'stack' rather than a plain 'raster'
+rgbStack <- rast(paste0(wd,"NEON_hyperspectral_tutorial_example_RGB_image.tif"))
 
-# Read in RGB image as a 'stack' rather than a plain 'raster'
-rgbStack <- stack(paste0(wd,"NEON_hyperspectral_tutorial_example_RGB_stack_image.tif"))
-
-# Plot as RGB image
+# plot as RGB image, with a linear stretch
 plotRGB(rgbStack,
         r=1,g=2,b=3, scale=300, 
         stretch = "lin")
 
 
+## ----dev-new, eval=FALSE--------------------------------------------------------------------------------------------------------------------
+## dev.new(noRStudioGD = TRUE)
+
 
-## ----click-to-select, eval=FALSE, comment=NA-------------------------
+## ----click-to-select, eval=FALSE, comment=NA------------------------------------------------------------------------------------------------
 
 # change plotting parameters to better see the points and numbers generated from clicking
-par(col="red", cex=3)
+par(col="red", cex=2)
+
+# use a histogram stretch in order to provide more contrast for selecting pixels
+plotRGB(rgbStack, r=1, g=2, b=3, scale=300, stretch = "hist") 
 
 # use the 'click' function
-c <- click(rgbStack, id=T, xy=T, cell=T, type="p", pch=16, col="magenta", col.lab="red")
+c <- click(rgbStack, n = 7, id=TRUE, xy=TRUE, cell=TRUE, type="p", pch=16, col="red", col.lab="red")
 
 
 
 
 
-## ----convert-cell-to-row-column--------------------------------------
+## ----convert-cell-to-row-column-------------------------------------------------------------------------------------------------------------
 # convert raster cell number into row and column (used to extract spectral signature below)
 c$row <- c$cell%/%nrow(rgbStack)+1 # add 1 because R is 1-indexed
 c$col <- c$cell%%ncol(rgbStack)
 
 
-## ----extract-spectral-signaures--------------------------------------
+## ----extract-spectral-signaures-------------------------------------------------------------------------------------------------------------
 
 # create a new dataframe from the band wavelengths so that we can add
 # the reflectance values for each cover type
-Pixel_df <- as.data.frame(wavelengths)
+pixel_df <- as.data.frame(wavelengths)
 
 # loop through each of the cells that we selected
 for(i in 1:length(c$cell)){
 # extract Spectra from a single pixel
-aPixel <- h5read(fhs,"/SJER/Reflectance/Reflectance_Data",
+aPixel <- h5read(h5_file,"/SJER/Reflectance/Reflectance_Data",
                  index=list(NULL,c$col[i],c$row[i]))
 
 # scale reflectance values from 0-1
@@ -76,73 +92,73 @@ b <- adply(aPixel,c(1))
 # rename the column that we just created
 names(b)[2] <- paste0("Point_",i)
 
-# add reflectance values for this pixel to our combined data.frame called Pixel_df
-Pixel_df <- cbind(Pixel_df,b[2])
+# add reflectance values for this pixel to our combined data.frame called pixel_df
+pixel_df <- cbind(pixel_df,b[2])
 }
 
 
 
-## ----plot-spectral-signatures, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the five different land cover types: Field, Tree, Roof, Soil, and Water. On the x-axis is wavelength in nanometers and on the y-axis is reflectance values."----
+## ----plot-spectral-signatures, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the six different land cover types: Water, Tree, Grass, Soil, Roof, and Road. The x-axis is wavelength in nanometers and the y-axis is reflectance."----
 # Use the melt() function to reshape the dataframe into a format that ggplot prefers
-Pixel.melt <- reshape2::melt(Pixel_df, id.vars = "wavelengths", value.name = "Reflectance")
+pixel.melt <- reshape2::melt(pixel_df, id.vars = "wavelengths", value.name = "Reflectance")
 
 # Now, let's plot some spectral signatures!
 ggplot()+
-  geom_line(data = Pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
-  scale_colour_manual(values = c("green2", "green4", "grey50","tan4","blue3"),
-                      labels = c("Field", "Tree", "Roof","Soil","Water"))+
+  geom_line(data = pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
+  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
+                      labels = c("Water", "Tree", "Grass","Soil","Roof","Road"))+
   labs(color = "Cover Type")+
   ggtitle("Land cover spectral signatures")+
   theme(plot.title = element_text(hjust = 0.5, size=20))+
   xlab("Wavelength")
 
 
-## ----mask-atmospheric-absorbtion-bands, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the five different land cover types: Field, Tree, Roof, Soil, and Water. Added to the plot are two rectangles in regions near 1400nm and 1850nm where the reflectance measurements are obscured by atmospheric absorption. On the x-axis is wavelength in nanometers and on the y-axis is reflectance values."----
+## ----mask-atmospheric-absorption-bands, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the six different land cover types: Water, Tree, Grass, Soil, Roof, and Road. Added to the plot are two greyed-out rectangles in regions near 1400nm and 1850nm where the reflectance measurements are obscured by atmospheric absorption. The x-axis is wavelength in nanometers and the y-axis is reflectance."----
 
-# grab Reflectance metadata (which contains absorption band limits)
+# grab reflectance metadata (which contains absorption band limits)
 reflMetadata <- h5readAttributes(fhs,"/SJER/Reflectance" )
 
 ab1 <- reflMetadata$Band_Window_1_Nanometers
 ab2 <- reflMetadata$Band_Window_2_Nanometers
 
 # Plot spectral signatures again with rectangles showing the absorption bands
 ggplot()+
-  geom_line(data = Pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
-  geom_rect(mapping=aes(ymin=min(Pixel.melt$Reflectance),ymax=max(Pixel.melt$Reflectance), xmin=ab1[1], xmax=ab1[2]), color="black", fill="grey40", alpha=0.8)+
-  geom_rect(mapping=aes(ymin=min(Pixel.melt$Reflectance),ymax=max(Pixel.melt$Reflectance), xmin=ab2[1], xmax=ab2[2]), color="black", fill="grey40", alpha=0.8)+
-  scale_colour_manual(values = c("green2", "green4", "grey50","tan4","blue3"),
-                      labels = c("Field", "Tree", "Roof","Soil","Water"))+
+  geom_line(data = pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
+  geom_rect(mapping=aes(ymin=min(pixel.melt$Reflectance),ymax=max(pixel.melt$Reflectance), xmin=ab1[1], xmax=ab1[2]), color="black", fill="grey40", alpha=0.8)+
+  geom_rect(mapping=aes(ymin=min(pixel.melt$Reflectance),ymax=max(pixel.melt$Reflectance), xmin=ab2[1], xmax=ab2[2]), color="black", fill="grey40", alpha=0.8)+
+  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
+                      labels = c("Water","Tree","Grass","Soil","Roof","Road"))+
   labs(color = "Cover Type")+
   ggtitle("Land cover spectral signatures")+
   theme(plot.title = element_text(hjust = 0.5, size=20))+
   xlab("Wavelength")
 
 
-## ----remove-absorbtion-band-reflectances, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the five different land cover types: Field, Tree, Roof, Soil, and Water. Values falling within the two rectangles from the previous image have been set to NA and ommited from the plot. On the x-axis is wavelength in nanometers and on the y-axis is reflectance values."----
+## ----remove-absorption-band-reflectances, fig.width=9, fig.height=6, fig.cap="Plot of spectral signatures for the six different land cover types. Values falling within the two rectangles from the previous image have been set to NA and ommited from the plot. The x-axis is wavelength in nanometers and the y-axis is reflectance."----
 
 # Duplicate the spectral signatures into a new data.frame
-Pixel.melt.masked <- Pixel.melt
+pixel.melt.masked <- pixel.melt
 
-# Mask out all values within each of the two atmospheric absorbtion bands
-Pixel.melt.masked[Pixel.melt.masked$wavelengths>ab1[1]&Pixel.melt.masked$wavelengths<ab1[2],]$Reflectance <- NA
-Pixel.melt.masked[Pixel.melt.masked$wavelengths>ab2[1]&Pixel.melt.masked$wavelengths<ab2[2],]$Reflectance <- NA
+# Mask out all values within each of the two atmospheric absorption bands
+pixel.melt.masked[pixel.melt.masked$wavelengths>ab1[1]&pixel.melt.masked$wavelengths<ab1[2],]$Reflectance <- NA
+pixel.melt.masked[pixel.melt.masked$wavelengths>ab2[1]&pixel.melt.masked$wavelengths<ab2[2],]$Reflectance <- NA
 
 # Plot the masked spectral signatures
 ggplot()+
-  geom_line(data = Pixel.melt.masked, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
-  scale_colour_manual(values = c("green2", "green4", "grey50","tan4","blue3"),
-                      labels = c("Field", "Tree", "Roof", "Soil", "Water"))+
+  geom_line(data = pixel.melt.masked, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
+  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
+                      labels = c("Water","Tree","Grass","Soil","Roof","Road"))+
   labs(color = "Cover Type")+
   ggtitle("Land cover spectral signatures")+
   theme(plot.title = element_text(hjust = 0.5, size=20))+
   xlab("Wavelength")
 
 
 
-## ----challenge-answer, echo=FALSE, eval=FALSE------------------------
+## ----challenge-answer, echo=FALSE, eval=FALSE-----------------------------------------------------------------------------------------------
 ## 
 ## # Challenge Answers - These challenge problems will depend on the specific
-## # pixels that you select, but here we can answer these questions in general.
+## # pixels that you select, but here we answer these questions generally.
 ## 
 ## # 1. Each vegetation class will likely have slightly different spectral signatures,
 ## # mostly distinguished by the amplitude of the near-IR bands. As we saw in this