<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Band Stacking, RGB & False Color Images Using Functions
Objectives
In this lesson, we will cover
Step1:  Use pwd (print working directory) to see what directory we're currently in.
Step2:  For this tutorial to run we need to be in the hyperspectral directory. If you're not there, you can use cd (change directory) to get there.
Step3:  Next we'll import some hyperspectral functions specific to the AOP data, from the neon_aop_hyperspectral.py module. We encourage you to look through the code in the module on your own time to gain a better understanding of how to read in the NEON hyperspectral hdf5 format. We'll start by reading in an h5 reflectance file using the aop_h5refl2array function, which pulls out the hyperspectral data cube (m x n x bands) and corresponding metadata that contains relevant information such as the spatial referencing. 
This function also cleans the reflectance data
Step4:  To see all the functions and modules we've imported use %whos
Step5:  As with standard python functions, we've included docstrings in these functions so that you can use help or ? to display details about how to use the function.
Step6:  Now that we have an idea of how this function works, let's try it out. First, set the variable serc_h5_tile to be reflectance data, specifying the relative path, or alterantively using r'Path\To\File.h5' or  os.path.join to create the full path to the data file. Note that if you want to run this notebook later on a different reflectance tile, you just have to change this variable.
Step7:  Now that we've specified our reflectance tile, we can call aop_h5refl2array to read in the reflectance tile as a Python Numpy array called refl , and the associated metadata into a dictionary metadata.
Step8:  We can use the shape method to see the dimensions of the array we read in. NEON tiles are (1000 x 1000 x # of bands), the number of bands may vary depending on the hyperspectral sensor used, but should be in the vicinity of 426.
Step9:  Plot a single band of reflectance data
Next we'll use the function plot_aop_refl to plot a single band of reflectance data. Read the Parameters section of the docstring to understand the required inputs & data type for each of these; only the band and spatial extent are required inputs, the rest are optional inputs that, if specified, allow you to set the range color values, specify the axis, add a title, colorbar, colorbar title, and change the colormap (default is to plot in greyscale).
Step10:  Now that we have loaded this function, let's extract a single band from the SERC reflectance array and plot it
Step11:  RGB Plots - Band Stacking
It is often useful to look at several bands together. We can extract and stack three reflectance bands in the red, green, and blue (RGB) spectrums to produce a color image that looks like what we see with our eyes; this is your typical camera image. In the next part of this tutorial, we will learn to stack multiple bands and make a geotif raster from the compilation of these bands. We can see that different combinations of bands allow for different visualizations of the remotely-sensed objects and also conveys useful information about the chemical makeup of the Earth's surface. 
For this exercise, we'll use the function stack_rgb to extract the bands we want to stack. This function uses splicing to extract the nth band from the reflectance array, and then uses the numpy function stack to create a new 3D array (1000 x 1000 x 3) consisting of only the three bands we want to visualize.
Step12:  First we will look at red, green, and blue bands, whos indices are defined below. To confirm that these band indices correspond to wavelengths in the expected portion of the spectrum, we can print out the wavelength values stored in metadata['wavelength']
Step13:  We selected these bands so that they fall within the visible range of the electromagnetic spectrum (400-700 nm)
Step14:  Plot an RGB band combination
Next, we can use the function plot_aop_refl, even though we have more than one band. This function only works for a single or 3-band array, so ensure the array you use has the proper dimensions before using. You do not need to specify the colorlimits as the matplotlib.pyplot automatically scales 3-band arrays to 8-bit color (256).
Step15:  You'll notice that this image is very dark; it is possible to make out some of the features (roads, buildings), but it is not ideal. Since colorlimits don't apply to 3-band images, we have to use some other image processing tools to enhance the visibility of this image. 
Image Processing -- Contrast Stretch & Histogram Equalization
We can also try out some image processing routines to better visualize the reflectance data using the ski-image package. 
Histogram equalization is a method in image processing of contrast adjustment using the image's histogram. Stretching the histogram can improve the contrast of a displayed image by eliminating very high or low reflectance values that skew the display of the image. 
<p>
<img src="http
Step16:  False Color Image - Color Infrared (CIR)
We can also create an image from bands outside of the visible spectrum. An image containing one or more bands outside of the visible range is called a false-color image. Here we'll use the green and blue bands as before, but we replace the red band with a near-infrared (NIR) band. 
For more information about non-visible wavelengths, false color images, and some frequently used false-color band combinations, refer to NASA's Earth Observatory page
Step17:  Demo
Step18:  Demo
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import numpy as np
import h5py, os, osr, copy
%matplotlib inline
Explanation: Band Stacking, RGB & False Color Images Using Functions
Objectives
In this lesson, we will cover:
1. pulling in the function aop_h5refl2array to read in a reflectance hdf5 file and convert it a python 3-dimensional array (rows x columns x bands)
2. plotting different band combinations using the functions plot_aop_refl and stack_rgb, and 
3. using widgets for more interactive data visualization. 
Background on RGB & False-Color Images
We can combine any three bands from the NEON reflectance data to make an RGB image that allows us to visualize spectral data collected at the Earth's surface. A natural-color image, made with bands from the red, green, and blue wavelengths looks close to what we would see with the naked eye, or what we would see in a typical photograph. We can also choose band combinations from other wavelengths, and map them to the red, blue, and green colors to highlight different spectral features. A false-color image is made with one or more bands from a non-visible portion of the electromagnetic spectrum that are mapped to red, green, and blue colors. These images can display other information about the landscape that is not easily seen with a natural-color image. 
The NASA Goddard Media Studio video Peeling Back Landsat's Layers of Data gives a good overview of natural and false color band combinations. Note that Landsat collects information from 11 wavelength bands, while NEON AOP hyperspectral data collects information from 426 bands!
<a href="https://www.youtube.com/watch?v=YP0et8l_bvY" target="_blank"><img src="./notebook_images//NASA_Goddard_RGB_Images.png" 
alt="IMAGE ALT TEXT HERE" width="400" height="300" border="10" /></a>
https://www.youtube.com/watch?v=YP0et8l_bvY
https://svs.gsfc.nasa.gov//vis/a010000/a011400/a011491/index.html
Further Reading: Check out the NASA Earth Observatory Article on "How to Interpret a False-Color Satellite Image"
https://earthobservatory.nasa.gov/Features/FalseColor/
Before we get started, let's import our standard packages, and set up our plot preferences to inline so our figures show up in the notebook as we go:
End of explanation
pwd
Explanation: Use pwd (print working directory) to see what directory we're currently in.
End of explanation
cd hyperspectral
Explanation: For this tutorial to run we need to be in the hyperspectral directory. If you're not there, you can use cd (change directory) to get there.
End of explanation
from neon_aop_hyperspectral import *
Explanation: Next we'll import some hyperspectral functions specific to the AOP data, from the neon_aop_hyperspectral.py module. We encourage you to look through the code in the module on your own time to gain a better understanding of how to read in the NEON hyperspectral hdf5 format. We'll start by reading in an h5 reflectance file using the aop_h5refl2array function, which pulls out the hyperspectral data cube (m x n x bands) and corresponding metadata that contains relevant information such as the spatial referencing. 
This function also cleans the reflectance data: it sets any no data values within the reflectance tile to NaN (not a number) and applies the reflectance scale factor so the final array that is returned represents unitless scaled reflectance, with values ranging between 0 and 1 (0-100%).
End of explanation
%whos
Explanation: To see all the functions and modules we've imported use %whos
End of explanation
help(aop_h5refl2array)
Explanation: As with standard python functions, we've included docstrings in these functions so that you can use help or ? to display details about how to use the function.
End of explanation
#h5_tile = ('../../../data/NEON_D16_MCRA_DP3_566300_4901000_reflectance.h5') 
h5_tile = ('../../../data/2019_JORN/Refl/NEON_D14_JORN_DP3_324000_3620000_reflectance.h5') 
Explanation: Now that we have an idea of how this function works, let's try it out. First, set the variable serc_h5_tile to be reflectance data, specifying the relative path, or alterantively using r'Path\To\File.h5' or  os.path.join to create the full path to the data file. Note that if you want to run this notebook later on a different reflectance tile, you just have to change this variable.
End of explanation
refl,metadata = aop_h5refl2array(h5_tile)
Explanation: Now that we've specified our reflectance tile, we can call aop_h5refl2array to read in the reflectance tile as a Python Numpy array called refl , and the associated metadata into a dictionary metadata.
End of explanation
refl.shape
list(metadata.items())
Explanation: We can use the shape method to see the dimensions of the array we read in. NEON tiles are (1000 x 1000 x # of bands), the number of bands may vary depending on the hyperspectral sensor used, but should be in the vicinity of 426.
End of explanation
plot_aop_refl?
Explanation: Plot a single band of reflectance data
Next we'll use the function plot_aop_refl to plot a single band of reflectance data. Read the Parameters section of the docstring to understand the required inputs & data type for each of these; only the band and spatial extent are required inputs, the rest are optional inputs that, if specified, allow you to set the range color values, specify the axis, add a title, colorbar, colorbar title, and change the colormap (default is to plot in greyscale).
End of explanation
b56 = refl[:,:,55]
plot_aop_refl(b56,
              metadata['spatial extent'],
              colorlimit=(0,0.3),
              title='Band 56 Reflectance',
              cmap_title='Reflectance',
              colormap='Greys_r') 
Explanation: Now that we have loaded this function, let's extract a single band from the SERC reflectance array and plot it:
End of explanation
stack_rgb?
Explanation: RGB Plots - Band Stacking
It is often useful to look at several bands together. We can extract and stack three reflectance bands in the red, green, and blue (RGB) spectrums to produce a color image that looks like what we see with our eyes; this is your typical camera image. In the next part of this tutorial, we will learn to stack multiple bands and make a geotif raster from the compilation of these bands. We can see that different combinations of bands allow for different visualizations of the remotely-sensed objects and also conveys useful information about the chemical makeup of the Earth's surface. 
For this exercise, we'll use the function stack_rgb to extract the bands we want to stack. This function uses splicing to extract the nth band from the reflectance array, and then uses the numpy function stack to create a new 3D array (1000 x 1000 x 3) consisting of only the three bands we want to visualize.
End of explanation
rgb_bands = (58,34,19)
print('Band 58 Center Wavelength = %.2f' %(metadata['wavelength'][57]),'nm')
print('Band 33 Center Wavelength = %.2f' %(metadata['wavelength'][33]),'nm')
print('Band 19 Center Wavelength = %.2f' %(metadata['wavelength'][18]),'nm')
Explanation: First we will look at red, green, and blue bands, whos indices are defined below. To confirm that these band indices correspond to wavelengths in the expected portion of the spectrum, we can print out the wavelength values stored in metadata['wavelength']:
End of explanation
rgb = stack_rgb(refl,rgb_bands)
rgb.shape
Explanation: We selected these bands so that they fall within the visible range of the electromagnetic spectrum (400-700 nm):
<img src="./notebook_images/NEON_AOP_RGB_Band_Wavelengths.png" style="width: 800px;">
Below we use stack_rgb to create an RGB array. Check that the dimensions of this array are as expected.
- TIP: checking the shape of arrays with the .shape method is a good habit to get into when creating your own workflows, and can be a handy tool for troubleshooting.
End of explanation
plot_aop_refl(rgb,
              metadata['spatial extent'],
              title='RGB Image',
              cbar='off') 
Explanation: Plot an RGB band combination
Next, we can use the function plot_aop_refl, even though we have more than one band. This function only works for a single or 3-band array, so ensure the array you use has the proper dimensions before using. You do not need to specify the colorlimits as the matplotlib.pyplot automatically scales 3-band arrays to 8-bit color (256).
End of explanation
plot_aop_rgb?
plot_aop_rgb(rgb,
             metadata['spatial extent'],
             ls_pct=1,
             plot_title = 'RGB')
Explanation: You'll notice that this image is very dark; it is possible to make out some of the features (roads, buildings), but it is not ideal. Since colorlimits don't apply to 3-band images, we have to use some other image processing tools to enhance the visibility of this image. 
Image Processing -- Contrast Stretch & Histogram Equalization
We can also try out some image processing routines to better visualize the reflectance data using the ski-image package. 
Histogram equalization is a method in image processing of contrast adjustment using the image's histogram. Stretching the histogram can improve the contrast of a displayed image by eliminating very high or low reflectance values that skew the display of the image. 
<p>
<img src="http://opencv-python-tutroals.readthedocs.io/en/latest/_images/histogram_equalization.png" style="width: 300px;"/>
</p>
Let's see what the image looks like using a 5% linear contrast stretch using the skimage module's exposure function.
End of explanation
CIRbands = (90,34,19)
print('Band 90 Center Wavelength = %.2f' %(metadata['wavelength'][89]),'nm')
print('Band 34 Center Wavelength = %.2f' %(metadata['wavelength'][33]),'nm')
print('Band 19 Center Wavelength = %.2f' %(metadata['wavelength'][18]),'nm')
cir = stack_rgb(refl,CIRbands)
plot_aop_rgb(cir,
             metadata['spatial extent'],
             ls_pct=2,
             plot_title='Color Infrared Image')
Explanation: False Color Image - Color Infrared (CIR)
We can also create an image from bands outside of the visible spectrum. An image containing one or more bands outside of the visible range is called a false-color image. Here we'll use the green and blue bands as before, but we replace the red band with a near-infrared (NIR) band. 
For more information about non-visible wavelengths, false color images, and some frequently used false-color band combinations, refer to NASA's Earth Observatory page:
https://earthobservatory.nasa.gov/Features/FalseColor/
End of explanation
from ipywidgets import *
array = refl.copy()
metadata = metadata.copy()
def RGBplot_widget(R,G,B):
    
    #Pre-allocate array  size
    rgbArray = np.zeros((array.shape[0],array.shape[1],3), 'uint8')
    
    Rband = array[:,:,R-1].astype(np.float)
    #Rband_clean = clean_band(Rband,Refl_md)
    
    Gband = array[:,:,G-1].astype(np.float)
    #Gband_clean = clean_band(Gband,Refl_md)
    
    Bband = array[:,:,B-1].astype(np.float)
    #Bband_clean = clean_band(Bband,Refl_md)
    
    rgbArray[..., 0] = Rband*256
    rgbArray[..., 1] = Gband*256
    rgbArray[..., 2] = Bband*256
    
    # Apply Adaptive Histogram Equalization to Improve Contrast:
    
    img_nonan = np.ma.masked_invalid(rgbArray) #first mask the image 
    img_adapteq = exposure.equalize_adapthist(img_nonan, clip_limit=0.10)
    
    plot = plt.imshow(img_adapteq,extent=metadata['spatial extent']); 
    plt.title('Bands: \nR:' + str(R) + ' (' + str(round(metadata['wavelength'][R-1])) +'nm)'
              + '\n G:' + str(G) + ' (' + str(round(metadata['wavelength'][G-1])) + 'nm)'
              + '\n B:' + str(B) + ' (' + str(round(metadata['wavelength'][B-1])) + 'nm)'); 
    ax = plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') 
    rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) 
    
interact(RGBplot_widget, R=(1,426,1), G=(1,426,1), B=(1,426,1));
Explanation: Demo: Exploring Band Combinations Interactively
Now that we have made a couple different band combinations, we can demo a Python widget to explore different combinations of bands in the visible and non-visible portions of the spectrum.
End of explanation
rgbArray = rgb.copy()
def linearStretch(percent):
    pLow, pHigh = np.percentile(rgbArray[~np.isnan(rgbArray)], (percent,100-percent))
    img_rescale = exposure.rescale_intensity(rgbArray, in_range=(pLow,pHigh))
    plt.imshow(img_rescale,extent=metadata['spatial extent'])
    plt.title('SERC RGB \n Linear ' + str(percent) + '% Contrast Stretch'); 
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False, style='plain') 
    rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) 
interact(linearStretch,percent=(0,20,1));
Explanation: Demo: Interactive Linear Stretch & Equalization
Here is another widget to play around with, demonstrating how to interactively visualize linear contrast stretches with a variable percent.
End of explanation
<END_TASK>