In this tutorial, we will learn how to plot spectral signatures of several different land cover types using an interactive click feature of the terra package.

Getting Started

First, we need to load our required packages and set the working directory.

# load required packages
library(rhdf5)
library(reshape2)
library(terra)
library(plyr)
library(ggplot2)
library(grDevices)

# set working directory, you can change the data path if desired

wd <- "~/data/" 
setwd(wd)

Download the reflectance tile, if you haven't already, using neonUtilities::byTileAOP:

byTileAOP(dpID = 'DP3.30006.001',
          site = 'SJER',
          year = '2021',
          easting = 257500,
          northing = 4112500,
          savepath = wd)

And then we can read in the hyperspectral hdf5 data. We will also collect a few other important pieces of information (band wavelengths and scaling factor) while we're at it.

# define filepath to the hyperspectral dataset
h5_file <- paste0(wd,"DP3.30006.001/neon-aop-products/2021/FullSite/D17/2021_SJER_5/L3/Spectrometer/Reflectance/NEON_D17_SJER_DP3_257000_4112000_reflectance.h5")

# read in the wavelength information from the HDF5 file
wavelengths <- h5read(h5_file,"/SJER/Reflectance/Metadata/Spectral_Data/Wavelength")

# grab scale factor from the Reflectance attributes
reflInfo <- h5readAttributes(h5_file,"/SJER/Reflectance/Reflectance_Data" )

scaleFact <- reflInfo$Scale_Factor

Now, we will read in the RGB image that we created in an earlier tutorial and plot it.

# read in RGB image as a 'stack' rather than a plain 'raster'
rgbStack <- rast(paste0(wd,"NEON_D17_SJER_DP3_257000_4112000_reflectance_2021_RGB.tif"))

# plot as RGB image, with a linear stretch
plotRGB(rgbStack,
        r=1,g=2,b=3, scale=300, 
        stretch = "lin")

Interactive click Function from the terra Package

Next, we use an interactive clicking function to identify the pixels that we want to extract spectral signatures for. To follow along with this tutorial, we suggest the following six cover types (exact locations shown in the image below).

1. Water
2. Tree canopy (avoid the shaded northwestern side of the tree)
3. Irrigated grass
4. Bare soil (baseball diamond infield)
5. Building roof (blue)
6. Road

As shown here:

For this next part, if you are following along in RStudio, you will need to enter these line below directly in the Console. dev.new(noRStudioGD = TRUE) will open up a separate window for plotting, which is where you will click the pixels to extract spectra, using the terra::click functionality.

dev.new(noRStudioGD = TRUE)

Now we can create our RGB plot, and start clicking on this in the pop-out Graphics window.

# change plotting parameters to better see the points and numbers generated from clicking
par(col="red", cex=2)

# use a histogram stretch in order to provide more contrast for selecting pixels
plotRGB(rgbStack, r=1, g=2, b=3, scale=300, stretch = "hist") 

# use the 'click' function
c <- click(rgbStack, n = 6, id=TRUE, xy=TRUE, cell=TRUE, type="p", pch=16, col="red", col.lab="red")

Once you have clicked your six points, the graphics window should close. If you want to choose new points, or if you accidentally clicked a point that you didn't intend to, run the previous 2 chunks of code again to re-start.

The click() function identifies the cell number that you clicked, but in order to extract spectral signatures, we need to convert that cell number into a row and column, as shown here:

# convert raster cell number into row and column (used to extract spectral signature below)
c$row <- c$cell%/%nrow(rgbStack)+1 # add 1 because R is 1-indexed
c$col <- c$cell%%ncol(rgbStack)

Extract Spectral Signatures from HDF5 file

Next, we will loop through each of the cells that and use the h5read() function to extract the reflectance values of all bands from the given pixel (row and column).

# create a new dataframe from the band wavelengths so that we can add the reflectance values for each cover type
pixel_df <- as.data.frame(wavelengths)

# loop through each of the cells that we selected

for(i in 1:length(c$cell)){
# extract spectral values from a single pixel
aPixel <- h5read(h5_file,"/SJER/Reflectance/Reflectance_Data",
                 index=list(NULL,c$col[i],c$row[i]))

# scale reflectance values from 0-1
aPixel <- aPixel/as.vector(scaleFact)

# reshape the data and turn into dataframe
b <- adply(aPixel,c(1))

# rename the column that we just created
names(b)[2] <- paste0("Point_",i)

# add reflectance values for this pixel to our combined data.frame called pixel_df
pixel_df <- cbind(pixel_df,b[2])
}

Plot Spectral signatures using ggplot2

Finally, we have everything that we need to plot the spectral signatures for each of the pixels that we clicked. In order to color our lines by the different land cover types, we will first reshape our data using the melt() function, then plot the spectral signatures.

# Use the melt() function to reshape the dataframe into a format that ggplot prefers
pixel.melt <- reshape2::melt(pixel_df, id.vars = "wavelengths", value.name = "Reflectance")

# Now, let's plot the spectral signatures!

ggplot()+
  geom_line(data = pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
                      labels = c("Water","Tree","Grass","Soil","Roof","Road"))+
  labs(color = "Cover Type")+
  ggtitle("Land cover spectral signatures")+
  theme(plot.title = element_text(hjust = 0.5, size=20))+
  xlab("Wavelength")

Nice! However, there seems to be something weird going on in the wavelengths near ~1400nm and ~1850 nm...

Atmospheric Absorption Bands

Those irregularities around 1400nm and 1850 nm are two major atmospheric absorption bands - regions where gasses in the atmosphere (primarily carbon dioxide and water vapor) absorb radiation, and therefore, obscure the reflected radiation that the imaging spectrometer measures. Fortunately, the lower and upper bound of each of those atmospheric absorption bands is specified in the HDF5 file. Let's read those bands and plot rectangles where the reflectance measurements are obscured by atmospheric absorption.

# grab reflectance metadata (which contains absorption band limits)
reflMetadata <- h5readAttributes(h5_file,"/SJER/Reflectance" )
ab1 <- reflMetadata$Band_Window_1_Nanometers
ab2 <- reflMetadata$Band_Window_2_Nanometers

# Plot spectral signatures again with grey rectangles highlighting the absorption bands
ggplot()+
  geom_line(data = pixel.melt, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
  geom_rect(mapping=aes(ymin=min(pixel.melt$Reflectance),ymax=max(pixel.melt$Reflectance), xmin=ab1[1], xmax=ab1[2]), color="black", fill="grey40", alpha=0.8)+
  geom_rect(mapping=aes(ymin=min(pixel.melt$Reflectance),ymax=max(pixel.melt$Reflectance), xmin=ab2[1], xmax=ab2[2]), color="black", fill="grey40", alpha=0.8)+
  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
                      labels = c("Water","Tree","Grass","Soil","Roof","Road"))+
  labs(color = "Cover Type")+
  ggtitle("Land cover spectral signatures with\n atmospheric absorption bands greyed out")+
  theme(plot.title = element_text(hjust = 0.5, size=20))+
  xlab("Wavelength")

Now we can clearly see that the noisy sections of each spectral signature are within the atmospheric absorption bands. For our final step, let's take all reflectance values from within each absorption band and set them to NA to remove the noisiest sections from the plot.

# Duplicate the spectral signatures into a new data.frame
pixel.melt.masked <- pixel.melt

# Mask out all values within each of the two atmospheric absorption bands
pixel.melt.masked[pixel.melt.masked$wavelengths>ab1[1]&pixel.melt.masked$wavelengths<ab1[2],]$Reflectance <- NA
pixel.melt.masked[pixel.melt.masked$wavelengths>ab2[1]&pixel.melt.masked$wavelengths<ab2[2],]$Reflectance <- NA

# Plot the masked spectral signatures
ggplot()+
  geom_line(data = pixel.melt.masked, mapping = aes(x=wavelengths, y=Reflectance, color=variable), lwd=1.5)+
  scale_colour_manual(values = c("blue3","green4","green2","tan4","grey50","black"),
                      labels = c("Water","Tree","Grass","Soil","Roof","Road"))+
  labs(color = "Cover Type")+
  ggtitle("Land cover spectral signatures with\n atmospheric absorption bands removed")+
  theme(plot.title = element_text(hjust = 0.5, size=20))+
  xlab("Wavelength")

There you have it, spectral signatures for six different land cover types, with the regions from the atmospheric absorption bands removed.

In this tutorial, we will subset an existing HDF5 file containing NEON hyperspectral data. The purpose of this exercise is to generate a smaller file for use in subsequent analyses to reduce the file transfer time and processing power needed.

Why subset your dataset?

There are many reasons why you may wish to subset your HDF5 file. Primarily, HDF5 files may contain a large amount of information that is not necessary for your purposes. By subsetting the file, you can reduce file size, thereby shrinking your storage needs, shortening file transfer/download times, and reducing your processing time for analysis. In this example, we will take a full HDF5 file of NEON hyperspectral reflectance data from the San Joaquin Experimental Range (SJER) that has a file size of ~652 Mb and make a new HDF5 file with a reduced spatial extent, and a reduced spectral resolution, yielding a file of only ~50.1 Mb. This reduction in file size will make it easier and faster to conduct your analysis and to share your data with others. We will then use this subsetted file in the Introduction to Hyperspectral Remote Sensing Data series.

If you find that downloading the full 652 Mb file takes too much time or storage space, you will find a link to this subsetted file at the top of each lesson in the Introduction to Hyperspectral Remote Sensing Data series.

Exploring the NEON hyperspectral HDF5 file structure

In order to determine what information that we want to put into our subset, we should first take a look at the full NEON hyperspectral HDF5 file structure to see what is included. To do so, we will load the required package for this tutorial (you can un-comment the middle two lines to load 'BiocManager' and 'rhdf5' if you don't already have it on your computer).

# Install rhdf5 package (only need to run if not already installed)
# install.packages("BiocManager")
# BiocManager::install("rhdf5")

# Load required packages
library(rhdf5)

Next, we define our working directory where we have saved the full HDF5 file of NEON hyperspectral reflectance data from the SJER site. Note, the filepath to the working directory will depend on your local environment. Then, we create a string (f) of the HDF5 filename and read its attributes.

# set working directory to ensure R can find the file we wish to import and where
# we want to save our files. Be sure to move the download into your working directory!
wd <- "~/Documents/data/" # This will depend on your local environment
setwd(wd)

# Make the name of our HDF5 file a variable
f_full <- paste0(wd,"NEON_D17_SJER_DP3_257000_4112000_reflectance.h5")

Next, let's take a look at the structure of the full NEON hyperspectral reflectance HDF5 file.

View(h5ls(f_full, all=T))

Wow, there is a lot of information in there! The majority of the groups contained within this file are Metadata, most of which are used for processing the raw observations into the reflectance product that we want to use. For demonstration and teaching purposes, we will not need this information. What we will need are things like the Coordinate_System information (so that we can georeference these data), the Wavelength dataset (so that we can match up each band with its appropriate wavelength in the electromagnetic spectrum), and of course the Reflectance_Data themselves. You can also see that each group and dataset has a number of associated attributes (in the 'num_attrs' column). We will want to copy over those attributes into the data subset as well. But first, we need to define each of the groups that we want to populate in our new HDF5 file.

Create new HDF5 file framework

In order to make a new subset HDF5 file, we first must create an empty file with the appropriate name, then we will begin to fill in that file with the essential data and attributes that we want to include. Note that the function h5createFile() will not overwrite an existing file. Therefore, if you have already created or downloaded this file, the function will throw an error! Each function should return 'TRUE' if it runs correctly.

# First, create a name for the new file
f <- paste0(wd, "NEON_hyperspectral_tutorial_example_subset.h5")

# create hdf5 file
h5createFile(f)

## [1] TRUE

# Now we create the groups that we will use to organize our data
h5createGroup(f, "SJER/")

## [1] TRUE

h5createGroup(f, "SJER/Reflectance")

## [1] TRUE

h5createGroup(f, "SJER/Reflectance/Metadata")

## [1] TRUE

h5createGroup(f, "SJER/Reflectance/Metadata/Coordinate_System")

## [1] TRUE

h5createGroup(f, "SJER/Reflectance/Metadata/Spectral_Data")

## [1] TRUE

Adding group attributes

One of the great things about HDF5 files is that they can contain data and attributes within the same group. As explained within the Introduction to Working with HDF5 Format in R series, attributes are a type of metadata that are associated with an HDF5 group or dataset. There may be multiple attributes associated with each group and/or dataset. Attributes come with a name and an associated array of information. In this tutorial, we will read the existing attribute data from the full hyperspectral tile using the h5readAttributes() function (which returns a list of attributes), then we loop through those attributes and write each attribute to its appropriate group using the h5writeAttribute() function.

First, we will do this for the low-level "SJER/Reflectance" group. In this step, we are adding attributes to a group rather than a dataset. To do so, we must first open a file and group interface using the H5Fopen and H5Gopen functions, then we can use h5writeAttribute() to edit the group that we want to give an attribute.

a <- h5readAttributes(f_full,"/SJER/Reflectance/")
fid <- H5Fopen(f)
g <- H5Gopen(fid, "SJER/Reflectance")

for(i in 1:length(names(a))){
  h5writeAttribute(attr = a[[i]], h5obj=g, name=names(a[i]))
}

# It's always a good idea to close the file connection immidiately
# after finishing each step that leaves an open connection.
h5closeAll()

Next, we will loop through each of the datasets within the Coordinate_System group, and copy those (and their attributes, if present) from the full tile to our subset file. The Coordinate_System group contains many important pieces of information for geolocating our data, so we need to make sure that the subset file has that information.

# make a list of all groups within the full tile file
ls <- h5ls(f_full,all=T)

# make a list of all of the names within the Coordinate_System group
cg <- unique(ls[ls$group=="/SJER/Reflectance/Metadata/Coordinate_System",]$name)

# Loop through the list of datasets that we just made above
for(i in 1:length(cg)){
  print(cg[i])
  
  # Read the inividual dataset within the Coordinate_System group
  d=h5read(f_full,paste0("/SJER/Reflectance/Metadata/Coordinate_System/",cg[i]))

  # Read the associated attributes (if any)
  a=h5readAttributes(f_full,paste0("/SJER/Reflectance/Metadata/Coordinate_System/",cg[i]))
    
  # Assign the attributes (if any) to the dataset
  attributes(d)=a
  
  # Finally, write the dataset to the HDF5 file
  h5write(obj=d,file=f,
          name=paste0("/SJER/Reflectance/Metadata/Coordinate_System/",cg[i]),
          write.attributes=T)
}

## [1] "Coordinate_System_String"
## [1] "EPSG Code"
## [1] "Map_Info"
## [1] "Proj4"

Spectral Subsetting

The goal of subsetting this dataset is to substantially reduce the file size, making it faster to download and process these data. While each AOP mosaic tile is not particularly large in terms of its spatial scale (1km by 1km at 1m resolution= 1,000,000 pixels, or about half as many pixels at shown on a standard 1080p computer screen), the 426 spectral bands available result in a fairly large file size. Therefore, we will reduce the spectral resolution of these data by selecting every fourth band in the dataset, which reduces the file size to 1/4 of the original!

Some circumstances demand the full spectral resolution file. For example, if you wanted to discern between the spectral signatures of similar minerals, or if you were conducting an analysis of the differences in the 'red edge' between plant functional types, you would want to use the full spectral resolution of the original hyperspectral dataset. Still, we can use far fewer bands for demonstration and teaching purposes, while still getting a good sense of what these hyperspectral data can do.

# First, we make our 'index', a list of number that will allow us to select every fourth band, using the "sequence" function seq()
idx <- seq(from = 1, to = 426, by = 4)

# We then use this index to select particular wavelengths from the full tile using the "index=" argument
wavelengths <- h5read(file = f_full, 
             name = "SJER/Reflectance/Metadata/Spectral_Data/Wavelength", 
             index = list(idx)
            )

# As per above, we also need the wavelength attributes
wavelength.attributes <- h5readAttributes(file = f_full, 
                       name = "SJER/Reflectance/Metadata/Spectral_Data/Wavelength")
attributes(wavelengths) <- wavelength.attributes

# Finally, write the subset of wavelengths and their attributes to the subset file
h5write(obj=wavelengths, file=f,
        name="SJER/Reflectance/Metadata/Spectral_Data/Wavelength",
        write.attributes=T)

Spatial Subsetting

Even after spectral subsetting, our file size would be greater than 100Mb. herefore, we will also perform a spatial subsetting process to further reduce our file size. Now, we need to figure out which part of the full image that we want to extract for our subset. It takes a number of steps in order to read in a band of data and plot the reflectance values - all of which are thoroughly described in the Intro to Working with Hyperspectral Remote Sensing Data in HDF5 Format in R tutorial. For now, let's focus on the essentials for our problem at hand. In order to explore the spatial qualities of this dataset, let's plot a single band as an overview map to see what objects and land cover types are contained within this mosaic tile. The Reflectance_Data dataset has three dimensions in the order of bands, columns, rows. We want to extract a single band, and all 1,000 columns and rows, so we will feed those values into the index= argument as a list. For this example, we will select the 58th band in the hyperspectral dataset, which corresponds to a wavelength of 667nm, which is in the red end of the visible electromagnetic spectrum. We will use NULL in the column and row position to indicate that we want all of the columns and rows (we agree that it is weird that NULL indicates "all" in this circumstance, but that is the default behavior for this, and many other, functions).

# Extract or "slice" data for band 58 from the HDF5 file
b58 <- h5read(f_full,name = "SJER/Reflectance/Reflectance_Data",
             index=list(58,NULL,NULL))
h5closeAll()

# convert from array to matrix
b58 <- b58[1,,]

# Make a plot to view this band
image(log(b58), col=grey(0:100/100))

As we can see here, this hyperspectral reflectance tile contains a school campus that is under construction. There are many different land cover types contained here, which makes it a great example! Perhaps the most unique feature shown is in the bottom right corner of this image, where we can see the tip of a small reservoir. Let's be sure to capture this feature in our spatial subset, as well as a few other land cover types (irrigated grass, trees, bare soil, and buildings).

While raster images count their pixels from the top left corner, we are working with a matrix, which counts its pixels from the bottom left corner. Therefore, rows are counted from the bottom to the top, and columns are counted from the left to the right. If we want to sample the bottom right quadrant of this image, we need to select rows 1 through 500 (bottom half), and columns 501 through 1000 (right half). Note that, as above, the index= argument in h5read() requires a list of three dimensions for this example - in the order of bands, columns, rows.

subset_rows <- 1:500
subset_columns <- 501:1000
# Extract or "slice" data for band 44 from the HDF5 file
b58 <- h5read(f_full,name = "SJER/Reflectance/Reflectance_Data",
             index=list(58,subset_columns,subset_rows))
h5closeAll()

# convert from array to matrix
b58 <- b58[1,,]

# Make a plot to view this band
image(log(b58), col=grey(0:100/100))

Perfect - now we have a spatial subset that includes all of the different land cover types that we are interested in investigating.

Extracting a subset

Now that we have determined our ideal spectral and spatial subsets for our analysis, we are ready to put both of those pieces of information into our h5read() function to extract our example subset out of the full NEON hyperspectral dataset. Here, we are taking every fourth band (using our idx variabe), columns 501:1000 (the right half of the tile) and rows 1:500 (the bottom half of the tile). The results in us extracting every fourth band of the bottom-right quadrant of the mosaic tile.

# Read in reflectance data.
# Note the list that we feed into the index argument! 
# This tells the h5read() function which bands, rows, and 
# columns to read. This is ultimately how we reduce the file size.
hs <- h5read(file = f_full, 
             name = "SJER/Reflectance/Reflectance_Data", 
             index = list(idx, subset_columns, subset_rows)
            )

As per the 'adding group attributes' section above, we will need to add the attributes to the hyperspectral data (hs) before writing to the new HDF5 subset file (f). The hs variable already has one attribute, $dim, which contains the actual dimensions of the hs array, and will be important for writing the array to the f file later. We will want to combine this attribute with all of the other Reflectance_Data group attributes from the original HDF5 file, f. However, some of the attributes will no longer be valid, such as the Dimensions and Spatial_Extent_meters attributes, so we will need to overwrite those before assigning these attributes to the hs variable to then write to the f file.

# grab the '$dim' attribute - as this will be needed 
# when writing the file at the bottom
hsd <- attributes(hs)

# We also need the attributes for the reflectance data.
ha <- h5readAttributes(file = f_full, 
                       name = "SJER/Reflectance/Reflectance_Data")

# However, some of the attributes are not longer valid since 
# we changed the spatial extend of this dataset. therefore, 
# we will need to overwrite those with the correct values.
ha$Dimensions <- c(500,500,107) # Note that the HDF5 file saves dimensions in a different order than R reads them
ha$Spatial_Extent_meters[1] <- ha$Spatial_Extent_meters[1]+500
ha$Spatial_Extent_meters[3] <- ha$Spatial_Extent_meters[3]+500
attributes(hs) <- c(hsd,ha)

# View the combined attributes to ensure they are correct
attributes(hs)

## $dim
## [1] 107 500 500
## 
## $Cloud_conditions
## [1] "For cloud conditions information see Weather Quality Index dataset."
## 
## $Cloud_type
## [1] "Cloud type may have been selected from multiple flight trajectories."
## 
## $Data_Ignore_Value
## [1] -9999
## 
## $Description
## [1] "Atmospherically corrected reflectance."
## 
## $Dimension_Labels
## [1] "Line, Sample, Wavelength"
## 
## $Dimensions
## [1] 500 500 107
## 
## $Interleave
## [1] "BSQ"
## 
## $Scale_Factor
## [1] 10000
## 
## $Spatial_Extent_meters
## [1]  257500  258000 4112500 4113000
## 
## $Spatial_Resolution_X_Y
## [1] 1 1
## 
## $Units
## [1] "Unitless."
## 
## $Units_Valid_range
## [1]     0 10000

# Finally, write the hyperspectral data, plus attributes, 
# to our new file 'f'.
h5write(obj=hs, file=f,
        name="SJER/Reflectance/Reflectance_Data",
        write.attributes=T)

## You created a large dataset with compression and chunking.
## The chunk size is equal to the dataset dimensions.
## If you want to read subsets of the dataset, you should testsmaller chunk sizes to improve read times.

# It's always a good idea to close the HDF5 file connection
# before moving on.
h5closeAll()

That's it! We just created a subset of the original HDF5 file, and included the most essential groups and metadata for exploratory analysis. You may consider adding other information, such as the weather quality indicator, when subsetting datasets for your own purposes.

If you want to take a look at the subset that you just made, run the h5ls() function:

View(h5ls(f, all=T))

This tutorial explores NEON geolocation data. The focus is on the locations of NEON observational sampling and sensor data; NEON remote sensing data are inherently spatial and have dedicated tutorials. If you are interested in connecting remote sensing with ground-based measurements, the methods in the vegetation structure and canopy height model tutorial can be generalized to other data products.

In planning your analyses, consider what level of spatial resolution is required. There is no reason to carefully map each measurement if precise spatial locations aren't required to address your hypothesis! For example, if you want to use the Vegetation structure data product to calculate a site-scale estimate of biomass and production, the spatial coordinates of each tree are probably not needed. If you want to explore relationships between vegetation and beetle communities, you will need to identify the sampling plots where NEON measures both beetles and vegetation, but finer-scale coordinates may not be needed. Finally, if you want to relate vegetation measurements to airborne remote sensing data, you will need very accurate coordinates for each measurement on the ground.

Setup R Environment

We'll need several R packages in this tutorial. Install the packages, if not already installed, and load the libraries for each. 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.

# run once to get the package, and re-run if you need to get updates

install.packages("ggplot2")

install.packages("dplyr")

install.packages("neonUtilities")

install.packages("neonOS")

install.packages("devtools")

devtools::install_github("NEONScience/NEON-geolocation/geoNEON")



library(ggplot2)

library(dplyr)

library(neonUtilities)

library(neonOS)

library(geoNEON)



token <- Sys.getenv("NEON_TOKEN")

Locations for observational data

Plot level locations

Both aquatic and terrestrial observational data downloads include spatial data in the downloaded files. The spatial data in the aquatic data files are the most precise locations available for the sampling events. The spatial data in the terrestrial data downloads represent the locations of the sampling plots. In some cases, the plot is the most precise location available, but for many terrestrial data products, more precise locations can be calculated for specific sampling events.

Here, we'll download the Vegetation structure (DP1.10098.001) data product, examine the plot location data in the download, then calculate the locations of individual trees. These steps can be extrapolated to other terrestrial observational data products; the specific sampling layout varies from data product to data product, but the methods for working with the data are similar.

First, let's download the vegetation structure data from one site, Wind River Experimental Forest (WREF). If downloading data using the neonUtilities package is new to you, check out the Download and Explore tutorial.

vst <- loadByProduct(dpID="DP1.10098.001", 
                     site="WREF",
                     include.provisional=T,
                     check.size=F,
                     token=token)

Data downloaded this way are stored in R as a large list. For this tutorial, we'll work with the individual dataframes within this large list. Alternatively, each dataframe can be assigned as its own object.

To find the spatial data for any given data product, view the variables files to figure out which data table the spatial data are contained in.

View(vst$variables_10098)

Looking through the variables, we can see that the spatial data (decimalLatitude and decimalLongitude, etc) are in the vst_perplotperyear table. Let's take a look at the table.

View(vst$vst_perplotperyear)

As noted above, the spatial data here are at the plot level; the latitude and longitude represent the centroid of the sampling plot. We can map these plots on the landscape using the easting and northing variables; these are the UTM coordinates. At this site, tower plots are 40 m x 40 m, and distributed plots are 20 m x 20 m; we can use the symbols() function to draw boxes of the correct size.

We'll also use the treesPresent variable to subset to only those plots where trees were found and measured.

# start by subsetting data to plots with trees

vst.trees <- vst$vst_perplotperyear |>
  filter(treesPresent=="Y")



# create plot size variable

vst.trees <- vst.trees |>
  mutate(plotsize = case_when(plotType=="tower" ~ 40,
                              plotType=="distributed" ~ 20,
                              TRUE ~ NA))



# create map of plots

symbols(vst.trees$easting,
        vst.trees$northing,
        squares=vst.trees$plotsize, 
        inches=F,
        xlab="Easting", 
        ylab="Northing")

We can see where the plots are located across the landscape, and we can see the denser cluster of plots in the area near the micrometeorology tower.

For many analyses, this level of spatial data may be sufficient. Calculating the precise location of each tree is only required for certain hypotheses; consider whether you need these data when working with a data product with plot-level spatial data.

Looking back at the variables_10098 table, notice that there is a table in this data product called vst_mappingandtagging, suggesting we can find mapping data there. Let's take a look.

View(vst$vst_mappingandtagging)

Here we see data fields for stemDistance and stemAzimuth. Looking back at the variables_10098 file, we see these fields contain the distance and azimuth from a pointID to a specific stem. To calculate the precise coordinates of each tree, we would need to get the locations of the pointIDs, and then adjust the coordinates based on distance and azimuth. The Data Product User Guide describes how to carry out these steps, and can be downloaded from the Data Product Details page.

However, carrying out these calculations yourself is not the only option! The geoNEON package contains a function that can do this for you, for the TOS data products with location data more precise than the plot level.

Sampling locations

The getLocTOS() function in the geoNEON package uses the NEON API to access NEON location data and then makes protocol-specific calculations to return precise locations for each sampling effort. This function works for a subset of NEON TOS data products. The list of tables and data products that can be entered is in the package documentation on GitHub.

For more information about the NEON API, see the API tutorial and the API web page. For more information about the location calculations used in each data product, see the Data Product User Guide for each product.

The getLocTOS() function requires three inputs:

• A data table that contains spatial data from a NEON TOS data product
• The NEON table name of that data table
• NEON API token

For vegetation structure locations, the function call looks like this. This function may take a while to download all the location data.

vst.loc <- getLocTOS(data=vst$vst_mappingandtagging,
                     dataProd="vst_mappingandtagging",
                     token=token)

What additional data are now available that were obtained and calculated by getLocTOS()? Let's see what new column names have been added to the table.

setdiff(names(vst.loc),
        names(vst$vst_mappingandtagging))

##  [1] "utmZone"                  "adjNorthing"              "adjEasting"               "adjCoordinateUncertainty"
##  [5] "adjDecimalLatitude"       "adjDecimalLongitude"      "adjElevation"             "adjElevationUncertainty" 
##  [9] "locationCurrent"          "locationStartDate"        "locationEndDate"

Now we have adjusted latitude, longitude, and elevation, and the corresponding easting and northing UTM data. We also have coordinate uncertainty data for these coordinates.

As we did with the plots above, we can use the easting and northing data to plot the locations of the individual trees.

plot(vst.loc$adjEasting, 
     vst.loc$adjNorthing,
     xlab="Easting", 
     ylab="Northing",
     pch=".")

We can see the mapped trees in the same plots we mapped above. We've plotted each individual tree as a ., so all we can see at this scale is the cluster of dots that make up each plot.

Let's zoom in on a single plot:

plot(subset(vst.loc, plotID=="WREF_085",
            select=c(adjEasting, adjNorthing)),
     pch=20, xlab="Easting", ylab="Northing")

Now we can see the location of each tree within the sampling plot WREF_085. This is interesting, but it would be more interesting if we could see more information about each tree. How are species distributed across the plot, for instance?

We can plot the tree species at each location using the text() function and the taxonID field.

vst.085 <- vst.loc |>
  filter(plotID=="WREF_085")



plot(vst.085$adjEasting,
     vst.085$adjNorthing,
     xlab="Easting", 
     ylab="Northing",
     type="n")

text(vst.085$adjEasting,
     vst.085$adjNorthing,
     labels=vst.085$taxonID,
     cex=0.5)

Almost all of the mapped trees in this plot are either Pseudotsuga menziesii or Tsuga heterophylla (Douglas fir and Western hemlock), not too surprising at Wind River.

But suppose we want to map the diameter of each tree? This is a very common way to present a stem map, it gives a visual as if we were looking down on the plot from overhead and had cut off each tree at its measurement height.

Other than taxon, the attributes of the trees, such as diameter, height, growth form, and canopy position, are found in the vst_apparentindividual table, not in the vst_mappingandtagging table. We'll need to join the two tables to get the tree attributes together with their mapped locations.

The neonOS package contains the function joinTableNEON(), which can be used to do this. See the tutorial for the neonOS package for more details about this function.

veg <- joinTableNEON(vst.loc, 
                     vst$vst_apparentindividual,
                     name1="vst_mappingandtagging",
                     name2="vst_apparentindividual")

Now we can use the symbols() function to plot the diameter of each tree, at its spatial coordinates, to create a correctly scaled map of boles in the plot. Note that stemDiameter is in centimeters, while easting and northing UTMs are in meters, so we divide by 100 to scale correctly.

veg.085 <- veg |>
  filter(plotID=="WREF_085")



symbols(veg.085$adjEasting, 
        veg.085$adjNorthing, 
        circles=veg.085$stemDiameter/100/2, 
        inches=F, 
        xlab="Easting", 
        ylab="Northing")

If you are interested in taking the vegetation structure data a step further, and connecting measurements of trees on the ground to remotely sensed Lidar data, check out the Vegetation Structure and Canopy Height Model tutorial.

If you are interested in working with other terrestrial observational (TOS) data products, the basic techniques used here to find precise sampling locations and join data tables can be adapted to other TOS data products. Consult the Data Product User Guide for each data product to find details specific to that data product.

Locations for sensor data

Downloads of instrument system (IS) data include a file called sensor_positions.csv. The sensor positions file contains information about the coordinates of each sensor, relative to a reference location.

While the specifics vary, techniques are generalizable for working with sensor data and the sensor_positions.csv file. For this tutorial, let's look at the sensor locations for soil temperature (PAR; DP1.00041.001) at
the NEON Treehaven site (TREE) in July 2018. To reduce our file size, we'll use the 30 minute averaging interval. Our final product from this section is to create a depth profile of soil temperature in one soil plot.

If downloading data using the neonUtilties package is new to you, check out the neonUtilities tutorial.

This function will download about 7 MB of data as written so we have check.size =F for ease of running the code.

soilT <- loadByProduct(dpID="DP1.00041.001", site="TREE",
                    startdate="2018-07", enddate="2018-07",
                    timeIndex=30, check.size=F,
                    token=token)

Sensor positions file

Now we can specifically look at the sensor positions file.

# create object for sensor positions file

pos <- soilT$sensor_positions_00041



# view column names

names(pos)

##  [1] "domainID"                         "siteID"                           "HOR.VER"                         
##  [4] "sensorLocationID"                 "sensorLocationDescription"        "positionStartDateTime"           
##  [7] "positionEndDateTime"              "referenceLocationID"              "referenceLocationIDDescription"  
## [10] "referenceLocationIDStartDateTime" "referenceLocationIDEndDateTime"   "xOffset"                         
## [13] "yOffset"                          "zOffset"                          "pitch"                           
## [16] "roll"                             "azimuth"                          "locationReferenceLatitude"       
## [19] "locationReferenceLongitude"       "locationReferenceElevation"       "eastOffset"                      
## [22] "northOffset"                      "xAzimuth"                         "yAzimuth"                        
## [25] "publicationDate"                  "release"

# view table

View(pos)

The sensor locations are indexed by the HOR.VER variable - see the file naming conventions page for more details.

Using unique() we can view all the location indices in this file.

unique(pos$HOR.VER)

##  [1] "001.501" "001.502" "001.503" "001.504" "001.505" "001.506" "001.507" "001.508" "001.509" "002.501" "002.502" "002.503"
## [13] "002.504" "002.505" "002.506" "002.507" "002.508" "002.509" "003.501" "003.502" "003.503" "003.504" "003.505" "003.506"
## [25] "003.507" "003.508" "003.509" "004.501" "004.502" "004.503" "004.504" "004.505" "004.506" "004.507" "004.508" "004.509"
## [37] "005.501" "005.502" "005.503" "005.504" "005.505" "005.506" "005.507" "005.508" "005.509"

Soil temperature data are collected in 5 instrumented soil plots inside the tower footprint. We see this reflected in the data where HOR = 001 to 005. Within each plot, temperature is measured at 9 depths, seen in VER = 501 to 509. At some sites, the number of depths may differ slightly.

The x, y, and z offsets in the sensor positions file are the relative distance, in meters, to the reference latitude, longitude, and elevation in the file.

The HOR and VER indices in the sensor positions file correspond to the verticalPosition and horizontalPosition fields in soilT$ST_30_minute.

Note that there are two sets of position data for soil plot 001, and that one set has a positionEndDateTime date in the file. This indicates sensors either moved or were relocated; in this case there was a frost heave incident. You can read about it in the issue log, which is displayed on the Data Product Details page, and also included as a table in the data download:

soilT$issueLog_00041 |>
  filter(grepl("TREE soil plot 1", locationAffected))

##     id parentIssueID            issueDate         resolvedDate       dateRangeStart         dateRangeEnd
## 1 9328            NA 2019-05-23T00:00:00Z 2019-05-23T00:00:00Z 2018-11-07T00:00:00Z 2019-04-19T00:00:00Z
##                                                                                                                         locationAffected
## 1 D05 TREE soil plot 1 measurement levels 1-9 (HOR.VER: 001.501, 001.502, 001.503, 001.504, 001.505, 001.506, 001.507, 001.508, 001.509)
##                                                                                                                                                                                                                          issue
## 1 Soil temperature sensors were pushed or pulled out of the ground by 3 cm over winter, presumably due to freeze-thaw action. The exact timing of this is unknown, but it occurred sometime between 2018-11-07 and 2019-04-19.
##                                                                                       resolution
## 1 Sensor depths were updated in the database with a start date of 2018-11-07 for the new depths.

Since we're working with data from July 2018, and the change in sensor locations is dated Nov 2018, we'll use the original locations. There are a number of ways to drop the later locations from the table; here, here, we find the rows in which the positionEndDateTime field is empty, indicating no end date, and the rows corresponding to soil plot 001, and drop all the rows that meet both criteria.

pos <- pos |>
  filter_out(grepl("001.", pos$HOR.VER) & 
               is.na(positionEndDateTime))

Our goal is to plot a time series of temperature, stratified by depth, so let's start by joining the data file and sensor positions file, to bring the depth measurements into the same data frame with the data.

# paste horizontalPosition and verticalPosition together

# to match HOR.VER in the sensor positions file

soilT$ST_30_minute <- soilT$ST_30_minute |>
  mutate(HOR.VER = paste(horizontalPosition,
                         verticalPosition,
                         sep="."))



# left join to keep all temperature records

soilTHV <- merge(soilT$ST_30_minute, pos, 
                 by="HOR.VER", all.x=T)

And now we can plot soil temperature over time for each depth. We'll use ggplot since it's well suited to this kind of stratification. Each soil plot is its own panel, and each depth is its own line:

gg <- ggplot(soilTHV, 
             aes(endDateTime, soilTempMean, 
                 group=zOffset, color=zOffset)) +
             geom_line() + 
        facet_wrap(~horizontalPosition)

gg

## Warning: Removed 1488 rows containing missing values or values outside the scale range (`geom_line()`).

We can see that as soil depth increases, temperatures become much more stable, while the shallowest measurement has a clear diurnal cycle. We can also see that something has gone wrong with one of the sensors in plot 002. To remove those data, use only values where the final quality flag passed, i.e. finalQF = 0

gg <- ggplot(subset(soilTHV, finalQF==0), 
             aes(endDateTime, soilTempMean, 
                 group=zOffset, color=zOffset)) +
             geom_line() + 
        facet_wrap(~horizontalPosition)

gg

This data tutorial provides an introduction to working with NEON eddy flux data, using the neonUtilities R package or the neonutilities Python package. If you are new to NEON data, we recommend starting with a more general tutorial, such as the Download and Explore tutorial. Some of the functions and techniques described there will be used here, as well as functions and data formats that are unique to the eddy flux system.

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.

This tutorial assumes general familiarity with eddy flux data and associated concepts.

install.packages('BiocManager')
BiocManager::install('rhdf5')
install.packages('neonUtilities')
install.packages('ggplot2')

library(neonUtilities)
library(ggplot2)
token <- Sys.getenv("NEON_TOKEN")

pip install neonutilities

import neonutilities as nu
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from datetime import datetime
import dotenv
import os
dotenv.load_dotenv()

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

zipsByProduct(dpID="DP4.00200.001", package="basic", 
              site=c("NIWO", "HARV"), 
              startdate="2018-06", enddate="2018-07",
              savepath=getwd(), 
              check.size=F,
              token=token)

nu.zips_by_product(dpid="DP4.00200.001", package="basic", 
              site=["NIWO","HARV"], 
              startdate="2018-06", enddate="2018-07",
              savepath=os.getcwd(), 
              check_size=False,
              token=token)

flux <- stackEddy(filepath=paste0(getwd(), "/filesToStack00200"),
                 level="dp04")

## Extracting data

## Joining data variables by file

## Stacking files by month

## Getting metadata tables

names(flux)

## [1] "HARV"                        "NIWO"                       
## [3] "variables"                   "objDesc"                    
## [5] "issueLog"                    "citation_00200_RELEASE-2026"

head(flux$NIWO)

##               timeBgn             timeEnd data.fluxCo2.nsae.flux
## 1 2018-06-01 00:00:00 2018-06-01 00:29:59              0.1895604
## 2 2018-06-01 00:30:00 2018-06-01 00:59:59              0.9429929
## 3 2018-06-01 01:00:00 2018-06-01 01:29:59              0.5101887
## 4 2018-06-01 01:30:00 2018-06-01 01:59:59              0.7980803
## 5 2018-06-01 02:00:00 2018-06-01 02:29:59              0.4345524
## 6 2018-06-01 02:30:00 2018-06-01 02:59:59              0.8440311
##   data.fluxCo2.stor.flux data.fluxCo2.turb.flux data.fluxH2o.nsae.flux
## 1            -0.06349451              0.2530549              16.428558
## 2             0.09749592              0.8454970               7.965026
## 3             0.02525844              0.4849303               5.387009
## 4             0.23341588              0.5646644               9.933795
## 5             0.17415575              0.2603966               2.963931
## 6            -0.01496958              0.8590007               4.488221
##   data.fluxH2o.stor.flux data.fluxH2o.turb.flux data.fluxMome.turb.veloFric
## 1             3.34458774              13.083971                   0.2025301
## 2            -1.51509573               9.480122                   0.1922362
## 3            -4.26233930               9.649348                   0.1197608
## 4             0.15083070               9.782964                   0.1190342
## 5            -0.01668207               2.980613                   0.1589262
## 6            -0.52273792               5.010959                   0.1115558
##   data.fluxTemp.nsae.flux data.fluxTemp.stor.flux data.fluxTemp.turb.flux
## 1               4.7823169              -1.4575674               6.2398843
## 2              -0.5981044               0.3404014              -0.9385058
## 3              -4.3428949               0.1870749              -4.5299698
## 4             -12.6516137              -2.4904875             -10.1611262
## 5              -5.0544306              -0.7514722              -4.3029584
## 6              -7.6999704              -1.9047428              -5.7952276
##   data.foot.stat.angZaxsErth data.foot.stat.distReso
## 1                    94.2262                    8.34
## 2                   355.4252                    8.34
## 3                   359.8013                    8.34
## 4                   137.7743                    8.34
## 5                   188.4800                    8.34
## 6                   183.1920                    8.34
##   data.foot.stat.veloYaxsHorSd data.foot.stat.veloZaxsHorSd
## 1                    0.7955893                    0.2713232
## 2                    0.8590177                    0.2300000
## 3                    1.2601763                    0.2300000
## 4                    0.7332641                    0.2300000
## 5                    0.7096286                    0.2300000
## 6                    0.3789859                    0.2300000
##   data.foot.stat.veloFric data.foot.stat.distZaxsMeasDisp
## 1               0.2025428                            8.34
## 2               0.2000000                            8.34
## 3               0.2000000                            8.34
## 4               0.2000000                            8.34
## 5               0.2000000                            8.34
## 6               0.2000000                            8.34
##   data.foot.stat.distZaxsRgh data.foot.stat.distObkv data.foot.stat.paraStbl
## 1                 0.04083968               -72.43385            -0.115139545
## 2                 0.29601180               965.24689             0.008640276
## 3                 0.23002973                27.47032             0.303600430
## 4                 0.83400000                10.29068             0.810442357
## 5                 0.83400000                58.07739             0.143601482
## 6                 0.83400000                15.03590             0.554672317
##   data.foot.stat.distZaxsAbl data.foot.stat.distXaxs90
## 1                       1000                    325.26
## 2                       1000                    258.54
## 3                       1000                    275.22
## 4                       1000                    208.50
## 5                       1000                    208.50
## 6                       1000                    208.50
##   data.foot.stat.distXaxsMax data.foot.stat.distYaxs90 qfqm.fluxCo2.nsae.qfFinl
## 1                     133.44                     25.02                        1
## 2                     108.42                     50.04                        1
## 3                     108.42                     75.06                        1
## 4                      83.40                     75.06                        1
## 5                      83.40                     66.72                        1
## 6                      83.40                     41.70                        1
##   qfqm.fluxCo2.stor.qfFinl qfqm.fluxCo2.turb.qfFinl qfqm.fluxH2o.nsae.qfFinl
## 1                        1                        1                        1
## 2                        1                        1                        1
## 3                        1                        1                        1
## 4                        1                        1                        1
## 5                        1                        1                        1
## 6                        1                        1                        1
##   qfqm.fluxH2o.stor.qfFinl qfqm.fluxH2o.turb.qfFinl qfqm.fluxMome.turb.qfFinl
## 1                        1                        1                         0
## 2                        0                        1                         0
## 3                        0                        1                         1
## 4                        0                        1                         1
## 5                        0                        1                         0
## 6                        0                        1                         0
##   qfqm.fluxTemp.nsae.qfFinl qfqm.fluxTemp.stor.qfFinl qfqm.fluxTemp.turb.qfFinl
## 1                         1                         0                         1
## 2                         1                         0                         1
## 3                         1                         0                         1
## 4                         1                         0                         1
## 5                         1                         0                         1
## 6                         1                         0                         1
##   qfqm.foot.turb.qfFinl
## 1                     0
## 2                     0
## 3                     0
## 4                     0
## 5                     0
## 6                     0

term <- unlist(strsplit(names(flux$NIWO), split=".", fixed=T))
flux$objDesc[which(flux$objDesc$X...Object %in% term),]

##           X...Object
## 22       angZaxsErth
## 84              data
## 91          distReso
## 92        distXaxs90
## 93       distXaxsMax
## 94        distYaxs90
## 95       distZaxsAbl
## 100 distZaxsMeasDisp
## 101      distZaxsRgh
## 125             flux
## 126          fluxCo2
## 128          fluxH2o
## 130         fluxMome
## 132         fluxTemp
## 133             foot
## 200             nsae
## 244           qfFinl
## 274             qfqm
## 451             stat
## 453             stor
## 474          timeBgn
## 475          timeEnd
## 477             turb
## 505         veloFric
## 512    veloYaxsHorSd
## 515    veloZaxsHorSd
##                                                                                                           Description
## 22                                                                                                     Wind direction
## 84                                                                                             Represents data fields
## 91                                                                      Resolution of vector of distances (units = m)
## 92                                     Distance of 90% contriubtion of footprint aligned with dominant wind direction
## 93                                        Max distance of footprint contriubtion aligned with dominant wind direction
## 94  Distance of 90% contriubtion of footprint aligned with direction that is perpendicular to dominant wind direction
## 95                                              Atmospheric boundary layer height used in flux footprint calculations
## 100                                                                     Measurement height minus displacement height 
## 101                                                               Roughness length used in flux footprint calculation
## 125                                                                                             General term for flux
## 126                                                                     CO2 flux (can be turbulent, storage, or nsae)
## 128                                      H2O flux (can be turbulent, storage, or nsae)  - relates to latent heat flux
## 130                                                                Momentum flux (can be turbulent, storage, or nsae)
## 132                             Temperature flux (can be turbulent, storage, or nsae) - relates to sensible heat flux
## 133                                                                                                         Footprint
## 200                                                       Net surface atmospher exchange (turbulent + storage fluxes)
## 244         The final quality flag indicating if the data are valid for the given aggregation period (1=fail, 0=pass)
## 274   Quality flag and quality metrics, represents quality flags and quality metrics that accompany the provided data
## 451                                                                                                        Statistics
## 453                                                                                                           Storage
## 474                                                                      The beginning time of the aggregation period
## 475                                                                            The end time of the aggregation period
## 477                                                                                                         Turbulent
## 505                                                                           Friction velocity/ustar (units = m s-1)
## 512                                        Standard deviation of horizontal wind speed in Y direction (units = m s-1)
## 515                                        Standard deviation of horizontal wind speed in X direction (units = m s-1)

flux$variables

##    category   system variable             stat           units
## 1      data  fluxCo2     nsae          timeBgn              NA
## 2      data  fluxCo2     nsae          timeEnd              NA
## 3      data  fluxCo2     nsae             flux umolCo2 m-2 s-1
## 4      data  fluxCo2     stor          timeBgn              NA
## 5      data  fluxCo2     stor          timeEnd              NA
## 6      data  fluxCo2     stor             flux umolCo2 m-2 s-1
## 7      data  fluxCo2     turb          timeBgn              NA
## 8      data  fluxCo2     turb          timeEnd              NA
## 9      data  fluxCo2     turb             flux umolCo2 m-2 s-1
## 10     data  fluxH2o     nsae          timeBgn              NA
## 11     data  fluxH2o     nsae          timeEnd              NA
## 12     data  fluxH2o     nsae             flux           W m-2
## 13     data  fluxH2o     stor          timeBgn              NA
## 14     data  fluxH2o     stor          timeEnd              NA
## 15     data  fluxH2o     stor             flux           W m-2
## 16     data  fluxH2o     turb          timeBgn              NA
## 17     data  fluxH2o     turb          timeEnd              NA
## 18     data  fluxH2o     turb             flux           W m-2
## 19     data fluxMome     turb          timeBgn              NA
## 20     data fluxMome     turb          timeEnd              NA
## 21     data fluxMome     turb         veloFric           m s-1
## 22     data fluxTemp     nsae          timeBgn              NA
## 23     data fluxTemp     nsae          timeEnd              NA
## 24     data fluxTemp     nsae             flux           W m-2
## 25     data fluxTemp     stor          timeBgn              NA
## 26     data fluxTemp     stor          timeEnd              NA
## 27     data fluxTemp     stor             flux           W m-2
## 28     data fluxTemp     turb          timeBgn              NA
## 29     data fluxTemp     turb          timeEnd              NA
## 30     data fluxTemp     turb             flux           W m-2
## 31     data     foot     stat          timeBgn              NA
## 32     data     foot     stat          timeEnd              NA
## 33     data     foot     stat      angZaxsErth             deg
## 34     data     foot     stat         distReso               m
## 35     data     foot     stat    veloYaxsHorSd           m s-1
## 36     data     foot     stat    veloZaxsHorSd           m s-1
## 37     data     foot     stat         veloFric           m s-1
## 38     data     foot     stat distZaxsMeasDisp               m
## 39     data     foot     stat      distZaxsRgh               m
## 40     data     foot     stat         distObkv               m
## 41     data     foot     stat         paraStbl               -
## 42     data     foot     stat      distZaxsAbl               m
## 43     data     foot     stat       distXaxs90               m
## 44     data     foot     stat      distXaxsMax               m
## 45     data     foot     stat       distYaxs90               m
## 46     qfqm  fluxCo2     nsae          timeBgn              NA
## 47     qfqm  fluxCo2     nsae          timeEnd              NA
## 48     qfqm  fluxCo2     nsae           qfFinl              NA
## 49     qfqm  fluxCo2     stor           qfFinl              NA
## 50     qfqm  fluxCo2     stor          timeBgn              NA
## 51     qfqm  fluxCo2     stor          timeEnd              NA
## 52     qfqm  fluxCo2     turb          timeBgn              NA
## 53     qfqm  fluxCo2     turb          timeEnd              NA
## 54     qfqm  fluxCo2     turb           qfFinl              NA
## 55     qfqm  fluxH2o     nsae          timeBgn              NA
## 56     qfqm  fluxH2o     nsae          timeEnd              NA
## 57     qfqm  fluxH2o     nsae           qfFinl              NA
## 58     qfqm  fluxH2o     stor           qfFinl              NA
## 59     qfqm  fluxH2o     stor          timeBgn              NA
## 60     qfqm  fluxH2o     stor          timeEnd              NA
## 61     qfqm  fluxH2o     turb          timeBgn              NA
## 62     qfqm  fluxH2o     turb          timeEnd              NA
## 63     qfqm  fluxH2o     turb           qfFinl              NA
## 64     qfqm fluxMome     turb          timeBgn              NA
## 65     qfqm fluxMome     turb          timeEnd              NA
## 66     qfqm fluxMome     turb           qfFinl              NA
## 67     qfqm fluxTemp     nsae          timeBgn              NA
## 68     qfqm fluxTemp     nsae          timeEnd              NA
## 69     qfqm fluxTemp     nsae           qfFinl              NA
## 70     qfqm fluxTemp     stor           qfFinl              NA
## 71     qfqm fluxTemp     stor          timeBgn              NA
## 72     qfqm fluxTemp     stor          timeEnd              NA
## 73     qfqm fluxTemp     turb          timeBgn              NA
## 74     qfqm fluxTemp     turb          timeEnd              NA
## 75     qfqm fluxTemp     turb           qfFinl              NA
## 76     qfqm     foot     turb          timeBgn              NA
## 77     qfqm     foot     turb          timeEnd              NA
## 78     qfqm     foot     turb           qfFinl              NA

flux = nu.stack_eddy(filepath=os.getcwd() + "/filesToStack00200",
                 level="dp04")

flux.keys()

## dict_keys(['HARV', 'NIWO', 'variables', 'objDesc', 'issueLog', 'citation_00200_RELEASE-2026'])

flux['NIWO'].head

## <bound method NDFrame.head of                        timeBgn  ... qfqm.foot.turb.qfFinl
## 0    2018-06-01 00:00:00+00:00  ...                     0
## 1    2018-06-01 00:30:00+00:00  ...                     0
## 2    2018-06-01 01:00:00+00:00  ...                     0
## 3    2018-06-01 01:30:00+00:00  ...                     0
## 4    2018-06-01 02:00:00+00:00  ...                     0
## ...                        ...  ...                   ...
## 2923 2018-07-31 21:30:00+00:00  ...                     0
## 2924 2018-07-31 22:00:00+00:00  ...                     0
## 2925 2018-07-31 22:30:00+00:00  ...                     0
## 2926 2018-07-31 23:00:00+00:00  ...                     0
## 2927 2018-07-31 23:30:00+00:00  ...                     0
## 
## [2928 rows x 36 columns]>

termlist = [f.split('.') for f in flux['NIWO'].keys()]
term = [t for sublist in termlist for t in sublist]
objDesc = flux['objDesc']
objDesc[objDesc['Object'].isin(term)]

##                Object                                        Description
## 21        angZaxsErth                                     Wind direction
## 83               data                             Represents data fields
## 90           distReso      Resolution of vector of distances (units = m)
## 91         distXaxs90  Distance of 90% contriubtion of footprint alig...
## 92        distXaxsMax  Max distance of footprint contriubtion aligned...
## 93         distYaxs90  Distance of 90% contriubtion of footprint alig...
## 94        distZaxsAbl  Atmospheric boundary layer height used in flux...
## 99   distZaxsMeasDisp      Measurement height minus displacement height 
## 100       distZaxsRgh  Roughness length used in flux footprint calcul...
## 124              flux                              General term for flux
## 125           fluxCo2      CO2 flux (can be turbulent, storage, or nsae)
## 127           fluxH2o  H2O flux (can be turbulent, storage, or nsae) ...
## 129          fluxMome  Momentum flux (can be turbulent, storage, or n...
## 131          fluxTemp  Temperature flux (can be turbulent, storage, o...
## 132              foot                                          Footprint
## 199              nsae  Net surface atmospher exchange (turbulent + st...
## 243            qfFinl  The final quality flag indicating if the data ...
## 273              qfqm  Quality flag and quality metrics, represents q...
## 450              stat                                         Statistics
## 452              stor                                            Storage
## 473           timeBgn       The beginning time of the aggregation period
## 474           timeEnd             The end time of the aggregation period
## 476              turb                                          Turbulent
## 504          veloFric            Friction velocity/ustar (units = m s-1)
## 511     veloYaxsHorSd  Standard deviation of horizontal wind speed in...
## 514     veloZaxsHorSd  Standard deviation of horizontal wind speed in...

flux['variables']

##    category    system variable     stat            units
## 0      data   fluxCo2     nsae  timeBgn               NA
## 1      data   fluxCo2     nsae  timeEnd               NA
## 2      data   fluxCo2     nsae     flux  umolCo2 m-2 s-1
## 3      data   fluxCo2     stor  timeBgn               NA
## 4      data   fluxCo2     stor  timeEnd               NA
## ..      ...       ...      ...      ...              ...
## 73     qfqm  fluxTemp     turb  timeEnd               NA
## 74     qfqm  fluxTemp     turb   qfFinl               NA
## 75     qfqm      foot     turb  timeBgn               NA
## 76     qfqm      foot     turb  timeEnd               NA
## 77     qfqm      foot     turb   qfFinl               NA
## 
## [78 rows x 5 columns]

plot(flux$NIWO$data.fluxCo2.nsae.flux~flux$NIWO$timeBgn, 
     pch=".", xlab="Date", ylab="CO2 flux",
     ylim=c(-20,20))

plot(flux$NIWO$data.fluxCo2.nsae.flux~flux$NIWO$timeBgn, 
     pch=20, xlab="Date", ylab="CO2 flux",
     xlim=c(as.POSIXct("2018-07-07", tz="GMT"),
            as.POSIXct("2018-07-09", tz="GMT")),
    ylim=c(-20,20), xaxt="n")
axis.POSIXct(1, x=flux$NIWO$timeBgn, 
             format="%Y-%m-%d %H:%M:%S")

plt.scatter(flux['NIWO'].timeBgn, 
            flux['NIWO']['data.fluxCo2.nsae.flux'],
           marker='.')
plt.ylim((-20,20))

## (-20.0, 20.0)

plt.show()

plt.scatter(flux['NIWO'].timeBgn, 
            flux['NIWO']['data.fluxCo2.nsae.flux'],
           marker='.')
plt.ylim((-20,20))

## (-20.0, 20.0)

plt.xlim((datetime.strptime('2018-07-07 00:00:00', '%Y-%m-%d %H:%M:%S'),
          datetime.strptime('2018-07-09 00:00:00', '%Y-%m-%d %H:%M:%S')))

## (np.float64(17719.0), np.float64(17721.0))

plt.show()

pr <- loadByProduct("DP1.00024.001", site="NIWO", 
                    timeIndex=30, package="basic", 
                    startdate="2018-06", enddate="2018-07",
                    check.size=F,
                    token=token)

pr = nu.load_by_product("DP1.00024.001", site="NIWO", 
                    timeindex=30, package="basic", 
                    startdate="2018-06", 
                    enddate="2018-07",
                    check_size=False,
                    token=token)

pr.top <- pr$PARPAR_30min[which(pr$PARPAR_30min$verticalPosition==
                                max(pr$PARPAR_30min$verticalPosition)),]

pr.top$timeBgn <- pr.top$startDateTime

fx.pr <- merge(pr.top, flux$NIWO, by="timeBgn")

plot(fx.pr$data.fluxCo2.nsae.flux~fx.pr$PARMean,
     pch=".", ylim=c(-20,20),
     xlab="PAR", ylab="CO2 flux")

par30 = pr['PARPAR_30min']
par30['timeBgn'] = par30['startDateTime']
prtop = par30[par30['verticalPosition']=='040']

fxpr = pd.merge(flux['NIWO'], prtop, how='outer', on='timeBgn')

plt.scatter(fxpr.PARMean, 
            fxpr['data.fluxCo2.nsae.flux'],
           marker='.')
plt.ylim((-20,20))

## (-20.0, 20.0)

plt.show()

prof <- stackEddy(filepath=paste0(getwd(), "/filesToStack00200"),
                 level="dp03")

head(prof$NIWO[,1:7])

##               timeBgn             timeEnd data.co2Stor.rateRtioMoleDryCo2.0.1 m
## 1 2018-06-01 00:00:00 2018-06-01 00:29:59                         -2.681938e-04
## 2 2018-06-01 00:30:00 2018-06-01 00:59:59                          1.912094e-05
## 3 2018-06-01 01:00:00 2018-06-01 01:29:59                          3.709964e-04
## 4 2018-06-01 01:30:00 2018-06-01 01:59:59                          2.374054e-03
## 5 2018-06-01 02:00:00 2018-06-01 02:29:59                          3.397046e-03
## 6 2018-06-01 02:30:00 2018-06-01 02:59:59                          3.156952e-03
##   data.co2Stor.rateRtioMoleDryCo2.0.2 m data.co2Stor.rateRtioMoleDryCo2.0.3 m
## 1                         -2.681938e-04                         -2.681938e-04
## 2                          1.912094e-05                          1.912094e-05
## 3                          3.709964e-04                          3.709964e-04
## 4                          2.374054e-03                          2.374054e-03
## 5                          3.397046e-03                          3.397046e-03
## 6                          3.156952e-03                          3.156952e-03
##   data.co2Stor.rateRtioMoleDryCo2.0.4 m data.co2Stor.rateRtioMoleDryCo2.0.5 m
## 1                         -2.681938e-04                         -2.681938e-04
## 2                          1.912094e-05                          1.912094e-05
## 3                          3.709964e-04                          3.709964e-04
## 4                          2.374054e-03                          2.374054e-03
## 5                          3.397046e-03                          3.397046e-03
## 6                          3.156952e-03                          3.156952e-03

prof = nu.stack_eddy(filepath=os.getcwd() + "/filesToStack00200",
                 level="dp03")

prof['NIWO'].head()

prof.l2 <- stackEddy(filepath=paste0(getwd(), "/filesToStack00200"),
                 level="dp02")

head(prof.l2$HARV)

##   horizontalPosition verticalPosition             timeBgn             timeEnd
## 1                000              010 2018-06-01 00:00:00 2018-06-01 00:29:59
## 2                000              010 2018-06-01 00:30:00 2018-06-01 00:59:59
## 3                000              010 2018-06-01 01:00:00 2018-06-01 01:29:59
## 4                000              010 2018-06-01 01:30:00 2018-06-01 01:59:59
## 5                000              010 2018-06-01 02:00:00 2018-06-01 02:29:59
## 6                000              010 2018-06-01 02:30:00 2018-06-01 02:59:59
##   data.co2Stor.rateRtioMoleDryCo2.mean data.h2oStor.rateRtioMoleDryH2o.mean
## 1                                  NaN                                  NaN
## 2                          0.002666576                                  NaN
## 3                         -0.011224223                                  NaN
## 4                          0.006133056                                  NaN
## 5                         -0.019554655                                  NaN
## 6                         -0.007855632                                  NaN
##   data.tempStor.rateTemp.mean qfqm.co2Stor.rateRtioMoleDryCo2.qfFinl
## 1                2.583333e-05                                      1
## 2               -2.008056e-04                                      1
## 3               -1.901111e-04                                      1
## 4               -7.419444e-05                                      1
## 5               -1.537083e-04                                      1
## 6               -1.874861e-04                                      1
##   qfqm.h2oStor.rateRtioMoleDryH2o.qfFinl qfqm.tempStor.rateTemp.qfFinl
## 1                                      1                             0
## 2                                      1                             0
## 3                                      1                             0
## 4                                      1                             0
## 5                                      1                             0
## 6                                      1                             0

profl2 = nu.stack_eddy(filepath=os.getcwd() + "/filesToStack00200",
                 level="dp02")

profl2['HARV'].head()

##   horizontalPosition  ... qfqm.tempStor.rateTemp.qfFinl
## 0                000  ...                             0
## 1                000  ...                             0
## 2                000  ...                             0
## 3                000  ...                             0
## 4                000  ...                             0
## 
## [5 rows x 10 columns]

vars <- getVarsEddy(paste0(getwd(),
                           "/filesToStack00200/NEON.D01.HARV.DP4.00200.001.nsae.2018-07.basic.20260115T064348Z.h5"))
head(vars)

##    site level category system hor ver tmi subsys       name       otype
## 5  HARV  dp01     data   amrs 000 060 01m   <NA> angNedXaxs H5I_DATASET
## 6  HARV  dp01     data   amrs 000 060 01m   <NA> angNedYaxs H5I_DATASET
## 7  HARV  dp01     data   amrs 000 060 01m   <NA> angNedZaxs H5I_DATASET
## 9  HARV  dp01     data   amrs 000 060 30m   <NA> angNedXaxs H5I_DATASET
## 10 HARV  dp01     data   amrs 000 060 30m   <NA> angNedYaxs H5I_DATASET
## 11 HARV  dp01     data   amrs 000 060 30m   <NA> angNedZaxs H5I_DATASET
##      dclass   dim  oth
## 5  COMPOUND 44640 <NA>
## 6  COMPOUND 44640 <NA>
## 7  COMPOUND 44640 <NA>
## 9  COMPOUND  1488 <NA>
## 10 COMPOUND  1488 <NA>
## 11 COMPOUND  1488 <NA>

iso <- stackEddy(filepath=paste0(getwd(), "/filesToStack00200"),
               level="dp01", var=c("rtioMoleDryCo2","rtioMoleDryH2o",
                                   "dlta13CCo2","dlta18OH2o"), avg=30)

head(iso$HARV[,1:8])

##   horizontalPosition verticalPosition             timeBgn             timeEnd
## 1                000              010 2018-06-01 00:00:00 2018-06-01 00:29:59
## 2                000              010 2018-06-01 00:30:00 2018-06-01 00:59:59
## 3                000              010 2018-06-01 01:00:00 2018-06-01 01:29:59
## 4                000              010 2018-06-01 01:30:00 2018-06-01 01:59:59
## 5                000              010 2018-06-01 02:00:00 2018-06-01 02:29:59
## 6                000              010 2018-06-01 02:30:00 2018-06-01 02:59:59
##   data.co2Stor.rtioMoleDryCo2.mean data.co2Stor.rtioMoleDryCo2.min
## 1                         509.3375                        451.4786
## 2                         502.2736                        463.5470
## 3                         521.6139                        442.8649
## 4                         469.6317                        432.6588
## 5                         484.7725                        436.2842
## 6                         476.8554                        443.7055
##   data.co2Stor.rtioMoleDryCo2.max data.co2Stor.rtioMoleDryCo2.vari
## 1                        579.3518                         845.0795
## 2                        533.6622                         161.3652
## 3                        563.0518                         547.9924
## 4                        508.7463                         396.8379
## 5                        537.4641                         662.9449
## 6                        515.6598                         246.6969

iso = nu.stack_eddy(filepath=os.getcwd() + "/filesToStack00200",
                 level="dp01", avg=30, 
                    var=["rtioMoleDryCo2","rtioMoleDryH2o",
                          "dlta13CCo2","dlta18OH2o"])

iso['HARV'].head()

##   horizontalPosition  ... ucrt.h2oTurb.rtioMoleDryH2o.se
## 0                000  ...                            NaN
## 1                000  ...                            NaN
## 2                000  ...                            NaN
## 3                000  ...                            NaN
## 4                000  ...                            NaN
## 
## [5 rows x 85 columns]

iso.d <- iso$HARV[grep("2018-06-25", iso$HARV$timeBgn, fixed=T),]
iso.d <- iso.d[-which(is.na(as.numeric(iso.d$verticalPosition))),]

g <- ggplot(iso.d, aes(y=verticalPosition)) + 
  geom_path(aes(x=data.co2Stor.rtioMoleDryCo2.mean, 
                group=timeBgn, col=timeBgn)) + 
  theme(legend.position="none") + 
  xlab("CO2") + ylab("Tower level")
g

g <- ggplot(iso.d, aes(y=verticalPosition)) + 
  geom_path(aes(x=data.isoCo2.dlta13CCo2.mean, 
                group=timeBgn, col=timeBgn)) + 
  theme(legend.position="none") + 
  xlab("d13C") + ylab("Tower level")
g

isoH = iso['HARV']
isoH6 = isoH[(isoH['timeBgn'] >= '2018-06-25') & (isoH['timeBgn'] < '2018-06-26')]
isod = isoH6[isoH6['horizontalPosition'] == '000']

fig, ax = plt.subplots()
groups = isod["timeBgn"].unique()
cmap = plt.cm.viridis
colors = cmap(np.linspace(0, 1, len(groups)))
for color, group_name in zip(colors, groups):
group = isod[isod["timeBgn"] == group_name]
ax.plot(group['data.co2Stor.rtioMoleDryCo2.mean'],
group['verticalPosition'],
color=color)
ax.set_xlabel("CO2")
ax.set_ylabel("Tower level")
plt.show()

fig, ax = plt.subplots()
for color, group_name in zip(colors, groups):
group = isod[isod["timeBgn"] == group_name]
ax.plot(group["data.isoCo2.dlta13CCo2.mean"],
group["verticalPosition"],
color=color)
ax.set_xlabel("13C")
ax.set_ylabel("Tower level")
plt.show()

An extremely common task for remote sensing researchers is to connect remote sensing data with data on the ground. This might be a research question in itself - what reflectance wavelengths are most closely correlated with a particular ground feature? - or a ground-truthing for extrapolation and prediction, or development or testing of a model.

This tutorial explores the relationship between a Lidar-derived canopy height model and tree heights measured from the ground, because these two datasets provide a straightforward introduction to thinking about how to relate different sources of data. They are a good exemplar for the two major components of connecting airborne and ground data:

1. The mechanics: Finding the remotely sensed pixels that correspond to a given ground measurement
2. The science: Considering biases in each data source, and the ways the measurements might differ even if neither is incorrect

We will explore these issues between these two datasets, and you can use what you learn here as a roadmap for making similar comparisons between different datasets.

The two NEON data products that estimate tree height:

• DP3.30015.001, Ecosystem structure, aka Canopy Height Model (CHM)
• DP1.10098.001, Vegetation structure

The CHM data are derived from the Lidar point cloud data collected by the remote sensing platform. The vegetation structure data are collected by by field staff on the ground.

We will be using data from the Wind River Experimental Forest (WREF) NEON field site located in Washington state. The predominant vegetation there is tall evergreen conifers.

install.packages("neonUtilities")
install.packages("neonOS")
install.packages("terra")
install.packages("dplyr")
install.packages("devtools")
devtools::install_github("NEONScience/NEON-geolocation/geoNEON")

library(terra)
library(dplyr)
library(neonUtilities)
library(neonOS)
library(geoNEON)
token <- Sys.getenv("NEON_TOKEN")
wd <- "~/data"

pip install neonutilities
pip install geopandas
pip install rasterio
pip install rioxarray
pip install rasterstats
pip install python-dotenv

import neonutilities as nu
import pandas as pd
import numpy as np
import rasterstats as rs
import geopandas as gpd
import rioxarray as rxr
import matplotlib.pyplot as plt
import matplotlib.collections
import rasterio
from rasterio import sample
from rasterio.enums import Resampling
import requests
import time
import os
import dotenv
dotenv.load_dotenv()

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

veglist <- loadByProduct(dpID="DP1.10098.001", 
                         site="WREF", 
                         package="basic", 
                         release="RELEASE-2024",
                         check.size = FALSE,
                         token=token)

veglist = nu.load_by_product(dpid="DP1.10098.001", 
                         site="WREF", 
                         package="basic", 
                         release="RELEASE-2024",
                         check_size = False,
                         token=token)

vegmap <- getLocTOS(veglist$vst_mappingandtagging, 
                          "vst_mappingandtagging")

## Please note locations will be calculated only for mapped woody individuals. To find subplot locations for unmapped individuals, use this function with the vst_apparentindividual, vst_non-woody, and/or vst_shrubgroup tables.

vegmapall = veglist["vst_mappingandtagging"]
vegmap = vegmapall.loc[vegmapall["pointID"] != ""]
vegmap = vegmap.reindex()
vegmap["points"] = vegmap["namedLocation"] + "." + vegmap["pointID"]
veg_points = list(set(list(vegmap["points"])))
easting = []
northing = []
coord_uncertainty = []
elev_uncertainty = []
drop_points = []

for i in veg_points:
    time.sleep(1)
    vres = requests.get("https://data.neonscience.org/api/v0/locations/"+i)
    vres_json = vres.json()
    if vres_json["data"] is not None:
        easting.append(vres_json["data"]["locationUtmEasting"])
        northing.append(vres_json["data"]["locationUtmNorthing"])
        props = pd.DataFrame.from_dict(vres_json["data"]["locationProperties"])
        cu = props.loc[props["locationPropertyName"] == "Value for Coordinate uncertainty"]["locationPropertyValue"]
        coord_uncertainty.append(cu[cu.index[0]])
        eu = props.loc[props["locationPropertyName"] == "Value for Elevation uncertainty"]["locationPropertyValue"]
        elev_uncertainty.append(eu[eu.index[0]])
    else:
        drop_points.append(i)

veg_points_clean = [v for v in veg_points if v not in drop_points]

ptdct = dict(points=veg_points_clean, 
             easting=easting, 
             northing=northing,
             coordinateUncertainty=coord_uncertainty,
             elevationUncertainty=elev_uncertainty)
ptfrm = pd.DataFrame.from_dict(ptdct)
ptfrm.set_index("points", inplace=True)

vegmap = vegmap.join(ptfrm, 
                     on="points", 
                     how="inner")

vegmap["adjEasting"] = (vegmap["easting"]
                        + vegmap["stemDistance"]
                        * np.sin(vegmap["stemAzimuth"]
                                   * np.pi / 180))
vegmap["adjNorthing"] = (vegmap["northing"]
+ vegmap["stemDistance"]
* np.cos(vegmap["stemAzimuth"]
* np.pi / 180))

vegmap["adjCoordinateUncertainty"] = vegmap["coordinateUncertainty"] + 0.6
vegmap["adjElevationUncertainty"] = vegmap["elevationUncertainty"] + 1

veg <- joinTableNEON(veglist$vst_apparentindividual, 
                     vegmap, 
                     name1="vst_apparentindividual",
                     name2="vst_mappingandtagging")

veglist["vst_apparentindividual"].set_index("individualID", inplace=True)
veg = vegmap.join(veglist["vst_apparentindividual"],
                  on="individualID",
                  how="inner",
                  lsuffix="_MAT",
                  rsuffix="_AI")

veg2017 <- veg |>
  filter(eventID.y=="vst_WREF_2017")
veg2017.075 <- veg2017 |>
filter(plotID=="WREF_075")
symbols(veg2017.075$adjEasting,
veg2017.075$adjNorthing,
circles=veg2017.075$stemDiameter/100/2,
xlab="Easting",
ylab="Northing",
inches=F)

symbols(veg2017.075$adjEasting, 
        veg2017.075$adjNorthing, 
        circles=veg2017.075$stemDiameter/100/2, 
        xlab="Easting", 
        ylab="Northing", 
        inches=F)
symbols(veg2017.075$adjEasting, 
        veg2017.075$adjNorthing, 
        circles=veg2017.075$adjCoordinateUncertainty, 
        inches=F, add=T, 
        fg="lightblue")

veg2017 = veg.loc[veg["eventID_AI"]=="vst_WREF_2017"]
veg75 = veg2017.loc[veg2017["plotID_AI"]=="WREF_075"]
fig, ax = plt.subplots()
xy = np.array(tuple(zip(veg75.adjEasting, veg75.adjNorthing)))
srad = veg75.stemDiameter/100/2
patches = [plt.Circle(center, size) for center, size in zip(xy, srad)]
coll = matplotlib.collections.PatchCollection(patches, facecolors="white", edgecolors="black")
ax.add_collection(coll)
ax.margins(0.1)
plt.show()

fig, ax = plt.subplots()
sunc = veg75.adjCoordinateUncertainty
patchunc = [plt.Circle(center, size) for center, size in zip(xy, sunc)]
coll = matplotlib.collections.PatchCollection(patches, facecolors="None", edgecolors="black")
collunc = matplotlib.collections.PatchCollection(patchunc, facecolors="None", edgecolors="lightblue")
ax.add_collection(coll)
ax.add_collection(collunc)
ax.margins(0.1)
plt.show()

byTileAOP(dpID="DP3.30015.001", site="WREF", year=2017, 
          easting=veg2017.075$adjEasting, 
          northing=veg2017.075$adjNorthing,
          check.size=FALSE, savepath=wd,
          token=token)
chm <- rast(paste0(wd, "/DP3.30015.001/neon-aop-products/2017/FullSite/D16/2017_WREF_1/L3/DiscreteLidar/CanopyHeightModelGtif/NEON_D16_WREF_DP3_580000_5075000_CHM.tif"))

plot(chm, col=topo.colors(5))

nu.by_tile_aop(dpid="DP3.30015.001", site="WREF", year="2017", 
          easting=list(veg75.adjEasting), 
          northing=list(veg75.adjNorthing),
          check_size=False, savepath=os.getcwd(),
          token=token)

chm = rasterio.open(os.getcwd() + "/DP3.30015.001/neon-aop-products/2017/FullSite/D16/2017_WREF_1/L3/DiscreteLidar/CanopyHeightModelGtif/NEON_D16_WREF_DP3_580000_5075000_CHM.tif")
chmx = rxr.open_rasterio(os.getcwd() + "/DP3.30015.001/neon-aop-products/2017/FullSite/D16/2017_WREF_1/L3/DiscreteLidar/CanopyHeightModelGtif/NEON_D16_WREF_DP3_580000_5075000_CHM.tif").squeeze()

plt.imshow(chm.read(1))
plt.show()

vegsub <- veg2017 |>
  filter(adjEasting >= ext(chm)[1] &
         adjEasting <= ext(chm)[2] &
         adjNorthing >= ext(chm)[3] & 
         adjNorthing <= ext(chm)[4])

vegsub = veg2017.loc[(veg2017["adjEasting"] >= chm.bounds[0]) &
                 (veg2017["adjEasting"] <= chm.bounds[1]) &
                 (veg2017["adjNorthing"] >= chm.bounds[2]) & 
                 (veg2017["adjNorthing"] <= chm.bounds[3])]
vegsub = vegsub.reset_index(drop=True)

valCHM <- extract(chm, 
                  cbind(vegsub$adjEasting,
                  vegsub$adjNorthing))
plot(valCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegsub$height, pch=20, xlab="Height",
ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(valCHM$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegsub$height, use="complete")

## [1] 0.3910897

valCHM = list(sample.sample_gen(chm, 
                                tuple(zip(vegsub["adjEasting"], 
                                          vegsub["adjNorthing"])),
                                masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegsub.height, valCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

CHMlist = np.array([c.tolist()[0] for c in valCHM])
idx = np.intersect1d(np.where(np.isfinite(vegsub.height)), 
                     np.where(CHMlist != None))
np.corrcoef(vegsub.height[idx], list(CHMlist[idx]))[0,1]

## np.float64(0.3913279146475449)

valCHMbuff <- extract(chm, 
                  buffer(vect(cbind(vegsub$adjEasting,
                  vegsub$adjNorthing)),
                  width=vegsub$adjCoordinateUncertainty),
                  fun=max)
plot(valCHMbuff$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegsub$height, pch=20, xlab="Height",
ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(valCHMbuff$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegsub$height, use="complete")

## [1] 0.3846956

vegloc = vegsub[["individualID","adjEasting","adjNorthing","adjCoordinateUncertainty"]]
v = vegloc.rename(columns={"individualID": "indID", "adjEasting": "easting",
               "adjNorthing": "northing", "adjCoordinateUncertainty": "coordUnc"},
               inplace=False)
gdf = gpd.GeoDataFrame(
v, geometry=gpd.points_from_xy(v.easting, v.northing))
gdf["geometry"] = gdf["geometry"].buffer(distance=gdf["coordUnc"])
gdf.to_file(os.getcwd() + "/trees_with_buffer.shp")
chm_height = rs.zonal_stats(os.getcwd() + "/trees_with_buffer.shp", chmx.values,
affine=chmx.rio.transform(),
nodata=-9999, stats="max")
valCHMbuff = [h["max"] for h in chm_height]

fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegsub.height, valCHMbuff, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

CHMbufflist = np.array(valCHMbuff)
idx = np.intersect1d(np.where(np.isfinite(vegsub.height)), 
                     np.where(CHMbufflist != None))
np.corrcoef(vegsub.height[idx], list(CHMbufflist[idx]))[0,1]

## np.float64(0.3848610594340662)

easting10 <- 10*floor(vegsub$adjEasting/10)
northing10 <- 10*floor(vegsub$adjNorthing/10)
vegsub <- cbind(vegsub, easting10, northing10)

vegbin <- stats::aggregate(vegsub, 
                           by=list(vegsub$easting10, 
                                   vegsub$northing10), 
                           FUN=max)

CHM10 <- terra::aggregate(chm, fact=10, fun=max)
plot(CHM10, col=topo.colors(5))

vegbin$easting10 <- vegbin$easting10 + 5
vegbin$northing10 <- vegbin$northing10 + 5
binCHM <- extract(CHM10, cbind(vegbin$easting10, 
                               vegbin$northing10))
plot(binCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
       vegbin$height, pch=20, 
     xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(binCHM$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegbin$height, use="complete")

## [1] 0.3593519

vegsub["easting10"] = 10*np.floor(vegsub.adjEasting/10)
vegsub["northing10"] = 10*np.floor(vegsub.adjNorthing/10)
vegsubloc = vegsub[["height","easting10","northing10"]]

vegbin = vegsubloc.groupby(["easting10", "northing10"]).max().add_suffix('_max').reset_index()

target_res = (10, 10)
with rasterio.open(os.getcwd() + "/DP3.30015.001/neon-aop-products/2017/FullSite/D16/2017_WREF_1/L3/DiscreteLidar/CanopyHeightModelGtif/NEON_D16_WREF_DP3_580000_5075000_CHM.tif") as src:
data, transform = rasterio.warp.reproject(source=src.read(),
src_transform=src.transform,
src_crs=src.crs,
dst_crs=src.crs,
dst_nodata=src.nodata,
dst_resolution=target_res,
resampling=Resampling.max)
profile = src.profile
profile.update(transform=transform, driver='GTiff',
height=data.shape[1], width=data.shape[2])
 with rasterio.open(os.getcwd() + '/CHM_10m.tif', 'w', **profile) as dst:
                dst.write(data)

chm10 = rasterio.open(os.getcwd() + '/CHM_10m.tif')

plt.imshow(chm10.read(1))
plt.show()

valCHM10 = list(sample.sample_gen(chm10, tuple(zip(vegbin["easting10"]+5,
                                                   vegbin["northing10"]+5)),
                                                   masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegbin.height_max, valCHM10, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

CHM10list = np.array([c.tolist()[0] for c in valCHM10])
idx = np.intersect1d(np.where(np.isfinite(vegbin.height_max)), 
                     np.where(CHM10list != None))
np.corrcoef(vegbin.height_max[idx], list(CHM10list[idx]))[0,1]

## np.float64(0.35904092842414137)

vegsub <- vegsub |>
  arrange(desc(height))

vegfil <- vegsub
for(i in 1:nrow(vegsub)) {
    if(is.na(vegfil$height[i]))
        next
    dist <- sqrt((vegsub$adjEasting[i]-vegsub$adjEasting)^2 + 
                (vegsub$adjNorthing[i]-vegsub$adjNorthing)^2)
    vegfil$height[which(dist<0.3*vegsub$height[i] & 
                        vegsub$height<vegsub$height[i])] <- NA
}
vegfil <- vegfil |>
filter(!is.na(height))

filterCHM <- extract(chm, 
                     cbind(vegfil$adjEasting, 
                           vegfil$adjNorthing))
plot(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
       vegfil$height, pch=20, 
     xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM,
    vegfil$height)

## [1] 0.8170048

vegfil = vegsub.sort_values(by="height", ascending=False, ignore_index=True)

height = vegfil.height.reset_index()
for i in vegfil.index:
    if height.height[i] is None:
        pass
    else:
        dist = np.sqrt(np.square(vegfil.adjEasting[i]-vegfil.adjEasting) + 
                       np.square(vegfil.adjNorthing[i]-vegfil.adjNorthing))
        idx = vegfil.index[(vegfil.height<height.height[i]) & (dist<0.3*height.height[i])]
        height.loc[idx, "height"] = None

filterCHM = list(sample.sample_gen(chm, tuple(zip(vegfil["adjEasting"],
                                                  vegfil["adjNorthing"])),
                                                  masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(height.height, filterCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

filCHMlist = np.array([c.tolist()[0] for c in filterCHM])
idx = np.intersect1d(np.where(np.isfinite(height.height)), 
                     np.where(filCHMlist != None))
np.corrcoef(height.height[idx], list(filCHMlist[idx]))[0,1]

## np.float64(0.8393670222898226)

vegfil <- vegfil |>
  filter(plantStatus=="Live")
filterCHM <- extract(chm,
cbind(vegfil$adjEasting,
vegfil$adjNorthing))
plot(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegfil$height, pch=20,
xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM,
    vegfil$height)

## [1] 0.9113093

idx = vegfil.index[vegfil.plantStatus!="Live"]
height.loc[idx, "height"] = None
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(height.height, filterCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

idx = np.intersect1d(np.where(np.isfinite(height.height)), 
                     np.where(filCHMlist != None))
np.corrcoef(height.height[idx], list(filCHMlist[idx]))[0,1]

## np.float64(0.9113648148226839)