In this tutorial, we will use the Spectral Python (SPy) package to run a KMeans unsupervised classification algorithm and then we will run Principal Component Analysis to reduce data dimensionality.

In this tutorial, we will use the Spectral Python (SPy) package to run KMeans unsupervised classification algorithm as well as Principal Component Analysis (PCA).

To learn more about the Spectral Python packages read:

• Spectral Python User Guide.
• Spectral Python Unsupervised Classification.

KMeans Clustering

KMeans is an iterative clustering algorithm used to classify unsupervised data (eg. data without a training set) into a specified number of groups. The algorithm begins with an initial set of randomly determined cluster centers. Each pixel in the image is then assigned to the nearest cluster center (using distance in N-space as the distance metric) and each cluster center is then re-computed as the centroid of all pixels assigned to the cluster. This process repeats until a desired stopping criterion is reached (e.g. max number of iterations).

Read more on KMeans clustering from Spectral Python.

To visualize how the algorithm works, it's easier look at a 2D data set. In the example below, watch how the cluster centers shift with progressive iterations,

Principal Component Analysis (PCA) - Dimensionality Reduction

Many of the bands within hyperspectral images are often strongly correlated. The principal components transformation represents a linear transformation of the original image bands to a set of new, uncorrelated features. These new features correspond to the eigenvectors of the image covariance matrix, where the associated eigenvalue represents the variance in the direction of the eigenvector. A very large percentage of the image variance can be captured in a relatively small number of principal components (compared to the original number of bands).

Read more about PCA with Spectral Python.

Let's get started! First, import the required packages.

First, import the required packages and set display preferences:

import h5py
import matplotlib
import neonutilities as nu
import numpy as np
import os
import requests
from spectral import *
from time import time
import dotenv

# Set the data download path, change this path if desired
data_dir = os.path.join(r'C:\data')

For this example, we will download a bidirectional surface reflectance data cube at the SERC site, collected in 2022.

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the python-dotenv package as follows, optionally specifying the path to the .env file in load_dotenv().

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

nu.by_tile_aop(dpid='DP3.30006.002',
               site='SERC',
               year='2022',
               easting=368005,
               northing=4306005,
               include_provisional=True,
               token=token,
               savepath=os.path.join(data_dir)) # save to the home directory under a 'data' subfolder

Provisional NEON data are included. To exclude provisional data, use input parameter include_provisional=False.


Continuing will download 2 NEON data files totaling approximately 659.9 MB. Do you want to proceed? (y/n)  y


Downloading 2 NEON data files totaling approximately 659.9 MB

  0%|                                                                                                                                 | 0/2 [00:00<?, ?it/s]C:\Users\bhass\AppData\Roaming\Python\Python313\site-packages\neonutilities\helper_mods\api_helpers.py:790: UserWarning: Filepaths on Windows are limited to 260 characters. Attempting to download a filepath that is 940 characters long. Set the working or savepath directory to be closer to the root directory or enable long path support in Windows.
  warnings.warn(
C:\Users\bhass\AppData\Roaming\Python\Python313\site-packages\neonutilities\helper_mods\api_helpers.py:790: UserWarning: Filepaths on Windows are limited to 260 characters. Attempting to download a filepath that is 935 characters long. Set the working or savepath directory to be closer to the root directory or enable long path support in Windows.
  warnings.warn(
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 117.22it/s]

Let's see what data were downloaded.

# iterate over directory recursively to show path of downloaded h5 file
for root, dirs, files in os.walk(data_dir):
    for name in files:
        if name.endswith('.h5'):
            h5_tile = os.path.join(root, name)
            print(h5_tile)  # printing file name

C:\data\DP3.30006.002\neon-aop-provisional-products\2022\FullSite\D02\2022_SERC_6\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5
C:\data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

# function to download data stored on the internet in a public url to a local file
def download_url(url,download_dir):
    if not os.path.isdir(download_dir):
        os.makedirs(download_dir)
    filename = url.split('/')[-1]
    r = requests.get(url, allow_redirects=True)
    file_object = open(os.path.join(download_dir,filename),'wb')
    file_object.write(r.content)

module_url = "https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/tutorials/Python/AOP/aop_python_modules/neon_aop_hyperspectral.py"
download_url(module_url,'../python_modules')
# os.listdir('../python_modules') #optionally show the contents of this directory to confirm the file downloaded

sys.path.insert(0, '../python_modules')
# import the neon_aop_hyperspectral module, the semicolon supresses an empty plot from displaying
import neon_aop_hyperspectral as neon_hs;

# read in the reflectance data using the aop_h5refl2array function, this may also take a bit of time
start_time = time()
refl, refl_metadata, wavelengths = neon_hs.aop_h5refl2array(h5_tile,'Reflectance')
print("--- It took %s seconds to read in the data ---" % round((time() - start_time),0))

Reading in  C:\data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5
--- It took 15.0 seconds to read in the data ---

The next few cells show how you can look at the contents, values, and dimensions of the refl_metadata, wavelengths, and refl variables, respectively.

refl_metadata

{'shape': (1000, 1000, 426),
 'no_data_value': -9999.0,
 'scale_factor': 10000.0,
 'bad_band_window1': array([1340, 1445], dtype=int32),
 'bad_band_window2': array([1790, 1955], dtype=int32),
 'projection': b'+proj=UTM +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs',
 'EPSG': 32618,
 'res': {'pixelWidth': 1.0, 'pixelHeight': 1.0},
 'extent': (368000.0, 369000.0, 4306000.0, 4307000.0),
 'ext_dict': {'xMin': 368000.0,
  'xMax': 369000.0,
  'yMin': 4306000.0,
  'yMax': 4307000.0},
 'source': 'C:\\data\\DP3.30006.002\\neon-aop-provisional-products\\2025\\FullSite\\D02\\2025_SERC_7\\L3\\Spectrometer\\Reflectance\\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5'}

print('First and last 5 center wavelengths, in nm:')
print(wavelengths[:5])
print(wavelengths[-5:])

First and last 5 center wavelengths, in nm:
[381.858398 386.868896 391.879395 396.889893 401.900391]
[2491.281494 2496.291992 2501.30249  2506.312988 2511.323486]

refl.shape

(1000, 1000, 426)

Next let's define a function to clean and subset the data.

def clean_neon_refl_data(data, metadata, wavelengths, subset_factor=1):
    """Clean h5 reflectance data and metadata
    1. set data ignore value (-9999) to NaN
    2. apply reflectance scale factor (10000)
    3. remove bad bands (water vapor band windows + last 10 bands): 
        Band_Window_1_Nanometers = 1340, 1445
        Band_Window_2_Nanometers = 1790, 1955
    4. if subset_factor, subset by that factor
    """
    
    # use copy so original data and metadata doesn't change
    data_clean = data.copy().astype(float)
    metadata_clean = metadata.copy()
    
    #set data ignore value (-9999) to NaN:
    if metadata['no_data_value'] in data:
        nodata_ind = np.where(data_clean==metadata['no_data_value'])
        data_clean[nodata_ind]=np.nan 
    
    #apply reflectance scale factor (divide by 10000)
    data_clean = data_clean/metadata['scale_factor']
    
    #remove bad bands 
    #1. define indices corresponding to min/max center wavelength for each bad band window:
    bb1_ind0 = np.max(np.where(np.asarray(wavelengths<float(metadata['bad_band_window1'][0]))))
    bb1_ind1 = np.min(np.where(np.asarray(wavelengths>float(metadata['bad_band_window1'][1]))))

    bb2_ind0 = np.max(np.where(np.asarray(wavelengths<float(metadata['bad_band_window2'][0]))))
    bb2_ind1 = np.min(np.where(np.asarray(wavelengths>float(metadata['bad_band_window2'][1]))))
    bb3_ind0 = len(wavelengths)-15
    
    #define valid band ranges from indices:
    vb1 = list(range(10,bb1_ind0)); 
    vb2 = list(range(bb1_ind1,bb2_ind0))
    vb3 = list(range(bb2_ind1,bb3_ind0))
    # combine them to get a list of the valid bands
    vbs = vb1 + vb2 + vb3
    # subset by subset_factor (if subset_factor = 1 this will return the original valid_bands list)
    valid_bands_subset = vbs[::subset_factor]

    # subset the reflectance data by the valid_bands_subset
    data_clean = data_clean[:,:,valid_bands_subset]

    # subset the wavelengths by the same valid_bands_subset
    wavelengths_clean =[wavelengths[i] for i in valid_bands_subset]
    
    return data_clean, wavelengths_clean

Now use this function to clean and subset the data, using a subset factor of 2 to start.

# clean the data - remove the band bands and subset
start_time = time()
refl_clean, wavelengths_clean = clean_neon_refl_data(refl, refl_metadata, wavelengths, subset_factor=2)
print("--- It took %s seconds to clean and subset the reflectance data ---" % round((time() - start_time),0))

--- It took 8.0 seconds to clean and subset the reflectance data ---

# Look at the dimensions of the data after cleaning:
print('Cleaned Data Dimensions:',refl_clean.shape)
print('Cleaned Wavelengths:',len(wavelengths_clean))

Cleaned Data Dimensions: (1000, 1000, 173)
Cleaned Wavelengths: 173

start_time = time()
# run kmeans with 5 clusters and 50 iterations
(m,c) = kmeans(refl_clean, 5, 50) 
print("--- It took %s minutes to run kmeans on the reflectance data ---" % round((time() - start_time)/60,1))

spectral:INFO: k-means iteration 1 - 490247 pixels reassigned.
k-means iteration 1 - 490247 pixels reassigned.
spectral:INFO: k-means iteration 2 - 89955 pixels reassigned.
k-means iteration 2 - 89955 pixels reassigned.
spectral:INFO: k-means iteration 3 - 30478 pixels reassigned.
k-means iteration 3 - 30478 pixels reassigned.
spectral:INFO: k-means iteration 4 - 24830 pixels reassigned.
k-means iteration 4 - 24830 pixels reassigned.
spectral:INFO: k-means iteration 5 - 22757 pixels reassigned.
k-means iteration 5 - 22757 pixels reassigned.
spectral:INFO: k-means iteration 6 - 23624 pixels reassigned.
k-means iteration 6 - 23624 pixels reassigned.
spectral:INFO: k-means iteration 7 - 25620 pixels reassigned.
k-means iteration 7 - 25620 pixels reassigned.
spectral:INFO: k-means iteration 8 - 27850 pixels reassigned.
k-means iteration 8 - 27850 pixels reassigned.
spectral:INFO: k-means iteration 9 - 32993 pixels reassigned.
k-means iteration 9 - 32993 pixels reassigned.
spectral:INFO: k-means iteration 10 - 43142 pixels reassigned.
k-means iteration 10 - 43142 pixels reassigned.
spectral:INFO: k-means iteration 11 - 49025 pixels reassigned.
k-means iteration 11 - 49025 pixels reassigned.
spectral:INFO: k-means iteration 12 - 48883 pixels reassigned.
k-means iteration 12 - 48883 pixels reassigned.
spectral:INFO: k-means iteration 13 - 44926 pixels reassigned.
k-means iteration 13 - 44926 pixels reassigned.
spectral:INFO: k-means iteration 14 - 39465 pixels reassigned.
k-means iteration 14 - 39465 pixels reassigned.
spectral:INFO: k-means iteration 15 - 33889 pixels reassigned.
k-means iteration 15 - 33889 pixels reassigned.
spectral:INFO: k-means iteration 16 - 28935 pixels reassigned.
k-means iteration 16 - 28935 pixels reassigned.
spectral:INFO: k-means iteration 17 - 24496 pixels reassigned.
k-means iteration 17 - 24496 pixels reassigned.
spectral:INFO: k-means iteration 18 - 20838 pixels reassigned.
k-means iteration 18 - 20838 pixels reassigned.
spectral:INFO: k-means iteration 19 - 17702 pixels reassigned.
k-means iteration 19 - 17702 pixels reassigned.
spectral:INFO: k-means iteration 20 - 14771 pixels reassigned.
k-means iteration 20 - 14771 pixels reassigned.
spectral:INFO: k-means iteration 21 - 12259 pixels reassigned.
k-means iteration 21 - 12259 pixels reassigned.
spectral:INFO: k-means iteration 22 - 10129 pixels reassigned.
k-means iteration 22 - 10129 pixels reassigned.
spectral:INFO: k-means iteration 23 - 8446 pixels reassigned.
k-means iteration 23 - 8446 pixels reassigned.
spectral:INFO: k-means iteration 24 - 6892 pixels reassigned.
k-means iteration 24 - 6892 pixels reassigned.
spectral:INFO: k-means iteration 25 - 5754 pixels reassigned.
k-means iteration 25 - 5754 pixels reassigned.
spectral:INFO: k-means iteration 26 - 4726 pixels reassigned.
k-means iteration 26 - 4726 pixels reassigned.
spectral:INFO: k-means iteration 27 - 4046 pixels reassigned.
k-means iteration 27 - 4046 pixels reassigned.
spectral:INFO: k-means iteration 28 - 3425 pixels reassigned.
k-means iteration 28 - 3425 pixels reassigned.
spectral:INFO: k-means iteration 29 - 3045 pixels reassigned.
k-means iteration 29 - 3045 pixels reassigned.
spectral:INFO: k-means iteration 30 - 2681 pixels reassigned.
k-means iteration 30 - 2681 pixels reassigned.
spectral:INFO: k-means iteration 31 - 2336 pixels reassigned.
k-means iteration 31 - 2336 pixels reassigned.
spectral:INFO: k-means iteration 32 - 2446 pixels reassigned.
k-means iteration 32 - 2446 pixels reassigned.
spectral:INFO: k-means iteration 33 - 2572 pixels reassigned.
k-means iteration 33 - 2572 pixels reassigned.
spectral:INFO: k-means iteration 34 - 2790 pixels reassigned.
k-means iteration 34 - 2790 pixels reassigned.
spectral:INFO: k-means iteration 35 - 3218 pixels reassigned.
k-means iteration 35 - 3218 pixels reassigned.
spectral:INFO: k-means iteration 36 - 3594 pixels reassigned.
k-means iteration 36 - 3594 pixels reassigned.
spectral:INFO: k-means iteration 37 - 4278 pixels reassigned.
k-means iteration 37 - 4278 pixels reassigned.
spectral:INFO: k-means iteration 38 - 5336 pixels reassigned.
k-means iteration 38 - 5336 pixels reassigned.
spectral:INFO: k-means iteration 39 - 7232 pixels reassigned.
k-means iteration 39 - 7232 pixels reassigned.
spectral:INFO: k-means iteration 40 - 11120 pixels reassigned.
k-means iteration 40 - 11120 pixels reassigned.
spectral:INFO: k-means iteration 41 - 18893 pixels reassigned.
k-means iteration 41 - 18893 pixels reassigned.
spectral:INFO: k-means iteration 42 - 30213 pixels reassigned.
k-means iteration 42 - 30213 pixels reassigned.
spectral:INFO: k-means iteration 43 - 37760 pixels reassigned.
k-means iteration 43 - 37760 pixels reassigned.
spectral:INFO: k-means iteration 44 - 35916 pixels reassigned.
k-means iteration 44 - 35916 pixels reassigned.
spectral:INFO: k-means iteration 45 - 28996 pixels reassigned.
k-means iteration 45 - 28996 pixels reassigned.
spectral:INFO: k-means iteration 46 - 21972 pixels reassigned.
k-means iteration 46 - 21972 pixels reassigned.
spectral:INFO: k-means iteration 47 - 16906 pixels reassigned.
k-means iteration 47 - 16906 pixels reassigned.
spectral:INFO: k-means iteration 48 - 13454 pixels reassigned.
k-means iteration 48 - 13454 pixels reassigned.
spectral:INFO: k-means iteration 49 - 10685 pixels reassigned.
k-means iteration 49 - 10685 pixels reassigned.
spectral:INFO: k-means iteration 50 - 8758 pixels reassigned.
k-means iteration 50 - 8758 pixels reassigned.
spectral:INFO: kmeans terminated with 5 clusters after 50 iterations.
kmeans terminated with 5 clusters after 50 iterations.


--- It took 2.7 minutes to run kmeans on the reflectance data ---

Note that the algorithm still had on the order of 10000 clusters reassigning, when the 50 iterations were reached. You may extend the # of iterations.

Let's try that now:

start_time = time()
# run kmeans with 5 clusters and 40 iterations
(m, c) = kmeans(refl_clean, 5, 40, start_clusters=c) 
print("--- It took %s minutes to run kmeans on the reflectance data ---" % round((time() - start_time)/60,1))

spectral:INFO: k-means iteration 1 - 791183 pixels reassigned.
k-means iteration 1 - 791183 pixels reassigned.
spectral:INFO: k-means iteration 2 - 6618 pixels reassigned.
k-means iteration 2 - 6618 pixels reassigned.
spectral:INFO: k-means iteration 3 - 5791 pixels reassigned.
k-means iteration 3 - 5791 pixels reassigned.
spectral:INFO: k-means iteration 4 - 5019 pixels reassigned.
k-means iteration 4 - 5019 pixels reassigned.
spectral:INFO: k-means iteration 5 - 4565 pixels reassigned.
k-means iteration 5 - 4565 pixels reassigned.
spectral:INFO: k-means iteration 6 - 4180 pixels reassigned.
k-means iteration 6 - 4180 pixels reassigned.
spectral:INFO: k-means iteration 7 - 3956 pixels reassigned.
k-means iteration 7 - 3956 pixels reassigned.
spectral:INFO: k-means iteration 8 - 3625 pixels reassigned.
k-means iteration 8 - 3625 pixels reassigned.
spectral:INFO: k-means iteration 9 - 3271 pixels reassigned.
k-means iteration 9 - 3271 pixels reassigned.
spectral:INFO: k-means iteration 10 - 2978 pixels reassigned.
k-means iteration 10 - 2978 pixels reassigned.
spectral:INFO: k-means iteration 11 - 2612 pixels reassigned.
k-means iteration 11 - 2612 pixels reassigned.
spectral:INFO: k-means iteration 12 - 2348 pixels reassigned.
k-means iteration 12 - 2348 pixels reassigned.
spectral:INFO: k-means iteration 13 - 2138 pixels reassigned.
k-means iteration 13 - 2138 pixels reassigned.
spectral:INFO: k-means iteration 14 - 1863 pixels reassigned.
k-means iteration 14 - 1863 pixels reassigned.
spectral:INFO: k-means iteration 15 - 1728 pixels reassigned.
k-means iteration 15 - 1728 pixels reassigned.
spectral:INFO: k-means iteration 16 - 1539 pixels reassigned.
k-means iteration 16 - 1539 pixels reassigned.
spectral:INFO: k-means iteration 17 - 1305 pixels reassigned.
k-means iteration 17 - 1305 pixels reassigned.
spectral:INFO: k-means iteration 18 - 1206 pixels reassigned.
k-means iteration 18 - 1206 pixels reassigned.
spectral:INFO: k-means iteration 19 - 1060 pixels reassigned.
k-means iteration 19 - 1060 pixels reassigned.
spectral:INFO: k-means iteration 20 - 982 pixels reassigned.
k-means iteration 20 - 982 pixels reassigned.
spectral:INFO: k-means iteration 21 - 913 pixels reassigned.
k-means iteration 21 - 913 pixels reassigned.
spectral:INFO: k-means iteration 22 - 795 pixels reassigned.
k-means iteration 22 - 795 pixels reassigned.
spectral:INFO: k-means iteration 23 - 732 pixels reassigned.
k-means iteration 23 - 732 pixels reassigned.
spectral:INFO: k-means iteration 24 - 671 pixels reassigned.
k-means iteration 24 - 671 pixels reassigned.
spectral:INFO: k-means iteration 25 - 600 pixels reassigned.
k-means iteration 25 - 600 pixels reassigned.
spectral:INFO: k-means iteration 26 - 533 pixels reassigned.
k-means iteration 26 - 533 pixels reassigned.
spectral:INFO: k-means iteration 27 - 469 pixels reassigned.
k-means iteration 27 - 469 pixels reassigned.
spectral:INFO: k-means iteration 28 - 393 pixels reassigned.
k-means iteration 28 - 393 pixels reassigned.
spectral:INFO: k-means iteration 29 - 333 pixels reassigned.
k-means iteration 29 - 333 pixels reassigned.
spectral:INFO: k-means iteration 30 - 305 pixels reassigned.
k-means iteration 30 - 305 pixels reassigned.
spectral:INFO: k-means iteration 31 - 266 pixels reassigned.
k-means iteration 31 - 266 pixels reassigned.
spectral:INFO: k-means iteration 32 - 215 pixels reassigned.
k-means iteration 32 - 215 pixels reassigned.
spectral:INFO: k-means iteration 33 - 171 pixels reassigned.
k-means iteration 33 - 171 pixels reassigned.
spectral:INFO: k-means iteration 34 - 128 pixels reassigned.
k-means iteration 34 - 128 pixels reassigned.
spectral:INFO: k-means iteration 35 - 117 pixels reassigned.
k-means iteration 35 - 117 pixels reassigned.
spectral:INFO: k-means iteration 36 - 101 pixels reassigned.
k-means iteration 36 - 101 pixels reassigned.
spectral:INFO: k-means iteration 37 - 94 pixels reassigned.
k-means iteration 37 - 94 pixels reassigned.
spectral:INFO: k-means iteration 38 - 97 pixels reassigned.
k-means iteration 38 - 97 pixels reassigned.
spectral:INFO: k-means iteration 39 - 77 pixels reassigned.
k-means iteration 39 - 77 pixels reassigned.
spectral:INFO: k-means iteration 40 - 82 pixels reassigned.
k-means iteration 40 - 82 pixels reassigned.
spectral:INFO: kmeans terminated with 5 clusters after 40 iterations.
kmeans terminated with 5 clusters after 40 iterations.


--- It took 2.2 minutes to run kmeans on the reflectance data ---

Passing the initial clusters in sped up the convergence considerably, the second time around.

Let's take a look at the new cluster centers c. In this case, these represent spectral signatures of the five clusters (classes) that the data were grouped into. First we can take a look at the shape:

print(c.shape)

(5, 173)

c contains 5 groups of spectral curves with 173 bands (the # of bands we've kept after subsetting and removing the water vapor windows, first 10 noisy bands and last 15 noisy bands). We can plot these spectral classes as follows:

import pylab
pylab.figure()
for i in range(c.shape[0]):
    pylab.plot(wavelengths_clean, c[i],'.')
pylab.show
pylab.title('Spectral Classes from K-Means Clustering')
pylab.xlabel('Wavelength (nm)')
pylab.ylabel('Reflectance');

Next, we can look at the classes in map view, as well as a true color image.

view = imshow(refl_clean, bands=(58,34,19),stretch=0.01, classes=m, extent=refl_metadata['extent'])
view.set_display_mode('overlay')
view.class_alpha = 1 #set transparency
view.show_data;

view = imshow(refl_clean, bands=(24,12,4), stretch=0.03, extent=refl_metadata['extent'])
view.show_data;

Challenge Questions: K-Means

1. What do you think the spectral classes in the figure you just created represent?
2. Try using a different number of clusters in the kmeans algorithm (e.g., 3 or 10) to see what spectral classes and classifications result.
3. Try using different (higher) subset_factor in the clean_neon_refl_data function, like 3 or 5. Does this factor change the final classes that are created in the kmeans algorithm? By how much can you subset the data by and still achieve similar classification results?

Principal Component Analysis (PCA)

This next section follows the Spectral Python Dimensionality Reduction section closely.

Many of the bands within hyperspectral images are often strongly correlated. The principal components transformation represents a linear transformation of the original image bands to a set of new, uncorrelated features. These new features correspond to the eigenvectors of the image covariance matrix, where the associated eigenvalue represents the variance in the direction of the eigenvector. A very large percentage of the image variance can be captured in a relatively small number of principal components (compared to the original number of bands) .

pc = principal_components(refl_clean)
pc_view = imshow(pc.cov, extent=refl_metadata['extent'])
xdata = pc.transform(refl_clean)

In the covariance matrix display, lighter values indicate strong positive covariance, darker values indicate strong negative covariance, and grey values indicate covariance near zero.

To reduce dimensionality using principal components, we can sort the eigenvalues in descending order and then retain enough eigenvalues (and corresponding eigenvectors) to capture a desired fraction of the total image variance. We then reduce the dimensionality of the image pixels by projecting them onto the remaining eigenvectors. We will choose to retain a minimum of 99.9% of the total image variance.

pc_999 = pc.reduce(fraction=0.999)

# How many eigenvalues are left?
print('# of eigenvalues:',len(pc_999.eigenvalues))

img_pc = pc_999.transform(refl_clean)
print(img_pc.shape)

v = imshow(img_pc[:,:,:3], stretch_all=True, extent=refl_metadata['extent']);

# of eigenvalues: 10
(1000, 1000, 10)

You can see that even though we've only retained a subset of the bands, a lot of the details about the scene are still visible.

If you had training data, you could use a Gaussian maximum likelihood classifier (GMLC) for the reduced principal components to train and classify against the training data.

Challenge Question: PCA

Run the k-means classification after running PCA and see if you get similar results. Does reducing the data dimensionality affect the classification results?

In this tutorial, we will calculate the Normalized Difference Vegetation Index (NDVI) using Python functions.

This tutorial works with the Level 3 Spectrometer orthorectified surface directional reflectance - mosaic data product.

Calculate NDVI & Extract Spectra with Masks

Background:

The Normalized Difference Vegetation Index (NDVI) is a standard band-ratio calculation frequently used to analyze ecological remote sensing data. NDVI indicates whether the remotely-sensed target contains live green vegetation. When sunlight strikes objects, certain wavelengths of the electromagnetic spectrum are absorbed and other wavelengths are reflected. The pigment chlorophyll in plant leaves strongly absorbs visible light (with wavelengths in the range of 400-700 nm) for use in photosynthesis. The cell structure of the leaves, however, strongly reflects near-infrared light (wavelengths ranging from 700 - 1100 nm). Plants reflect up to 60% more light in the near infrared portion of the spectrum than they do in the green portion of the spectrum. By calculating the ratio of Near Infrared (NIR) to Visible (VIS) bands in hyperspectral data, we can obtain a metric of vegetation density and health.

The formula for NDVI is: $$NDVI = \frac{(NIR - VIS)}{(NIR+ VIS)}$$

Start by setting plot preferences and loading the neon_aop_hyperspectral.py module:

import dotenv
import os, sys
from copy import copy
import requests
import neonutilities as nu
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

This next function provides a handy way to download the Python module that we will use in this lesson.

# function to download data stored on the internet in a public url to a local file
def download_url(url,download_dir):
    if not os.path.isdir(download_dir):
        os.makedirs(download_dir)
    filename = url.split('/')[-1]
    r = requests.get(url, allow_redirects=True)
    file_object = open(os.path.join(download_dir,filename),'wb')
    file_object.write(r.content)

Download the module from its location on GitHub, add the python_modules to the path and import the neon_aop_hyperspectral.py module as neon_hs.

# download the neon_aop_hyperspectral.py module from GitHub
module_url = "https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/tutorials/Python/AOP/aop_python_modules/neon_aop_hyperspectral.py"
download_url(module_url,'../python_modules')

# add the python_modules to the path and import the python neon download and hyperspectral functions
sys.path.insert(0, '../python_modules')

# import the neon_aop_hyperspectral module
import neon_aop_hyperspectral as neon_hs;

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the dotenv package as follows, optionally specifying the path to the .env file. Adjust the savepath variable to point to your desired location; we recommend keeping this close to the root directory since the download path to the data file will be nested.

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

nu.by_tile_aop('DP3.30006.002',
               'SERC',
               2025,
               easting=368000,
               northing=4306000,
               token=token,
               include_provisional=True,
               savepath='C:/Data')

Provisional NEON data are included. To exclude provisional data, use input parameter include_provisional=False.


Continuing will download 2 NEON data files totaling approximately 661.3 MB. Do you want to proceed? (y/n)  y


Downloading 2 NEON data files totaling approximately 661.3 MB

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [01:20<00:00, 40.26s/it]

Click y when prompted to download the h5 data. Once the progress bar shows 100%, the reflectance data tile will be downloaded to the 'C:/NEON_Data/DP3.30006.002' directory. You can use the code cell below to walk through all the directories and display where the .h5 file was downloaded.

# display .h5 data in the savepath
for root, dirs, files in os.walk(r'C:\Data\DP3.30006.002'):
    for file in files:
        if file.endswith(".h5"):
             h5_tile = os.path.join(root, file)
             print(h5_tile)

C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

Read in SERC Reflectance Tile

# read the h5 reflectance file (including the full path) to the variable h5_file_name
print(f'h5_tile: {h5_tile}')

h5_tile: C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

serc_refl, serc_refl_md, wavelengths = neon_hs.aop_h5refl2array(h5_tile,'Reflectance')

Reading in  C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

Extract Visible and Near Infrared Bands

Now that we have uploaded all the required functions, we can calculate NDVI and plot it. Below we print the center wavelengths of the visible band (57) and near-infrared band (89):

print('band 58 center wavelength (nm): ', wavelengths[57])
print('band 90 center wavelength (nm) : ', wavelengths[89])

band 58 center wavelength (nm):  667.457214
band 90 center wavelength (nm) :  827.793396

Calculate NDVI and Plot NDVI Maps

Here we see that band 58 represents red visible light, while band 90 is in the NIR portion of the spectrum. Let's extract these two bands from the reflectance array and calculate the ratio using the numpy.true_divide which divides arrays element-wise. This also handles a case where the denominator = 0, which would otherwise throw a warning or error.

vis = serc_refl[:,:,57]
nir = serc_refl[:,:,89]

# handle a divide by zero by setting the numpy errstate as follows
with np.errstate(divide='ignore', invalid='ignore'):
    ndvi = np.true_divide((nir-vis),(nir+vis))
    ndvi[ndvi == np.inf] = 0
    ndvi = np.nan_to_num(ndvi)

Let's take a look at the min, mean, and max values of NDVI that we calculated:

print(f'NDVI Min: {round(ndvi.min(),2)}')
print(f'NDVI Mean: {round(ndvi.mean(),2)}')
print(f'NDVI Max: {ndvi.max()}')

NDVI Min: -0.93
NDVI Mean: 0.62
NDVI Max: 1.0

We can use the function plot_aop_refl to plot this, and choose the seismic color pallette to highlight the difference between positive and negative NDVI values. Since this is a normalized index, the values should range from -1 to +1.

neon_hs.plot_aop_refl(ndvi,serc_refl_md['extent'],
                      colorlimit = (np.min(ndvi),np.max(ndvi)),
                      title='SERC Subset NDVI \n (VIS = Band 58, NIR = Band 90)',
                      cmap_title='NDVI',
                      colormap='seismic')

You can see that the water bodies have negative NDVI values, roads and buildings have NDVI values around 0, and vegetation has NDVI > 0. On your own, try out different color maps to see more nuances within the positive NDVI values.

Extract Spectra Using Masks

In the second part of this tutorial, we will learn how to extract the average spectra of pixels whose NDVI exceeds a specified threshold value. There are several ways to do this using numpy, including the mask functions numpy.ma, as well as numpy.where and finally using boolean indexing.

To start, lets copy the NDVI calculated above and use booleans to create an array only containing NDVI > 0.6.

# make a copy of ndvi
ndvi_gtpt6 = ndvi.copy()
#set all pixels with NDVI < 0.6 to nan, keeping only values > 0.6
ndvi_gtpt6[ndvi<0.6] = np.nan  
print('Mean NDVI > 0.6:',round(np.nanmean(ndvi_gtpt6),2))

Mean NDVI > 0.6: 0.85

Now let's plot the values of NDVI after masking out values < 0.6.

neon_hs.plot_aop_refl(ndvi_gtpt6,
                      serc_refl_md['extent'],
                      colorlimit=(0.6,1),
                      title='SERC Subset NDVI > 0.6 \n (VIS = Band 58, NIR = Band 90)',
                      cmap_title='NDVI',
                      colormap='RdYlGn')

Calculate the mean spectra, thresholded by NDVI

Below we will demonstrate how to calculate statistics on arrays where you have applied a mask numpy.ma. In this example, the function calculates the mean spectra for values that remain after masking out values by a specified threshold.

import numpy.ma as ma
def calculate_mean_masked_spectra(refl_array,ndvi,ndvi_threshold,ineq='>'):
    mean_masked_refl = np.zeros(refl_array.shape[2])
    for i in np.arange(refl_array.shape[2]):
        refl_band = refl_array[:,:,i]
        if ineq == '>':
            ndvi_mask = ma.masked_where((ndvi<=ndvi_threshold) | (np.isnan(ndvi)),ndvi)
        elif ineq == '<':
            ndvi_mask = ma.masked_where((ndvi>=ndvi_threshold) | (np.isnan(ndvi)),ndvi)   
        else:
            print('ERROR: Invalid inequality. Enter < or >')
        masked_refl = ma.MaskedArray(refl_band,mask=ndvi_mask.mask)
        mean_masked_refl[i] = ma.mean(masked_refl)
    return mean_masked_refl

We can test out this function for various NDVI thresholds. We'll test two together, and you can try out different values on your own. Let's look at the average spectra for healthy vegetation (NDVI > 0.6), and for a lower threshold (NDVI < 0.3).

serc_ndvi_gtpt6 = calculate_mean_masked_spectra(serc_refl,ndvi,0.6)
serc_ndvi_ltpt3 = calculate_mean_masked_spectra(serc_refl,ndvi,0.3,ineq='<') 

Finally, we can create a pandas dataframe of the wavelengths to plot the mean spectra.

#Remove water vapor bad band windows & last 10 bands 
w = wavelengths.copy()
w[((w >= 1340) & (w <= 1445)) | ((w >= 1790) & (w <= 1955))]=np.nan
w[-10:]=np.nan;  

nan_ind = np.argwhere(np.isnan(w))

serc_ndvi_gtpt6[nan_ind] = np.nan
serc_ndvi_ltpt3[nan_ind] = np.nan

#Create dataframe with masked NDVI mean spectra, scale by the reflectance scale factor
serc_ndvi_df = pd.DataFrame()
serc_ndvi_df['wavelength'] = w
serc_ndvi_df['mean_refl_ndvi_gtpt6'] = serc_ndvi_gtpt6/serc_refl_md['scale_factor']
serc_ndvi_df['mean_refl_ndvi_ltpt3'] = serc_ndvi_ltpt3/serc_refl_md['scale_factor']

Let's take a look at the first 5 values of this new dataframe:

serc_ndvi_df.head()

Plot the masked NDVI dataframe to display the mean spectra for NDVI values that exceed 0.6 and that are less than 0.3:

ax = plt.gca();
serc_ndvi_df.plot(ax=ax,x='wavelength',y='mean_refl_ndvi_gtpt6',color='green',
                  edgecolor='none',kind='scatter',label='Mean Spectra where NDVI > 0.6',legend=True);
serc_ndvi_df.plot(ax=ax,x='wavelength',y='mean_refl_ndvi_ltpt3',color='red',
                  edgecolor='none',kind='scatter',label='Mean Spectra where NDVI < 0.3',legend=True);
ax.set_title('Mean Spectra of Reflectance Masked by NDVI')
ax.set_xlim([np.nanmin(w),np.nanmax(w)]);
ax.set_xlabel("Wavelength, nm"); ax.set_ylabel("Reflectance")
ax.grid('on'); 

This tutorial covers how to read in a NEON lidar Canopy Height Model (CHM) geotiff file into a Python rasterio object, shows some basic information about the raster data, and then ends with classifying the CHM into height bins.

In this tutorial, we will work with the NEON AOP L3 LiDAR ecoysystem structure (Canopy Height Model) data product. For more information about NEON data products and the CHM product DP3.30015.001, see the Ecosystem structure data product page on NEON's Data Portal.

First, let's import the required packages and set our plot display to be in-line:

import dotenv
import os
import copy
import neonutilities as nu
import numpy as np
import rasterio as rio
from rasterio.plot import show, show_hist
import matplotlib.pyplot as plt

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the python-dotenv package as follows, optionally specifying the path to the .env file in load_dotenv().

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

Next, let's download a single tile (1 km x 1 km CHM file) using nu.by_tile_aop().

nu.by_tile_aop(dpid='DP3.30015.001',
               site='TEAK',
               year='2024',
               easting=320000,
               northing=4092000,
               token=token,
               savepath=r'C:\NEON_Data') # change if desired

Provisional NEON data are not included. To download provisional data, use input parameter include_provisional=True.


Continuing will download 2 NEON data files totaling approximately 2.9 MB. Do you want to proceed? (y/n)  y


Downloading 2 NEON data files totaling approximately 2.9 MB

100%|███████████████████████████████████████| 2/2 [00:00<00:00,  2.22it/s]

# iterate over directory recursively to show path of downloaded CHM.tif file
for root, dirs, files in os.walk(r'C:\NEON_Data\DP3.30015.001'):
    for name in files:
        if name.endswith('.tif'):
            chm_tile = os.path.join(root, name)
            print(chm_tile) 

C:\NEON_Data\DP3.30015.001\neon-aop-products\2024\FullSite\D17\2024_TEAK_7\L3\DiscreteLidar\CanopyHeightModelGtif\NEON_D17_TEAK_DP3_320000_4092000_CHM.tif

Open a GeoTIFF with rasterio

Let's look at the TEAK Canopy Height Model (CHM) to start. We can open and read this in Python using the rasterio.open function:

# read the chm file to the variable chm_dataset
chm_dataset = rio.open(chm_tile)

Now we can look at a few properties of this dataset to start to get a feel for the rasterio object:

print('chm_dataset:\n',chm_dataset)
print('\nshape:\n',chm_dataset.shape)
print('\nno data value:\n',chm_dataset.nodata)
print('\nspatial extent:\n',chm_dataset.bounds)
print('\ncoordinate information (crs):\n',chm_dataset.crs)

chm_dataset:
 <open DatasetReader name='C:\NEON_Data\DP3.30015.001\neon-aop-products\2024\FullSite\D17\2024_TEAK_7\L3\DiscreteLidar\CanopyHeightModelGtif\NEON_D17_TEAK_DP3_320000_4092000_CHM.tif' mode='r'>

shape:
 (1000, 1000)

no data value:
 -9999.0

spatial extent:
 BoundingBox(left=320000.0, bottom=4092000.0, right=321000.0, top=4093000.0)

coordinate information (crs):
 PROJCS["WGS 84 / UTM zone 11N",GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

Plot the Canopy Height Map and Histogram

We can use rasterio's built-in functions show and show_hist to plot and visualize the CHM tile. It is often useful to plot a histogram of the geotiff data in order to get a sense of the range and distribution of values.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
show(chm_dataset, ax=ax1);

show_hist(chm_dataset, bins=50, histtype='stepfilled',
          lw=0.0, stacked=False, alpha=0.3, ax=ax2);
ax2.set_xlabel("Canopy Height (meters)");
ax2.get_legend().remove()

plt.show();

On your own, adjust the number of bins, and range of the y-axis to get a better sense of the distribution of the canopy height values. We can see that a large portion of the values are zero. These correspond to bare ground. Let's look at a histogram and plot the data without these zero values which are dominating the frequency distribution. To do this, we'll remove all values > 2 m. Due to the vertical range resolution of the lidar sensor, data collected with the older Optech Gemini sensor can only resolve the ground to within 2 m, so anything below that height would be rounded down to zero. Our newer sensors (Riegl Q780 and Optech Galaxy Prime) have a higher range resolution, so the ground can be resolved to within ~0.7 m. To see which lidar sensor collected a given site, refer to the table at the bottom of the Flight Schedules and Coverage page (https://www.neonscience.org/data-collection/flight-schedules-coverage).

chm_data = chm_dataset.read(1)
valid_data = chm_data[chm_data>2]
plt.hist(valid_data.flatten(),bins=30);

From the histogram we can see that the majority of the trees are < 60m. The frequency of tall trees rapidly drops off.

Threshold Based Raster Classification

Next, we will create a classified raster object. To do this, we will use the numpy.where function to create a new raster based off boolean classifications. Let's classify the canopy height into five groups:

• Class 1: CHM = 0 m
• Class 2: 0m < CHM <= 15m
• Class 3: 10m < CHM <= 30m
• Class 4: 20m < CHM <= 45m
• Class 5: CHM > 45m

We can use np.where to find the indices where the specified criteria is met.

chm_reclass = chm_data.copy()
chm_reclass[np.where(chm_data==0)] = 1 # CHM = 0 : Class 1
chm_reclass[np.where((chm_data>0) & (chm_data<=15))] = 2 # 0m < CHM <= 10m - Class 2
chm_reclass[np.where((chm_data>15) & (chm_data<=30))] = 3 # 10m < CHM <= 20m - Class 3
chm_reclass[np.where((chm_data>30) & (chm_data<=45))] = 4 # 20m < CHM <= 30m - Class 4
chm_reclass[np.where(chm_data>45)] = 5 # CHM > 30m - Class 5

When we look at this variable, we can see that it is now populated with values between 1-5:

chm_reclass

array([[1., 1., 1., ..., 3., 3., 3.],
       [2., 2., 2., ..., 3., 3., 3.],
       [1., 2., 2., ..., 3., 3., 3.],
       ...,
       [3., 1., 4., ..., 2., 2., 2.],
       [3., 1., 1., ..., 2., 2., 2.],
       [1., 1., 1., ..., 2., 2., 1.]], shape=(1000, 1000))

Lastly we can use matplotlib to display this re-classified CHM. We will define our own colormap to plot these discrete classifications, and create a custom legend to label the classes. First, to include the spatial information in the plot, create a new variable called ext that pulls from the rasterio "bounds" field to create the extent in the expected format for plotting.

ext = [chm_dataset.bounds.left,
       chm_dataset.bounds.right,
       chm_dataset.bounds.bottom,
       chm_dataset.bounds.top]
ext

[320000.0, 321000.0, 4092000.0, 4093000.0]

import matplotlib.colors as colors
plt.figure(); 
cmap_chm = colors.ListedColormap(['lightblue','yellow','orange','green','red'])
plt.imshow(chm_reclass,extent=ext,cmap=cmap_chm)
plt.title('TEAK CHM Classification')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees

# Create custom legend to label the four canopy height classes:
import matplotlib.patches as mpatches
class1 = mpatches.Patch(color='lightblue', label='0 m')
class2 = mpatches.Patch(color='yellow', label='0-15 m')
class3 = mpatches.Patch(color='orange', label='15-30 m')
class4 = mpatches.Patch(color='green', label='30-45 m')
class5 = mpatches.Patch(color='red', label='>30 m')

ax.legend(handles=[class1,class2,class3,class4,class5],
          handlelength=0.7,bbox_to_anchor=(1.05, 0.4),loc='lower left',borderaxespad=0.);

In this tutorial, we will learn how to extract and plot a spectral reflectance profile (or spectral signature) from a single pixel of a reflectance band in a NEON hyperspectral HDF5 file.

This tutorial works with the Level 3 Spectrometer orthorectified surface bidirectional reflectance - mosaic.

In this lesson, we will cover how to extract and plot a spectral profile from a single pixel of a reflectance band in a NEON hyperspectral hdf5 file. To do this, we will use the aop_h5refl2array function to read in and clean our h5 reflectance data, and Python pandas to create a dataframe for the reflectance and associated wavelength data. We will end with an option example showing how to interactively view spectra from any pixel in a reflectance h5 tile.

Spectral Signatures

A spectral signature is a plot of the amount of light energy reflected by an object throughout the range of wavelengths in the electromagnetic spectrum. The spectral signature of an object conveys useful information about its structural and chemical composition. We can use these signatures to identify and classify different objects from a spectral image.

For example, vegetation has a distinct spectral signature.

Vegetation has a unique spectral signature characterized by high reflectance in the near infrared wavelengths, and much lower reflectance in the green portion of the visible spectrum. For more details, refer to Vegetation Analysis: Using Vegetation Indices in ENVI. We can extract reflectance values in the NIR and visible spectrums from hyperspectral data in order to map vegetation on the earth's surface. You can also use spectral curves as a proxy for vegetation health. We will explore this concept more in the next lesson, where we will calculate vegetation indices.

Let's get started. First import the required packages.

import os
import dotenv
import neonutilities as nu
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
import sys

This next function provides a handy way to download the Python module that we will use in this lesson. This uses the requests package.

# function to download data stored on the internet in a public url to a local file
def download_url(url,download_dir):
    if not os.path.isdir(download_dir):
        os.makedirs(download_dir)
    filename = url.split('/')[-1]
    r = requests.get(url, allow_redirects=True)
    file_object = open(os.path.join(download_dir,filename),'wb')
    file_object.write(r.content)

Download the module from its location on GitHub, add the python_modules to the path and import the neon_aop_hyperspectral.py module.

module_url = "https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/tutorials/Python/AOP/aop_python_modules/neon_aop_hyperspectral.py"
download_url(module_url,'../python_modules')
# os.listdir('../python_modules') #optionally show the contents of this directory to confirm the file downloaded

sys.path.insert(0, '../python_modules')

# import the neon_aop_hyperspectral module
import neon_aop_hyperspectral as neon_hs;

Now that we've imported the required packages and the hyperspectral module, we can download a reflectance dataset using neonutilities and start to explore it. We'll download data from the NEON site Smithsonian Environmental Research Center (SERC). First, use nu.list_available_dates to find what years of data are available for the bidirectional reflectance dataset at SERC.

nu.list_available_dates('DP3.30006.002','SERC')

PROVISIONAL Available Dates: 2022-05, 2025-06

Here we can see the available dates for this dataset. Let's use data from 2025.

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the dotenv package as follows, optionally specifying the path to the .env file.

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

Before downloading, we can check the extents in UTM x,y coordinates for the site as follows:

serc_2025_refl_exts = nu.get_aop_tile_extents('DP3.30006.002','SERC',2025,token=token)

Easting Bounds: (358000, 370000)
Northing Bounds: (4298000, 4312000)

For this example, download the reflectance tile with southwest coordinates of 368000, 4306000 using nu.by_tile_aop. Click y to continue the download after verifying the size (around 660 MB).

nu.by_tile_aop('DP3.30006.002',
               'SERC',
               2025,
               easting=368000,
               northing=4306000,
               token=token,
               include_provisional=True,
               savepath='C:/Data')

Provisional NEON data are included. To exclude provisional data, use input parameter include_provisional=False.


Continuing will download 2 NEON data files totaling approximately 661.3 MB. Do you want to proceed? (y/n)  y

The reflectance data tile is now downloaded into the 'C:/NEON_Data/DP3.30006.002' directory. You can use the code cell below to walk through all the directories and display where the .h5 file was downloaded.

# display .h5 data in the savepath
for root, dirs, files in os.walk(r'C:\Data\DP3.30006.002'):
    for file in files:
        if file.endswith(".h5"):
             h5_tile = os.path.join(root, file)
             print(h5_tile)

C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

# read in the data using the neon_hs module
serc_refl, serc_refl_md, wavelengths = neon_hs.aop_h5refl2array(h5_tile,'Reflectance')

Reading in  C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

Optionally, you can view the data stored in the metadata dictionary, and print the minimum, maximum, and mean reflectance values in the tile. In order to ignore NaN values, use numpy.nanmin/nanmax/nanmean.

for item in sorted(serc_refl_md):
    print(item + ':',serc_refl_md[item])

print('\nSERC Tile Reflectance Stats:')
print('min:',np.nanmin(serc_refl))
print('max:',round(np.nanmax(serc_refl),2))
print('mean:',round(np.nanmean(serc_refl),2))

EPSG: 32618
bad_band_window1: [1340 1445]
bad_band_window2: [1790 1955]
ext_dict: {'xMin': 368000.0, 'xMax': 369000.0, 'yMin': 4306000.0, 'yMax': 4307000.0}
extent: (368000.0, 369000.0, 4306000.0, 4307000.0)
no_data_value: -9999.0
projection: b'+proj=UTM +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
res: {'pixelWidth': 1.0, 'pixelHeight': 1.0}
scale_factor: 10000.0
shape: (1000, 1000, 426)
source: C:\Data\DP3.30006.002\neon-aop-provisional-products\2025\FullSite\D02\2025_SERC_7\L3\Spectrometer\Reflectance\NEON_D02_SERC_DP3_368000_4306000_bidirectional_reflectance.h5

SERC Tile Reflectance Stats:
min: -100
max: 15599
mean: 1219.84

For reference, plot the red band of the tile, using splicing, and the plot_aop_refl function:

sercb56 = serc_refl[:,:,55]/serc_refl_md['scale_factor']

neon_hs.plot_aop_refl(sercb56,
                      serc_refl_md['extent'],
                      colorlimit=(0,0.3),
                      title='SERC Tile Band 56',
                      cmap_title='Reflectance',
                      colormap='gist_earth')

We can use pandas to create a dataframe containing the wavelength and reflectance values for a single pixel - in this example, we'll look at the center pixel of the tile (500,500). To extract all reflectance values from a single pixel, use splicing as we did before to select a single band, but now we need to specify (y,x) and select all bands (using :).

serc_pixel_df = pd.DataFrame()
serc_pixel_df['reflectance'] = serc_refl[500,500,:]/serc_refl_md['scale_factor']
serc_pixel_df['wavelengths'] = wavelengths

We can preview the first and last five values of the dataframe using head and tail:

print(serc_pixel_df.head(5))
print(serc_pixel_df.tail(5))

   reflectance  wavelengths
0       0.0126   381.858398
1       0.0255   386.868896
2       0.0260   391.879395
3       0.0232   396.889893
4       0.0216   401.900391
     reflectance  wavelengths
421       0.5233  2491.281494
422       0.2784  2496.291992
423       0.1952  2501.302490
424       0.7879  2506.312988
425       1.4948  2511.323486

We can now plot the spectra, stored in this dataframe structure. pandas has a built in plotting routine, which can be called by typing .plot at the end of the dataframe.

serc_pixel_df.plot(x='wavelengths',y='reflectance',kind='scatter',edgecolor='none')
plt.title('Spectral Signature for SERC Pixel (500,500)')
ax = plt.gca() 
ax.set_xlim([np.min(serc_pixel_df['wavelengths']),np.max(serc_pixel_df['wavelengths'])])
ax.set_ylim(0,0.6)
ax.set_xlabel("Wavelength, nm")
ax.set_ylabel("Reflectance")
ax.grid('on')

Water Vapor Band Windows

We can see from the spectral profile above that there are spikes in reflectance around ~1400nm and ~1800nm. These result from water vapor which absorbs light between wavelengths 1340-1445 nm and 1790-1955 nm. The atmospheric correction that converts radiance to reflectance subsequently results in a spike at these two bands. The wavelengths of these water vapor bands is stored in the reflectance attributes, which is saved in the reflectance metadata dictionary created with h5refl2array:

bbw1 = serc_refl_md['bad_band_window1']; 
bbw2 = serc_refl_md['bad_band_window2']; 
print('Bad Band Window 1:',bbw1)
print('Bad Band Window 2:',bbw2)

Bad Band Window 1: [1340 1445]
Bad Band Window 2: [1790 1955]

Below we repeat the plot we made above, but this time draw in the edges of the water vapor band windows that we need to remove.

serc_pixel_df.plot(x='wavelengths',y='reflectance',kind='scatter',edgecolor='none');
plt.title('Spectral Signature for SERC Pixel (500,500)')
ax1 = plt.gca(); ax1.grid('on')
ax1.set_xlim([np.min(serc_pixel_df['wavelengths']),np.max(serc_pixel_df['wavelengths'])]); 
ax1.set_ylim(0,0.5)
ax1.set_xlabel("Wavelength, nm"); ax1.set_ylabel("Reflectance")

#Add in red dotted lines to show boundaries of bad band windows:
ax1.plot((1340,1340),(0,1.5), 'r--');
ax1.plot((1445,1445),(0,1.5), 'r--');
ax1.plot((1790,1790),(0,1.5), 'r--');
ax1.plot((1955,1955),(0,1.5), 'r--');

We can now set these bad band windows to nan, along with the last 10 bands, which are also often noisy (as seen in the spectral profile plotted above). First make a copy of the wavelengths so that the original metadata doesn't change.

w = wavelengths.copy() #make a copy to deal with the mutable data type
w[((w >= 1340) & (w <= 1445)) | ((w >= 1790) & (w <= 1955))]=np.nan #can also use bbw1[0] or bbw1[1] to avoid hard-coding in
w[-10:]=np.nan;  # the last 10 bands sometimes have noise - best to eliminate
#print(w) #optionally print wavelength values to show that -9999 values are replaced with nan

Interactive Spectra Visualization

Finally, we can create a widget to interactively view the spectra of different pixels along the reflectance tile. Run the cell below, and select different pixel_x and pixel_y values to gain a sense of what the spectra look like for different materials on the ground.

from ipywidgets import interact

# --- Data Initialization Setup ---
# define refl_band, refl, and metadata, as copies of the original serc_refl data
refl_band = sercb56
refl = serc_refl.copy()
metadata = serc_refl_md.copy()

# pre-extract coordinate mapping info
ext = metadata['extent']

def interactive_spectra_plot(pixel_x, pixel_y):
    # extract reflectance vector for the targeted pixel
    reflectance = refl[pixel_y, pixel_x, :]
    
    # make pixel dataframe
    pixel_df = pd.DataFrame({
        'reflectance': reflectance,
        'wavelengths': w
    })

    # figure set-up
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))

    # --- Plot 1: Spectra Scatter ---
    pixel_df.plot(
        ax=ax1, 
        x='wavelengths', 
        y='reflectance', 
        kind='scatter', 
        edgecolor='none'
    )
    
    ax1.set_title(f"Spectra of Pixel ({pixel_x}, {pixel_y})")
    ax1.set_xlabel("Wavelength, nm")
    ax1.set_ylabel("Reflectance")
    ax1.set_xlim([np.min(wavelengths), np.max(wavelengths)])
    ax1.set_ylim([np.min(reflectance), np.max(reflectance) * 1.1])
    ax1.grid(True)  # FIX: Changed from 'on' to True for modern Python/Matplotlib

    # --- Plot 2: Pixel Location Map ---
    plot = ax2.imshow(refl_band, extent=ext, cmap='gist_earth', clim=(0, 0.1))
    ax2.set_title('Pixel Location')
    
    # colorbar
    cbar = fig.colorbar(plot, ax=ax2, aspect=20)
    cbar.set_label('Reflectance', rotation=90, labelpad=20)
    
    # tick adjustments
    ax2.ticklabel_format(useOffset=False, style='plain')
    ax2.tick_params(axis='x', labelrotation=90) # FIX: Avoided old plt.setp logic
    
    # calculate coordinate markers
    marker_x = ext[0] + pixel_x
    marker_y = ext[3] - pixel_y
    ax2.plot(marker_x, marker_y, marker='s', markersize=5, color='red')
    
    ax2.set_xlim(ext[0], ext[1])
    ax2.set_ylim(ext[2], ext[3])
    
    plt.show() 

# --- Run Interactive UI ---
interact(
    interactive_spectra_plot, 
    pixel_x=(0, refl.shape[1] - 1, 1),
    pixel_y=(0, refl.shape[0] - 1, 1)
);

In this tutorial, we will learn how to extract and plot a spectral profile from a single pixel of a reflectance band in a NEON hyperspectral HDF5 file.

This tutorial uses the mosaiced or tiled NEON data product. For a tutorial using the flightline data, please see Plot a Spectral Signature in Python - Flightline Data.

In this exercise, we will learn how to extract and plot a spectral profile from a single pixel of a reflectance band in a NEON hyperspectral hdf5 file. To do this, we will use the aop_h5refl2array function to read in and clean our h5 reflectance data, and the Python package pandas to create a dataframe for the reflectance and associated wavelength data.

Spectral Signatures

A spectral signature is a plot of the amount of light energy reflected by an object throughout the range of wavelengths in the electromagnetic spectrum. The spectral signature of an object conveys useful information about its structural and chemical composition. We can use these signatures to identify and classify different objects from a spectral image.

Vegetation has a unique spectral signature characterized by high reflectance in the near infrared wavelengths, and much lower reflectance in the green portion of the visible spectrum. We can extract reflectance values in the NIR and visible spectrums from hyperspectral data in order to map vegetation on the earth's surface. You can also use spectral curves as a proxy for vegetation health. We will explore this concept more in the next lesson, where we will caluclate vegetation indices.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
import warnings
warnings.filterwarnings('ignore') #don't display warnings

Import the hyperspectral functions file that you downloaded into the variable neon_hs (for neon hyperspectral):

import os

# Note: you will need to update this filepath according to your local machine
os.chdir("/Users/olearyd/Git/data/")
import neon_aop_hyperspectral as neon_hs

# Note: you will need to update this filepath according to your local machine
sercRefl, sercRefl_md = neon_hs.aop_h5refl2array('/Users/olearyd/Git/data/NEON_D02_SERC_DP3_368000_4306000_reflectance.h5')

Optionally, you can view the data stored in the metadata dictionary, and print the minimum, maximum, and mean reflectance values in the tile. In order to handle any nan values, use Numpy nanmin nanmax and nanmean.

for item in sorted(sercRefl_md):
    print(item + ':',sercRefl_md[item])

print('SERC Tile Reflectance Stats:')
print('min:',np.nanmin(sercRefl))
print('max:',round(np.nanmax(sercRefl),2))
print('mean:',round(np.nanmean(sercRefl),2))

For reference, plot the red band of the tile, using splicing, and the plot_aop_refl function:

sercb56 = sercRefl[:,:,55]

neon_hs.plot_aop_refl(sercb56,
                      sercRefl_md['spatial extent'],
                      colorlimit=(0,0.3),
                      title='SERC Tile Band 56',
                      cmap_title='Reflectance',
                      colormap='gist_earth')

We can use pandas to create a dataframe containing the wavelength and reflectance values for a single pixel - in this example, we'll look at the center pixel of the tile (500,500).

import pandas as pd

To extract all reflectance values from a single pixel, use splicing as we did before to select a single band, but now we need to specify (y,x) and select all bands (using :).

serc_pixel_df = pd.DataFrame()
serc_pixel_df['reflectance'] = sercRefl[500,500,:]
serc_pixel_df['wavelengths'] = sercRefl_md['wavelength']

We can preview the first and last five values of the dataframe using head and tail:

print(serc_pixel_df.head(5))
print(serc_pixel_df.tail(5))

   reflectance  wavelengths
0       0.0860   383.534302
1       0.0667   388.542206
2       0.0531   393.550110
3       0.0434   398.558014
4       0.0375   403.565887
     reflectance  wavelengths
421       0.7394  2491.863037
422       0.2232  2496.870850
423       0.5458  2501.878906
424       1.4881  2506.886719
425       1.4882  2511.894531

We can now plot the spectra, stored in this dataframe structure. pandas has a built in plotting routine, which can be called by typing .plot at the end of the dataframe.

serc_pixel_df.plot(x='wavelengths',y='reflectance',kind='scatter',edgecolor='none')
plt.title('Spectral Signature for SERC Pixel (500,500)')
ax = plt.gca() 
ax.set_xlim([np.min(serc_pixel_df['wavelengths']),np.max(serc_pixel_df['wavelengths'])])
ax.set_ylim([np.min(serc_pixel_df['reflectance']),np.max(serc_pixel_df['reflectance'])])
ax.set_xlabel("Wavelength, nm")
ax.set_ylabel("Reflectance")
ax.grid('on')

Water Vapor Band Windows

We can see from the spectral profile above that there are spikes in reflectance around ~1400nm and ~1800nm. These result from water vapor which absorbs light between wavelengths 1340-1445 nm and 1790-1955 nm. The atmospheric correction that converts radiance to reflectance subsequently results in a spike at these two bands. The wavelengths of these water vapor bands is stored in the reflectance attributes, which is saved in the reflectance metadata dictionary created with h5refl2array:

bbw1 = sercRefl_md['bad band window1']; 
bbw2 = sercRefl_md['bad band window2']; 
print('Bad Band Window 1:',bbw1)
print('Bad Band Window 2:',bbw2)

Bad Band Window 1: [1340 1445]
Bad Band Window 2: [1790 1955]

Below we repeat the plot we made above, but this time draw in the edges of the water vapor band windows that we need to remove.

serc_pixel_df.plot(x='wavelengths',y='reflectance',kind='scatter',edgecolor='none');
plt.title('Spectral Signature for SERC Pixel (500,500)')
ax1 = plt.gca(); ax1.grid('on')
ax1.set_xlim([np.min(serc_pixel_df['wavelengths']),np.max(serc_pixel_df['wavelengths'])]); 
ax1.set_ylim(0,0.5)
ax1.set_xlabel("Wavelength, nm"); ax1.set_ylabel("Reflectance")

#Add in red dotted lines to show boundaries of bad band windows:
ax1.plot((1340,1340),(0,1.5), 'r--')
ax1.plot((1445,1445),(0,1.5), 'r--')
ax1.plot((1790,1790),(0,1.5), 'r--')
ax1.plot((1955,1955),(0,1.5), 'r--')

[<matplotlib.lines.Line2D at 0x81aaccb70>]

We can now set these bad band windows to nan, along with the last 10 bands, which are also often noisy (as seen in the spectral profile plotted above). First make a copy of the wavelengths so that the original metadata doesn't change.

import copy
w = copy.copy(sercRefl_md['wavelength']) #make a copy to deal with the mutable data type
w[((w >= 1340) & (w <= 1445)) | ((w >= 1790) & (w <= 1955))]=np.nan #can also use bbw1[0] or bbw1[1] to avoid hard-coding in
w[-10:]=np.nan;  # the last 10 bands sometimes have noise - best to eliminate
#print(w) #optionally print wavelength values to show that -9999 values are replaced with nan

Interactive Spectra Visualization

Finally, we can create a widget to interactively view the spectra of different pixels along the reflectance tile. Run the two cells below, and interact with them to gain a better sense of what the spectra look like for different materials on the ground.

#define index corresponding to nan values:
nan_ind = np.argwhere(np.isnan(w))

#define refl_band, refl, and metadata 
refl_band = sercb56
refl = copy.copy(sercRefl)
metadata = copy.copy(sercRefl_md)

from IPython.html.widgets import *

def spectraPlot(pixel_x,pixel_y):

    reflectance = refl[pixel_y,pixel_x,:]
    reflectance[nan_ind]=np.nan
    
    pixel_df = pd.DataFrame()
    pixel_df['reflectance'] = reflectance
    pixel_df['wavelengths'] = w

    fig = plt.figure(figsize=(15,5))
    ax1 = fig.add_subplot(1,2,1)

    # fig, axes = plt.subplots(nrows=1, ncols=2)
    pixel_df.plot(ax=ax1,x='wavelengths',y='reflectance',kind='scatter',edgecolor='none');
    ax1.set_title('Spectra of Pixel (' + str(pixel_x) + ',' + str(pixel_y) + ')')
    ax1.set_xlim([np.min(metadata['wavelength']),np.max(metadata['wavelength'])]); 
    ax1.set_ylim([np.min(pixel_df['reflectance']),np.max(pixel_df['reflectance']*1.1)])
    ax1.set_xlabel("Wavelength, nm"); ax1.set_ylabel("Reflectance")
    ax1.grid('on')

    ax2 = fig.add_subplot(1,2,2)
    plot = plt.imshow(refl_band,extent=metadata['spatial extent'],clim=(0,0.1)); 
    plt.title('Pixel Location'); 
    cbar = plt.colorbar(plot,aspect=20); plt.set_cmap('gist_earth'); 
    cbar.set_label('Reflectance',rotation=90,labelpad=20); 
    ax2.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
    rotatexlabels = plt.setp(ax2.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
    
    ax2.plot(metadata['spatial extent'][0]+pixel_x,metadata['spatial extent'][3]-pixel_y,'s',markersize=5,color='red')
    ax2.set_xlim(metadata['spatial extent'][0],metadata['spatial extent'][1])
    ax2.set_ylim(metadata['spatial extent'][2],metadata['spatial extent'][3])
    
interact(spectraPlot, pixel_x = (0,refl.shape[1]-1,1),pixel_y=(0,refl.shape[0]-1,1))

This tutorial introduces NEON's Level 3 (mosaicked) RGB camera images, Data Product (DP3.30010.001) and uses the Python package rasterio to read in and plot the camera data in Python. In this lesson, we will read in an RGB camera tile collected over the NEON Smithsonian Environmental Research Center (SERC) site and plot the mutliband image, as well as the individual bands. This lesson was adapted from the rasterio plotting documentation.

Background

As part of the NEON Airborne Operation Platform's suite of remote sensing instruments, the digital camera produces high-resolution (<= 10 cm) photographs of the earth’s surface. The camera records light energy that has reflected off the ground in the visible portion (red, green and blue) of the electromagnetic spectrum. Often the camera images are used to provide context for the hyperspectral and LiDAR data, but they can also be used for research purposes in their own right. One such example is the tree-crown mapping work by Weinstein et al. - see the links below for more information!

• Individual Tree-Crown Detection in RGB Imagery Using Semi-Supervised Deep Learning Neural Networks
• A remote sensing derived data set of 100 million individual tree crowns for the National Ecological Observatory Network
• DeepForest: A Python package for RGB deep learning tree crown delineation

For more interactive notebooks showing examples of working with airborne camera imagery, including with the DeepForest package (and other environmental applications), check out:

• Environmental Data Science (EDS) Book and the EDS Notebook Gallery

In this lesson we will keep it simple and show how to read in and plot a single camera file (1km x 1km ortho-mosaicked tile) - a first step in any research incorporating the AOP camera data (in Python).

Tip: To run a code chunk (cell) in Jupyter Notebook you can either select Cell > Run Cells with your cursor placed in the cell you want to run, or use the shortcut key Shift + Enter. For more handy shortcuts, refer to the tab Help > Keyboard Shortcuts.

Import required packages

First let's import the packages that we'll be using in this lesson.

import os
import dotenv
import neonutilities as nu
import rasterio as rio
from rasterio.plot import show, show_hist
import matplotlib.pyplot as plt

Next, let's download a single camera file (1 km x 1 km tile).

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the python-dotenv package as follows, optionally specifying the path to the .env file in load_dotenv().

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

# download the RGB Camera data to the C:/data directory - change this if desired
nu.by_tile_aop(dpid='DP3.30010.001',
               site='SERC',
               year=2021,
               easting=368000,
               northing=4306000,
               token=token,
               savepath=r'C:\data')

Provisional NEON data are not included. To download provisional data, use input parameter include_provisional=True.


Continuing will download 2 NEON data files totaling approximately 68.8 MB. Do you want to proceed? (y/n)  y


Downloading 2 NEON data files totaling approximately 68.8 MB

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:01<00:00,  1.05it/s]

Display the RGB tile that you've downloaded:

rgb_dir = os.path.expanduser(r"C:\data\DP3.30010.001")

for root, dirs, files in os.walk(rgb_dir):
    for file in files:
        if file.endswith('.tif'):
            rgb_file = os.path.join(root, file)
            print(rgb_file)

C:\data\DP3.30010.001\neon-aop-products\2021\FullSite\D02\2021_SERC_5\L3\Camera\Mosaic\2021_SERC_5_368000_4306000_image.tif

Open the Camera RGB data with rasterio

We can open and read this RGB data that we downloaded in Python using the rasterio.open function:

# read the RGB file (including the full path) to the variable rgb_dataset
rgb_dataset = rio.open(rgb_file)

Let's look at a few properties of this dataset to get a sense of the information stored in the rasterio object:

print('rgb_dataset:\n',rgb_dataset)
print('\nshape:\n',rgb_dataset.shape)
print('\nspatial extent:\n',rgb_dataset.bounds)
print('\ncoordinate information (crs):\n',rgb_dataset.crs)

rgb_dataset:
 <open DatasetReader name='C:\data\DP3.30010.001\neon-aop-products\2021\FullSite\D02\2021_SERC_5\L3\Camera\Mosaic\2021_SERC_5_368000_4306000_image.tif' mode='r'>

shape:
 (10000, 10000)

spatial extent:
 BoundingBox(left=368000.0, bottom=4306000.0, right=369000.0, top=4307000.0)

coordinate information (crs):
 PROJCS["WGS 84 / UTM zone 18N",GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

Unlike the other AOP data products, camera imagery is generated at 10cm resolution, so each 1km x 1km tile will contain 10000 pixels (other 1m resolution data products will have 1000 x 1000 pixels per tile, where each pixel represents 1 meter).

Plot the RGB multiband image

We can use rasterio's built-in functions show to plot the CHM tile.

show(rgb_dataset);

Plot each band of the RGB image

We can also plot each band (red, green, and blue) individually as follows:

fig, (axr, axg, axb) = plt.subplots(1,3, figsize=(21,7))
show((rgb_dataset, 1), ax=axr, cmap='Reds', title='red channel')
show((rgb_dataset, 2), ax=axg, cmap='Greens', title='green channel')
show((rgb_dataset, 3), ax=axb, cmap='Blues', title='blue channel')
plt.show()

That's all for this example! Most of the other AOP raster data are all single band images so you can't make a 3-band composite like for the camera. You can make RGB composites using different bands of the hyperspectral data. In summary, rasterio is a handy Python package for working with any geotiff files. You can download and visualize the lidar and spectrometer derived raster images similarly.

There are myriad resources out there to learn programming in R. After linking to a tutorial on how to install R and RStudio on your computer, we then outline a few different paths to learn R basics depending on how you enjoy learning, and finally we include a few resources for intermediate and advanced learning.

Setting Up your Computer

Start out by installing R and, we recommend, RStudio, on your computer. RStudio is an Interactive Development Environment (IDE) for the R program. It is optional, but recommended when working with R. Directions for installing can be found within the tutorial
Install Git, Bash Shell, R & RStudio. You will need administrator permissions on your computer.

Pathways to Learning the Basics of R

In-person trainings

If you prefer to learn through in-person trainings, consider local workshops from The Carpentries Software Carpentry or Data Carpentry (generally ~$25 for a 2-day workshop), courses offered by a local college or university (prices vary), or organize your colleagues to meet regularly to learn R together (free!).

Online interactive courses

If you prefer to learn in a semi-structured online environment, there are a wide variety of online courses for learning R including Data Camp, Coursera, edX, and Lynda.com. Many of these options include free introductory lessons or trial periods as well as paid courses. We do not have personal experience with these courses and do not recommend or specifically promote any course.

In program interactive course

Swirl is guided introduction to R where you code along with the instructions in R. You get direct feedback when you type a command incorrectly. To use this package, once you have R or RStudio open and running, use the following commands to start the first lesson.

install.packages("swirl")

library(swirl)

swirl()

Online tutorials

If you prefer a less structured online environment, these tutorial series may be better suited for you.

• Software Carpentry’s Programming with R
  • Learn R with a focus on tools needed for effective programming. Beyond the basics, it covers functions, loops, command line, and other key skills
• Data Carpentry’s R for data analysis and visualization of Ecological Data
  • Learn R with a focus on data analysis. Beyond the basics, it covers dyplr for data aggregation & manipulation, ggplot2 for plotting, and touches on interacting with an SQL database. Designed to be taught by an instructor but the materials also work for independent learning online.
• Ethan White’s Data Carpentry for Biologists Semester Course (online content)
  • This comprehensive course contains an R section. While the overall focus is on data science skills, learning R is a portion of it (note, this is an extensive course).
• RStudio’s list
  • RStudio links to many other learning opportunities. Start with the 'Beginners' learning path.

Video tutorials

A blend of having an instructor and self-paced, video tutorials may also be of interest. New stand-alone video tutorials are out each day, so we aren’t going to recommend a specific series. Find what works for you by searching “R Programming video tutorials” on YouTube.

Books

Books are still a great way to learn R (and other languages). Many books are available at local libraries (university or community) or online, if you want to try them out before buying. Below are a few of the many, many books that data scientists working on the NEON project have found useful.

• Michael Crawley’s The R Book is a classic that takes you from beginning steps to analyses and modelling.
• Grolemun and Wickham’s R for Data Science focuses on using R in data science applications using Hadley Wickham’s “tidyverse”. It does assume some basic familiarity with R. Bonus: it is available online or in book format! (If you are completely new, they recommend starting with Hands-on Programming with R).

Beyond the Basics

There are many intermediate and advanced courses, lessons, and tutorials linked in the above resources. For example, the Swirl package offers intermediate and advanced courses on specific topics, as does RStudio's list. See courses here; development is ongoing so new courses may be added.

However, once the basics are handled, you will find that much of your learning will happen through solving individual problems you encounter. To solve these problems, your favorite search engine is your friend. Paste the error (without specifics to your file/data) into the search menu and find answers from those who have had similar questions.

For more on working with NEON data in particular, be sure to check out the other NEON data tutorials.

This tutorial provides the basics on how to set up Docker on one's local computer and then connect to an eddy4R Docker container in order to use the eddy4R R package. There are no specific skills needed for this tutorial, however, you will need to know how to access the command line tool for your operating system (basic instructions given).

In the tutorial below, we give the very barest of information to get Docker set up for use with the NEON R package eddy4R. For more information on using Docker, consider reading through the content from CyVerse's Container Camp's Intro to Docker.

Install Docker

To work with the eddy4R–Docker image, you first need to sign up for an account at DockerHub.

Once logged in, getting Docker up and running on your favorite operating system (Mac/Windows/Linux) is very easy. The "getting started" guide on Docker has detailed instructions for setting up Docker. Unless you plan on being a very active user and devoloper in Docker, we recommend starting with the stable channel (not edge channel) as you may encounter fewer problems.

• Mac
• Windows
• Linux

If you're using Docker for Windows make sure you have shared your drive.

If you're using an older version of Windows or MacOS, you may need to use Docker Machine instead.

Test Docker installation

Once you are done installing Docker, test your Docker installation by running the following command to make sure you are using version 1.13 or higher.

You will need an open shell window (Linux; Mac=Terminal) or the Docker Quickstart Terminal (Windows).

docker --version
	

When run, you will see which version of Docker you are currently running.

Note: If you run just the word docker you should see a whole bunch of lines showing the different options available with docker. Alternatively you can test your installation by running the following:

docker run hello-world

Notice that the first line states that the image can't be found locally. The next few lines are pulling the image, so if you were to run the hello-world prompt again, it would already be local and you'd see the message start at "Hello from Docker!".

If these steps work, you are ready to go on to access the eddy4R-Docker image that houses the suite of eddy4R R packages. If these steps have not worked, follow the installation instructions a second time.

Accessing eddy4R

Download of the eddy4R–Docker image and subsequent creation of a local container can be performed by two simple commands in an open shell (Linux; Mac = Terminal) or the Docker Quickstart Terminal (Windows).

The first command docker login will prompt you for your DockerHub ID and password.

The second command docker run -d -p 8787:8787 -e PASSWORD=YOURPASSWORD stefanmet/eddy4r:1.0.0 will download the latest eddy4R–Docker image and start a Docker container that utilizes port 8787 for establishing a graphical interface via web browser.

• docker run: docker will preform some process on an isolated container
• -d: the container will start in a detached mode, which means the container run in the background and will print the container ID
• -p: publish a container to a specified port (which follows)
• 8787:8787: specify which port you want to use. The default 8787:8787 is great if you are running locally. The first 4 digits are the port on your machine, the last 4 digits are the port communicating with RStudio on Docker. You can change the first 4 digits if you want to use a different port on your machine, or if you are running many containers or are on a shared network, but the last 4 digits need to be 8787.
• -e PASSWORD=YOURPASSWORD: define a password environmental variable to use upon login to the Rstudio instance. YOURPASSWORD can be anything you want.
• stefanmet/eddy4r:1.0.0: finally, which container do you want to run.

Now try it.

docker login 

docker run -d -p 8787:8787 -e PASSWORD=YOURPASSWORD stefanmet/eddy4r:1.0.0

This last command will run a specified release version (eddy4r:1.0.0) of the Docker image. Alternatively you can use eddy4r:latest to get the most up-to-date development image of eddy4r.

If you are using data stored on your local machine, rather than cloud hosting, a physical file system location on the host computer (local/dir) can be mounted to a file system location inside the Docker container (docker/dir). This is achieved with the Docker run option -v local/dir:docker/dir.

Access RStudio session

Now you can access the interactive RStudio session for using eddy4r by using any web browser and going to http://host-ip-address:8787 where host-ip-address is the internal IP address of the Docker host. For example, if your host IP address is 10.100.90.169 then you should type http://10.100.90.169:8787 into your browser.

To determine the IP address of your Docker host, follow the instructions below for your operating system.

Windows

Depending on the version of Docker, older Docker Toolbox versus the newer Docker Desktop for Windows, there are different way to get the docker machine IP address:

• Docker Toolbox - Type docker-machine ip default into cmd.exe window. The output will be your local IP address for the docker machine.
• Docker Desktop for Windows - Type ipconfig into cmd.exe window. The output will include either DockerNAT IPv4 address or vEthernet IPv4 address that docker uses to communicate to the internet, which in most cases will be 10.0.75.1.

Mac

Type ifconfig | grep "inet " | grep -v 127.0.0.1 into your Terminal window. The output will be one or more local IP addresses for the docker machine. Use the numbers after the first inet output.

Linux

Type localhost in a shell session and the local IP will be the output.

Once in the web browser you can log into this instance of the RStudio session with the username as rstudio and password as defined by YOURPASSWORD. Once complete you are now in a RStudio user interface with eddy4R installed and ready to use.

Additional information about the use of RStudio and eddy4R packages in Docker containers can be found on the rocker-org/rocker website and the eddy4RWiki pages.

Using eddy4R

To learn to use the eddy4R package to calculate fluxes, please visit the eddy4R vignette (link pending).