Merging AOP L3 Tiles in R into Full-Site Rasters

Most AOP L3 data are provided in 1km x 1km tiled geotiff rasters. A single site may contain hundreds of separate raster data files for each data product. This tutorial shows you how to merge these tiles into a single raster, covering the full site, which enables more efficient regional-scale analysis.

In this tutorial, you will import the aop_mosaic_rasters.R script and use the function makeFullSiteMosaics to download AOP L3 tiled raster (geotiff) data sets for a given data product, site, and year. This function uses the neonUtilities R package byFileAOP function to download the data. Once downloaded, the tiles are merged into a single raster covering the full site (or whatever coverage was obtained that year). Finally, this full-site raster is saved as a geotiff and cloud-optimized geotiff (COG) file in an output folder specified in the function.

Once you have generated full-site geotiffs, it is much simpler to visualize and conduct further analysis of the data in your preferred geospatial application, eg. ArcGIS, QGIS, or Google Earth Engine.

This tutorial illustrates the function for the NEON aquatic site MCRA (McRae Creek) in Domain 16, the Pacific Northwest. This site was chosen for demonstration because it is one of the smallest AOP sites, and therefore is quicker to download. Let's get started!

First, clone the git repository locally and set the working directory to where you cloned the data skills repository. You will need to modify the file path below to reflect where you have cloned the repository, or saved the aop_merge_raster_functions.R script.

# change the wd depending on your local environment
wd <- file.path(path.expand("~"),"GitHubRepos","NEON-Data-Skills","tutorials","R","AOP","AOP-L3-rasters") 
setwd(wd)

First let's use source to pull in the aop_merge_raster_functions.R script, which is located in the working directory. For more details on what source does, please refer to the CU Earth Lab Tutorial on Source Functions in R

source("./aop_merge_raster_functions.R")

Let's take a look a the documentation for the function makeFullSiteMosaics

help(makeFullSiteMosaics)

You should see the function documentation pop up in the RStudio Help window, similar to below:

makeFullSiteMosaics

{TempPackage} R Documentation

Download all AOP files for a given site, year, and L3 product, mosaic the files, and save the full site mosaic to a tiff and cloud-optimized geotiff.

Description

This function 1) Runs the neonUtilities byFileAOP function to download NEON AOP data by site, year, and product (see byFileAOP documention for additional details). 2) mosaics the data product tiles into a full-site mosaic, as well as the associated error tifs, where applicable, and 3) saves the full site mosaics to a tif and cloud-optimized geotiff.

Usage

makeFullSiteMosaics( dpID, year, siteCode, dataRootDir, outFileDir, apiToken = NULL)

Arguments

dpID

The identifier of the AOP data product to download, in the form DP3.PRNUM.REV, e.g. DP3.30011.001. This works for all AOP L3 rasters except L3 reflectance. If an invalid data product ID is provided, the code will show an error message and display the valid dpIDs.

year

The four-digit year to search for data.

siteCode

The four-letter code of a single NEON site, e.g. 'MCRA'.

dataRootDir

The file path to download data to.

outFileDir

The file path where the full-site mosaic geotiffs and cloud-optimized geotiffs are saved.

apiToken

User specific API token (generated within neon.datascience user accounts). If not provided, no API token is used.

TIP: We recommend using an API token when downloading NEON data programmatically. The function will work without a token, if you leave it out (as described in the documentation), but it is best to get in the practice of using the token! Here I source the token from a file, where the file consists of a single line of code defining the token (called NEON_TOKEN). To set up an API token, please refer to the R tutorial on generating tokens, linked in the requirements section at the beginning of this lesson.

NEON_TOKEN <- "MY TOKEN"

We can read in the token using the source function as follows. This assumes the token is saved in the working directory, but you can also set the path to the token explicitly if you've saved it elsewhere.

source(paste0(wd,"/neon_token.R"))

Now that we have a general idea of how this function works, from the documentation, let's go ahead and run it, including the imported token. For this example, I set the download folder to 'c:/neon/data' and the output folder to 'c:/neon/outputs/2021_MCRA/CHM/'. Modify these paths as desired according to your project structure.

Please heed the following warnings before you run the code:

WARNING: This function is currently set so that it does not check the file size before downloading. You can change this, if desired, by either removing the check.size parameter in the the makeFullSiteMosaics function, or by changing the value for that setting to True (check.size=T). This will pause the function and prompt you to type y/n (yes or no) to continue with the download. We recommend changing this parameter if you have limited storage space on your computer, but note that doing so will require a manual input to complete the execution of the function.

WARNING: We recommend extending the timeout value when downloading large AOP sites so the connection doesn't stall out before the files finishe downloading. To change this setting, you can the R options, but note that this will change the computer's settings outside of this R environment. When you do set the timeout, be sure to change it back at the end. The default timeout is 60 seconds, we recommend changing it to 300 seconds (5 minutes).

timeout0 <- getOption('timeout')
print(timeout0)

## [1] 60

Set the timeout to 300 seconds (5 minutes):

options(timeout=300)
getOption('timeout')

## [1] 300

Set the folder to download the data to, and the output folder to save the final CHM rasters:

download_folder <- 'c:/neon/data'
chm_output_folder <- 'c:/neon/outputs/2021_MCRA/CHM/'


makeFullSiteMosaics('DP3.30015.001','2021','MCRA',download_folder,chm_output_folder,NEON_TOKEN)

Display the output files generated:

list.files(chm_output_folder)

## [1] "2021_MCRA_2_CHM.tif"     "2021_MCRA_2_CHM_COG.tif"

Now we can read in and plot the full-site CHM tif that we generated using the raster package:

MCRA_CHM <- raster(file.path(chm_output_folder,'2021_MCRA_2_CHM.tif'))
plot(MCRA_CHM,main="2021_MCRA_2 Canopy Height Model") # add title with main)

Now let's run the function to generate full-site mosaics for the AOP Veg Indices (DP3.30026.001). The function handles unzipping the zipped folders, and plotting each separate index tif and associated error tif. We will use the same download folder as before (note that the files will be downloaded to a different sub-directory), and specify a new output folder.

veg_indices_output_folder<-'c:/neon/outputs/2021_MCRA/VegIndices/'


makeFullSiteMosaics('DP3.30026.001','2021','MCRA',download_folder,veg_indices_output_folder,NEON_TOKEN)

Let's see the full-site Veg Indices files that were generated using list.files:

list.files(veg_indices_output_folder)

##  [1] "2021_MCRA_2_ARVI.tif"           "2021_MCRA_2_ARVI_COG.tif"       "2021_MCRA_2_ARVI_error.tif"    
##  [4] "2021_MCRA_2_ARVI_error_COG.tif" "2021_MCRA_2_EVI.tif"            "2021_MCRA_2_EVI_COG.tif"       
##  [7] "2021_MCRA_2_EVI_error.tif"      "2021_MCRA_2_EVI_error_COG.tif"  "2021_MCRA_2_NDVI.tif"          
## [10] "2021_MCRA_2_NDVI_COG.tif"       "2021_MCRA_2_NDVI_error.tif"     "2021_MCRA_2_NDVI_error_COG.tif"
## [13] "2021_MCRA_2_PRI.tif"            "2021_MCRA_2_PRI_COG.tif"        "2021_MCRA_2_PRI_error.tif"     
## [16] "2021_MCRA_2_PRI_error_COG.tif"  "2021_MCRA_2_SAVI.tif"           "2021_MCRA_2_SAVI_COG.tif"      
## [19] "2021_MCRA_2_SAVI_error.tif"     "2021_MCRA_2_SAVI_error_COG.tif"

Now we can read in the NDVI and NDVI error tifs:

MCRA_NDVI <- raster(file.path(veg_indices_output_folder,'2021_MCRA_2_NDVI.tif'))
MCRA_NDVI_error <- raster(file.path(veg_indices_output_folder,'2021_MCRA_2_NDVI_error.tif'))

And, finally, let's plot NDVI and NDVI error:

plot(MCRA_NDVI,main="2021_MCRA_2 NDVI")

plot(MCRA_NDVI_error,main="2021_MCRA_2 NDVI Error") 

Looks good! Go ahead and test out this function on different data products, years, or sites on your own. Note that larger sites will take more time to download the data, and may require significant memory resources on your computer.

Last but not least, change the timeout back to the original value:

options(timeout=timeout0)
getOption('timeout')

## [1] 60