Integrating airborne remote sensing and field campaigns for ecology and Earth system science

In recent years, the availability of airborne imaging spectroscopy (hyperspectral) data has expanded dramatically. The high spatial and spectral resolution of these data uniquely enable spatially explicit ecological studies including species mapping, assessment of drought mortality and foliar trait distributions. However, we have barely begun to unlock the potential of these data to use direct mapping of vegetation characteristics to infer subsurface properties of the critical zone. To assess their utility for Earth systems research, imaging spectroscopy data acquisitions require integration with large, coincident ground‐based datasets collected by experts in ecology and environmental and Earth science. Without coordinated, well‐planned field campaigns, potential knowledge leveraged from advanced airborne data collections could be lost. Despite the growing importance of this field, documented methods to couple such a wide variety of disciplines remain sparse. We coordinated the first National Ecological Observatory Network Airborne Observation Platform (AOP) survey performed outside of their core sites, which took place in the Upper East River watershed, Colorado. Extensive planning for sample tracking and organization allowed field and flight teams to update the ground‐based sampling strategy daily. This enabled collection of an extensive set of physical samples to support a wide range of ecological, microbiological, biogeochemical and hydrological studies. We present a framework for integrating airborne and field campaigns to obtain high‐quality data for foliar trait prediction and document an archive of coincident physical samples collected to support a systems approach to ecological research in the critical zone. This detailed methodological account provides an example of how a multi‐disciplinary and multi‐institutional team can coordinate to maximize knowledge gained from an airborne survey, an approach that could be extended to other studies. The coordination of imaging spectroscopy surveys with appropriately timed and extensive field surveys, along with high‐quality processing of these data, presents a unique opportunity to reveal new insights into the structure and dynamics of the critical zone. To our knowledge, this level of co‐aligned sampling has never been undertaken in tandem with AOP surveys and subsequent studies utilizing this archive will shed considerable light on the breadth of applications for which imaging spectroscopy data can be leveraged.


| INTRODUC TI ON
Spatially explicit and integrated measurements of ecological, biogeochemical and hydrological processes are increasingly important to test ecological theory and understand critical zone (CZ) evolution at higher spatial resolutions and increased levels of complexity. The critical zone concept, the integrated envelope of hydrobiogeochemical functioning from the top of vegetation canopies to the base of weathered bedrock, and its importance in our understanding of ecological processes, is a natural outgrowth of a long tradition of considering life and the life-supporting components of the planet as inextricable from one another. This integrated perspective has roots in Vernadsky's (1926) Biosphere, Tansley's (1935) ecosystem, Cole's (1958) ecosphere, Jenny's state factor framework for ecosystems (Amundson & Jenny, 1997), and Troll's geoecology or geoecosystems (Huggett, 1995;Troll, 1971). Current focus on exploring these interrelationships can be seen in a variety of efforts, from international Critical Zone Observatory (CZO) networks (Banwart et al., 2012;Brantley et al., 2017), to the grand challenge put forth by the U.S. National Research Council to explore the coevolution of landscapes and ecosystems (National Research Council, 2010;Troch et al., 2015), to efforts to incorporate hillslope hydrology into Earth system models (ESMs; Fan et al., 2019). Central to development of spatially explicit characterizations are remote sensing methods, paired with high-quality ground calibration data and supported by rapidly increasing computational capacity. These data can quantify surface and subsurface properties of the Earth system over extents and resolutions not possible to assess via ground-based sampling alone and have the potential to extend into areas that are remote and inaccessible. An understanding of these relationships and processes taking place across the CZ is essential for predicting ecosystem functioning, water resource availability and Earth system resilience to global change.
High fidelity imaging spectroscopy data, also known as hyperspectral data, have shown great promise for improving assessment of land surface characteristics over large areas and at fine spatial resolutions. NASA's Jet Propulsion Laboratory (JPL) initiated this effort with the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) programme and the AVIRIS Classic sensor, an airborne visible to shortwave infrared (VSWIR) imaging spectrometer that collects contiguous spectral bands from 350 to 2,500 nm using whiskbroom (across track) scanning (Vane et al., 1993) at 10 nm wavelength intervals. Building from AVIRIS Classic, a new generation of very high signal-to-noise, 5 nm pushbroom (along track) spectrometers including AVIRIS-NG run by the AVIRIS programme, the Global Airborne Observatory (GAO, formerly Carnegie Airborne Observatory; Asner et al., 2012) run by Arizona State University and the National Ecological Observatory Network's Airborne Observation Platform (NEON's AOP; Kampe, Johnson, Kuester, & Keller, 2009), are significantly expanding the quality and extent of available imaging spectroscopy data. In addition to a VSWIR spectrometer, both the GAO and AOP platforms include an integrated lidar sensor, which sample tracking and organization allowed field and flight teams to update the ground-based sampling strategy daily. This enabled collection of an extensive set of physical samples to support a wide range of ecological, microbiological, biogeochemical and hydrological studies.
3. We present a framework for integrating airborne and field campaigns to obtain high-quality data for foliar trait prediction and document an archive of coincident physical samples collected to support a systems approach to ecological research in the critical zone. This detailed methodological account provides an example of how a multi-disciplinary and multi-institutional team can coordinate to maximize knowledge gained from an airborne survey, an approach that could be extended to other studies.
4. The coordination of imaging spectroscopy surveys with appropriately timed and extensive field surveys, along with high-quality processing of these data, presents a unique opportunity to reveal new insights into the structure and dynamics of the critical zone. To our knowledge, this level of co-aligned sampling has never been undertaken in tandem with AOP surveys and subsequent studies utilizing this archive will shed considerable light on the breadth of applications for which imaging spectroscopy data can be leveraged.

K E Y W O R D S
airborne remote sensing, field surveys, foliar traits, hyperspectral imaging, imaging spectroscopy, metadata, NEON AOP, sample tracking collects co-aligned structural and topographic information. This sensor integration provides great insight into the structure, composition and function of ecosystems at meter-scale resolution across landscape and regional scales, currently not possible from spaceborne instruments. These integrated sensor packages have been utilized for, among other things, the mapping of aboveground carbon density, tree species, foliar traits, biodiversity and mineral composition of exposed surfaces (e.g. Clark et al., 2006;Colgan, Baldeck, Féret, & Asner, 2012;Féret & Asner, 2014;Smith et al., 2002;Wang et al., 2019).
While these data have primarily been used to elucidate patterns in the surface of the CZ, they may also provide a window into deeper CZ functioning. Because vegetation acts as a surface manifestation of less visible CZ components and processes, there is the potential to leverage the co-variability between vegetation and subsurface characteristics to map critical, but difficult to observe, soil properties. Early studies have indicated this methodology can provide insight into soil microbial functioning (Madritch et al., 2014;Soper et al., 2018), soil nutrient status Osborne et al., 2017) and the role of geomorphic processes in ecosystem organization (Chadwick & Asner, 2020).
Given that vegetation often obscures direct measurement of soil properties, linking land surface and belowground characteristics is a particularly promising area of research. However, many uncertainties need to be resolved to determine the potential of these relationships to assess landscape scale biogeochemical gradients and CZ geography.
Most ongoing research is focused on integration of field and airborne datasets to estimate land surface characteristics. The Spectranomics programme, associated with the GAO, has conducted extensive collections of foliar samples and species information which has supported both species and trait mapping in tropical forest systems around the world (Asner & Martin, 2009). NEON has expanded the range of ecosystem characteristics being sampled to include soil biogeochemistry and microbial samples within the footprint of AOP surveys. However, these additional biogeochemical and microbial samples are not always associated in time with the flights, or in close spatial proximity with the foliar samples, as required for directly linking above-and below-ground ecosystem properties (Hinckley et al., 2016;Kao et al., 2012).
High-quality, well-timed, and extensively organized ground sampling surveys that are targeted to calibrate, validate and extend site-specific models are necessary to assess the potential of imaging spectroscopy data to provide a window to the CZ. These efforts require close collaboration among airborne survey teams, field scientists and stakeholders in designing sampling schemes that ensure the greatest utility of the collected data. The scope of these sampling efforts exceeds those of many traditional field approaches, and thus require deeper levels of organization, collaboration and coordination between researchers across disciplines.
The methodology we present here provides the potential to transform an AOP dataset from a single investigator-based research project on surface distributions of vegetation characteristics, into interdisciplinary research on CZ and watershed science, supporting a number of research questions through the development of coordinated vegetation and soil archives, and geophysical datasets.
We developed an approach that leverages the NEON AOP to address questions regarding watershed-scale distributions of vegetation traits and communities, soil physical/chemical properties and microbial functional traits within the East River, Colorado, community watershed observatory (Hubbard et al., 2018). An important focus of our methodology is the collection of relevant field samples and data, co-aligned in time and space with the VSWIR and lidar data collected by the AOP over 330 km 2 across four watersheds in the Upper East River watershed, a headwater of the Colorado River.
This type of approach could be utilized to address a large range of research questions. Here, our central objective was to design sampling schemes that will allow us to explore three specific CZ focused questions: Can foliar composition and geomorphic position predict soil organic matter content? To what extent do subsurface structure and bedrock/soil properties, measured via geophysical imaging, dictate plant community distribution? Can imaging spectroscopy and landscape position predict the distribution of microbial functional traits throughout a watershed? Accordingly, we designed contemporaneous field surveys to collect (a) plant species and trait data and (b) a physical vegetation tissue and soil sample archive, which will be used to gain a vertically integrated, critical zone perspective.
In this manuscript, we provide a framework for others interested in pursuing integrated field-airborne research projects with the AOP or other airborne sensor platforms. While the details of these airborne and field collections are specific to this location and project goals, the approach is broadly applicable across a variety of sites with appropriate modifications to address different scientific questions. We detail the process of planning for an integrated airborne and field-based survey, highlight considerations for VSWIR data processing, report on field survey outcomes and provide a case study examining the generation of foliar trait maps. In the main text, we highlight the overall procedure and case study. In Supporting Information, we present detailed planning documents, code, input data and details of statistical modelling to demonstrate the process from planning to products. Accompanying the methodological descriptions, we also provide publicly available data products for individuals and teams.
Given the promise of this technology and the scale of investment by government agencies to support these instruments, we hope to promote the continued utilization of these methods by a broader scientific community. As we expand our knowledge of the CZ across sites through the National Science Foundation's Critical Zone Collaboration Network and gain a greatly extended imaging spectroscopy dataset through the development and deployment of NASA's Surface Biology and Geology (SBG) mission (National Academies of Sciences Engineering Medicine, 2018), studies exploring how imaging spectroscopy data reveal CZ structure and geography are essential to provide insight into ecosystem ecology and Earth system functioning.

| Study region
The data collections presented here were completed in the summer of 2018 across four watersheds in the Upper East River basin of Colorado ( Figure 1). The study region is home to the Rocky Mountain Biological Laboratory (RMBL) and the Department of Energy's (DOE) Watershed Function Scientific Focus Area (WF SFA) community watershed observatory and is described in detail by Hubbard et al. (2018). Briefly, the elevation range spans 2,800-4,000 m above sea level and the area is among the coldest in the Rocky Mountain region (Clow, 2010). The domain spans lower montane (2,700 m) to alpine (3,500 m) systems, which are home to a mosaic of vegetation types: multiple conifer species, quaking aspen Populus tremuloides stands, riparian shrubs and diverse meadows dominated by perennial grasses and forbs, interspersed with woody shrubs at lower elevations (Langenheim, 1962). This ecosystem has distinct seasonal phases that affect vegetation phenology across the mountainous landscape: snowmelt, pre-monsoon, monsoon and post-monsoon. Snowmelt and leaf out occur during April and May over most of the watershed, when air temperatures rise and a large pulse of water enters the system. Many plants in these mountainous systems grow rapidly to maximize a short growing season. Late June and early July is typically a dry period following snowmelt and pre-monsoon when most vegetated land area in the watershed is between green-up and senescence ( Figure 2). During this pre-monsoon period, there is infrequent rainfall and a period of soil drying where water limitation depends on slope, aspect and soil properties, in addition to time since snowmelt (Berdanier & Klein, 2011;Harte et al., 1995).
This particular year (2018) was characterized by a low snowpack, and consequently snowmelt (at the proximal Butte SNOTEL station) was 15 days earlier than the 1981-2010 median ( Figure 2).
Monsoon normally onsets in July, and the dry down and senescence of vegetation occurs post-monsoon in August and September. The seasonal dynamics of plant growth and impact of pre-monsoon drying can vary greatly across elevation, aspect and hillslope position (Wainwright et al., 2020). This site-specific information was essential when considering appropriate timing for imaging spectroscopy data acquisition and field survey collections, as detailed below.

| Airborne survey and field collection considerations
The NEON AOP survey, and collection of co-aligned ground-based datasets, required a combination of advanced planning and realtime decision-making. We first outline planning that was required in advance of the data collections to determine timing, sample organization procedures and identify sampling areas-the general landscapes-where field sampling could occur on a given day. We then describe methodologies implemented by the field and flight teams during the surveys, including local selection of sample sitesmeadow plots or individual plants-that were selected within a day's sampling area.

| Airborne planning and preparation
Scheduling the AOP survey during spring of 2018 required selecting a flight window using snowmelt trajectories from prior years with similar snowpack, because vegetation phenology is significantly influenced by the amount of snow and snowmelt timing (Ernakovich et al., 2014;Harte, Saleska, & Levy, 2015). Coordination with the flight team to determine flight timing, duration, and specifications began in February.
We used the Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat records of normalized difference vegetation index (NDVI) at selected sites over a range of elevations within the study extent to determine the approximate window between the onset of peak greenness and senescence during recorded years. We because generating foliar trait maps was one of the primary interests of many researchers, and peak greenness and high LAI conditions are best for these predictions (Asner & Martin, 2008). Based on estimated timing of peak NDVI, a desire to hedge on the early side of greenness to avoid the chance of early onset monsoonal rains, and aircraft availability, the airborne survey was scheduled for 13 June-3 July. Using historical weather patterns of the survey area, the NEON team built in three back-up days for every one day required to complete the full survey area to allow for down days when cloud cover would prevent data collection. Low clouds directly obscure the land surface and diffuse cloud cover and aerosols have nonlinear impacts through the VSWIR spectrum that are challenging to correct for in processing (Thompson et al., 2014). As a result, species and trait mapping is most accurate when data are collected in cloud-free and low aerosol conditions. Flight plan specifications included allowable time of flight, orientation of flight paths and height above-ground level (a.g.l.) of data collections. To achieve best illumination conditions, we implemented the AOP programme standard of restricting data collection to periods when the solar elevation angle is greater than 40° (10 a.m.-2 p.m. in this location). We considered multiple flight orientations to minimize aircraft flight time spent turning (which can significantly affect total flight time and therefore cost) and prioritize collections at lower elevations first to follow phenological change. Ultimately, we decided on a standard north-south orientation with 30% flightline overlap. Finally, our target ground resolution was 1-m resolution for both VSWIR and lidar DEM pixels once processed, for a nominal survey height of 1,000 m a.g.l. However, the complex topography of the Upper East River meant that despite efforts to adapt flight altitude to match terrain variation, the path length between ground and sensor ultimately varied by up to 1,500 m about the 1,000 m a.g.l. average in a small number of flightlines. We took steps to address these path length issues during atmospheric correction in post-flight data processing (described below).
Additional details regarding the AOP collection can be found in the NEON team's post-survey report ).

| Fieldwork planning and preparation
In early spring, we began planning our field sampling scope, methodologies and organization scheme. In addition to the research F I G U R E 2 (a) MODIS normalized difference vegetation index (NDVI) for the study region over the growing season during recent low snowpack years used for determination of peak greenness window. NDVI for 2018, along with the window of the actual Airborne Observation Platform (AOP) survey, are also shown. (b) Snow telemetry (SNOTEL) site data from the nearest location, where historical low snowpack years are coloured and the year 2018 is displayed in bold. These data are plotted relative to the day of a standardized water year starting on 1 October interests of our group, we surveyed researchers active within the study area to understand what data products or inferences from the AOP would be useful to their work. We received feedback indicating many researchers within the ecological community were interested in foliar trait and species maps, and many in the critical zone and watershed science communities were interested in extrapolating soil C, N and microbial trait characterizations.
Based on this input, we designed an integrated sampling campaign to collect vegetation, soil physical/chemical, microbial and geophysical data at co-located sampling sites to coincide with the timing of the AOP survey. We planned to collect samples within 72 hr of overflight, and without intervening rainfall, to the extent possible, given the rapid and sometimes unknown rates of change in processes of interest (Schmidt et al., 2007;Sorensen et al., 2020).
Once we identified the sample types of interest and analyses that we wanted to be able to complete with this sampling effort, we established a systematic plan for sample tracking using unique identifiers, organization, and proper storage and processing, summarized in Figure 3. This avoided confusion, mislabelling or mishandling when collecting and subsetting samples across a series of field teams. We made prelabelled field sampling packs and standardized field metadata sheets for each sampling team 2 weeks in advance of collections to avoid any duplication of sample names and ensure all supplies would be available and appropriately labelled. In addition, all supplies necessary for microbial biomass analysis and preservation for soil DNA extraction were prelabelled and shipped to RMBL.
We established appropriate sampling methods for each type of physical sample and developed sampling flowcharts to document field methods for sampling. These methodologies are described briefly in the following section and detailed protocols and flowcharts are provided in Supporting Information 1. The number of sampling sites was determined based on the amount of canopy samples desirable for sufficient ground reference data for developing statistical models to map foliar traits. Based on sample sizes in successful canopy-scale trait modelling projects (Chadwick & Asner, 2016;Singh, Serbin, McNeil, Kingdon, & Townsend, 2015), previous modelling efforts (Asner et al., 2011), and logistical realities, we aimed to collect samples from 400 sites across the 330 km 2 area: 200 meadows, 100 trees and 100 shrubs.
To achieve this, we planned for 10 days of collections during the AOP survey, with daily collections of 20 meadow plots (1 × 1 m), 10 trees and 10 shrubs (total of 40 sites) per day within a sampling area.
We identified 15 potential sampling areas in advance that could be selected for field survey during the course of AOP flights. When identifying areas for potential sampling, we selected areas spanning the elevation range of ecosystems within the watersheds, representing different geology, and were co-located with existing field data collections or monitoring when possible. An important consideration was working with RMBL to ensure we could obtain access F I G U R E 3 Flowchart detailing the handling of samples in the field, preprocessing work in the field-based laboratory, and storage and shipment. Each individual tree, shrub and meadow plot location was tagged and assigned an internal identifier during collection and an International Geo Sample Number (IGSN) identifier was assigned to each associated physical sample after processing and storage and permission for the type of sampling we planned, and therefore land ownership was also an important criterion. Finally, we wanted to ensure that field sampling occurred in locations that would span multiple days of data collection, to ensure that statistical model development would incorporate reflectance data from a variety of flight conditions and consequently lead to more robust and general models.
We organized a field team of 32 individuals (15-20 per day) to participate in the field collections and process samples as needed for preservation in the laboratory at RMBL. To standardize the sampling procedures for each area and account for turnover in field team members, we conducted training in advance, distributed detailed protocols and flowcharts to each sampling team, and scheduled people to maintain continuity of at least one team member ( Figure 3; Supporting Information 1). We planned for four teams to conduct physical sample collections each day, two for meadow plots and one each for tree and shrub individuals. In addition to physical sample collection, we planned for one team to rotate across the sampling area, visiting all sampling sites to measure soil moisture and depth, and to collect high-accuracy GPS data using a Real-Time Kinematic (RTK) enabled GPS receiver. In addition, across a day's sampling area, this rotating team collected geophysical parameters to gain a more spatially resolved dataset on subsurface soil physical characteristics including soil texture and moisture. Finally, we planned for post-airborne survey follow-up visits during the subsequent 3 months to all sampling sites for collection of additional subsurface soil samples. These visits allowed us to capture greater depth intervals than were possible during the initial field campaign for analysis of temporally stable soil properties.

| Field and airborne survey implementation
Daily decisions on whether flights would take place were coordinated between the flight team and the field team leads. Prior to field sampling, the flight team assessed weather conditions and forecasts for the day from their base of operations in Boulder, CO. They then contacted the field team lead by phone to discuss the conditions at the site and make a final determination on go/no-go for that day. Our selection of daily sampling areas was determined dynamically from our 15 pre-identified sites across the study region based on the location of the previous day's AOP survey.
To ensure that sampled vegetation was not shaded by topography, clouds or adjacent vegetation during the time of overflights, the AOP team provided preliminarily processed data within 24-36 hr of the day's survey. These data included a coarse orthorectification and processing of radiance data to four-band, red, green, blue and near-IR images. Our teams were equipped with iPads with these data loaded to enable notation of approximate locations of sampling and to photograph the site, to augment the centimetre-scale accuracy GPS location collected independently (process documented in Supporting Information 1).
When selecting individual sampling sites at the day's sampling area, each team chose plots or individuals with a high cover of photosynthetically active vegetation and homogeneous characteristics relative to the immediately surrounding area (within ~1.5 m) to decrease the confounding impact of neighbouring pixel reflectance (see Supporting Information 1 for additional details).
In addition, within a sampling area we selected a set of sites that represented a range of species and species assemblages, which was relatively easy in this system with few dominant tree species and patchy meadow environments. Leaf samples for determining foliar traits and species identities were collected for statistical model development. The tree and shrub sampling methods were analogous to those developed and utilized by the Spectranomics group (Asner & Martin, 2009. For meadow sites, we decided on foliar sampling procedures similar to clip strips, 0.1 × 2 m foliar sampling areas, implemented at core NEON sites (Weintraub & Hinckley, 2018), though we chose to sample the top 10 cm of green vegetation across the full 1 m 2 plot rather than the strip orientation.
In addition to field data collected to build statistical relationships with airborne acquisitions, additional samples were collected at each location to develop a sample archive which will expand the inferences that we may be able to draw from a single AOP flight.
These samples included, litter, roots, bulk density sample from 2 to 7 cm depth, loose soil from 0 to 10 cm and associated microbial samples. The details of the collection methodologies for each of these sample types are documented in Supporting Information 1. To our knowledge, this level of co-aligned sampling has never been undertaken in tandem with AOP surveys and subsequent studies utilizing this archive will shed considerable light on the breadth of applications for which imaging spectroscopy data can be leveraged.

| AIRBORNE AND FIELD SURVE Y COLLEC TION RE SULTS
The integrated aerial and field surveys resulted in the development of coordinated and complementary, high-resolution datasets that span the full extent of the study region. In this section, we document the data and sample archive generated directly from the paired surveys.

| Remote sensing datasets
The AOP survey was completed over 7 flight days between June 12 and 26, with a total of 102 individual flight lines collected. This collection resulted in the 72 lines required to span the survey area being surveyed under <10% cloud cover conditions, which allowed us to make a wall-to-wall mosaic of VSWIR data for the study area.
NEON provided VSWIR data that had been processed through their standard processing chain . In the following section, we discuss additional measures we took, customizing the processing of these data to improve them based on our site conditions.

| Ground reference data
We exceeded our field survey goal, collecting samples from 437 sample sites within 72 hr of AOP survey across 12 of the 15 possible study areas within the sampling domain. Sampling areas are shown in Figure 1a and an example of the distributions of sites within a sampling area are shown in Figure 4 insets. These collections took place over the course of 3 weeks in June of 2018. In addition, we collected follow-up samples from 40 conifer and riparian willows sites across the domain within the two subsequent weeks, because these vegetation types were less susceptible to the rapid phenological change experienced by meadow systems. These additional samples increased the total sampling areas to 13 in total and, in all, resulted in 487 foliar samples and sites with documented species cover . These site data, combined with data from foliar chemical analysis  and aggregated leaf area and mass data , were used as ground reference data for the development of the models for trait mapping discussed in the following section. Details of foliar laboratory analysis and LMA calculations are found in Supporting Information 1. The substantial upfront effort to develop a coordinated sample identification system, combined with procedural flowcharts (Figures S3-S5), allowed us to collect all specified samples from 437 locations with minimal sample loss or mislabelling.

| Sample archive development
Careful sample tracking and management resulted in a geolocated sample archive of plant foliar tissues, litter, roots and soil samples.
The surface soil samples from the primary 437 sites were split for immediate analysis of microbial biomass, air drying for further chemical analysis, and subsampled and archived at −80℃ for rhizosphere and bulk soil microbial DNA analyses. A bulk density sample was collected from the majority of sites, except where not possible due to rocky or other conditions, for both density calculation and soil texture analysis. The follow-up subsurface soil sampling, which took place through October, resulted in a total of 1,660 soil samples, which were air-dried and are currently being analysed. In addition, we collected data on soil moisture and depth at the time of sampling for inclusion as metadata.

| Geophysical sampling
Local scale (~100 m extent) geophysical datasets including soil apparent electrical conductivity and soil moisture were also successfully collected at the 12 sampling areas (Supporting Information 1). The soil apparent electrical conductivity within multiple depth intervals between 0 and 1.8 m was collected using a multi-coil electromagnetic induction (EMI)-based system. The soil electrical conductivity is primarily influenced by soil water content, fluid EC and grain surface conductivity. These data can be combined with soil samples to provide spatially continuous estimates of soil texture and moisture at multiple depths. These datasets encompass the sampling sites and immediately surrounding areas, and will provide a valuable

| Custom imaging spectroscopy and lidar data processing
To correct and mosaic the imaging spectroscopy data with a transparent and reproducible process, we summarize our procedure here with additional information and specifications on all files referenced here provided in Supporting Information 2. We used raw space (non-orthorectified) radiance data as our starting dataset . The prior processing steps are extensively documented by NEON (Gallery & Leisso, 2014). In brief, to generate these data, NEON processed the VSWIR data from digital number readings from the imaging spectrometer into radiance values using rigorously determined calibration parameters. Using a best-estimate of the three-space position and orientation of each sensor on the VSWIR focal plane along with a surface model from the lidar, they produced input geometry (IGM) and geometric lookup table (GLT) files for translation from raw to orthorectified space. In addition, observation conditions (OBS) files were provided that documented path length (sensor to ground in m), ground-to-sensor azimuth angle, ground-to-sensor zenith angle, ground-to-sun azimuth angle, ground-to-sun zenith angle, phase, surface slope, surface aspect, cosine i and gps time for every VSWIR pixel (Kampe, Gallery, Goulden, Leisso, & Krause, 2016).
Radiance data provide a full-spectrum reading of light received at the sensor during flight. We atmospherically corrected these radiance data to generate estimates of the apparent reflectance spectra that would have been observed at the land surface without atmospheric effects. To do this, we used the Atmospheric CORrection Now (ACORN-6LX, Imspec LLC, Glendale) software package, which uses MODTRAN lookup tables and estimates of atmospheric and surface water derived from wavelengths in the 950 and 1,150 nm regions to estimate surface reflectance. The simultaneous estimation of water vapour and surface liquid water improves the water vapour accuracy, and therefore surface reflectance estimation ( Figure S7). ACORN does not utilize input data for individual pixels on path length, elevation, timing and view-angle geometry, nor does it automatically estimate the aerosol density within the atmosphere.
Because these values can change significantly across a flightline we utilized a custom method that calculated reflectance in 200 m kernels with all local topography and view-angle geometry information, effectively accounting for variable flight conditions within each line that would highly influence the estimates of surface reflectance.
This process is an improvement over standard implementation of atmospheric algorithms and is detailed in Supporting Information 2. Implementation code repositories are referenced in Supporting Information 4, though ACORN, as with most atmospheric correction software, is proprietary.
We note that this procedure worked well at this site due to relatively mild atmospheric conditions (low water vapour and aerosols), and because of the persistent presence of vegetation which we could leverage to estimate aerosol depth. In other regions with greater aerosol content or high levels of humidity, more sophisticated atmospheric correction approaches may be needed (Thompson et al., 2018). Additional steps, including corrections for bidirectional reflectance distribution function (BRDF) effects or spectral smoothing may be applied (Thompson et al., 2019) as desired, but were not in this case. The outputs of this process are estimates of apparent surface reflectance, water vapour (wtrv), estimated atmospheric visibility in km and liquid water (wtrl, also called canopy water content or CWC over vegetation). The GLT files are subsequently used to convert these images, as well as images containing the observation parameters (OBS), from raw space to map space.
Once converted to map space, we selected flightlines to prioritize for inclusion in the complete mosaic. Some areas were flown multiple times if visibility conditions were suboptimal during an initial collection, or as part of flight operations. We chose the flightlines that aligned as closely in time to the field collections as possible and such that no rainfall events occurred between the time of sampling and the airborne surveys. In addition, when possible, we chose lines collected during the highest solar elevation angles to decrease the amount of topographic or vegetation shading within the line. In cases where clouds or cloud shade was visible within the line, we manually delineated that area and masked it out, using data from other lines to backfill. The flight plans were designed with approximately 30% overlap from one line to the next, which allowed us to mosaic the lines together using rules that specify the selection of the minimum angle between the sensor and the sun (minimum phase angle) at the time of data collection.
Finally, we generated shade masks using a custom ray tracing algorithm that estimates what pixels would have been shaded at the time of AOP data collection at the site (Figure 4). This is generated using the solar and sensor view angle geometry data (along with time of year and specific geolocation information) in the OBS data, IGM data and the digital surface elevation model . This concept was first implemented by the GAO programme (Asner et al., 2007) and calculates a path traced between the sun, each point on the ground, and the sensor, and any point that intersects the surface is considered shaded. Additional details are provided in Supporting Information 3 and links to code repository in Supporting Information 4. We found that this model estimates high levels of shading in forested regions, possibly due to the beam divergence of the lidar used (0.8 mRad), and therefore we consider it to be a conservative mask. Given that this mask potentially removes pixels that are partially sunlit, it may not be appropriate for all applications, but it does help ensure well-illuminated pixels can be identified for model development.
To make these data available to the community in a way that allows for open access and viewing without specialized software and without needing to download these large data volumes, we have made these datasets available in their processed and mosaiced form on Google Earth Engine (GEE; . This is the first Assignable Asset project to utilize this platform for dissemination, but the capacity exists for others to do so as well, and this platform could be more heavily leveraged by the community as a mechanism of free data dissemination. This platform allows for visualization, selection, extraction and download of portions of this dataset without access to high-performance computing, though the mosaic is also available through the DOE's ESS-DIVE repository (Varadharajan et al., 2019).
Accessing these data through GEE also facilitates rapid aggregation of additional remotely sensed datasets that can be used to extrapolate predictions made using NEON AOP data in either time or space.

| Trait model development
We present a case study for the types of maps that can be generated utilizing the NEON AOP platform in this type of biome, with an overview of the methods we implemented here and detailed steps and code references in Supporting Informations 3 and 4.
With the combined use of high-accuracy GPS data and GPS points recorded in the iPad software based on the preliminary imagery that we received from the AOP team, we generated polygons that represent the spatial extent of each sampling site . For meadow plots, we collected high-accuracy GPS points from each corner, and we defined a polygon for the plot based on the pixels that had the most spatial coverage within the corner points. For shrubs and trees, we outlined the extent of the crown for each individual that was sampled after field work was complete and the spectral data were processed through full orthorectification ( Figure 4). This method allows for the selection of all pixels that are associated with the individual that was sampled, rather than selecting an arbitrary distance from the GPS point that was collected.
Utilizing an estimated crown diameter to automatically define pixels is problematic because of the non-uniform nature of tree canopies.
This procedure allowed us to circumvent the challenging problem of absolute geospatial accuracy of both field and aircraft data by utilizing field and expert judgement to determine relative alignment between datasets and identify pixels for extraction. We extracted spectral information from each crown with a centroid-based extraction method, code reference in Supporting Information 4.
Once the spectral data were extracted from the sampling site polygons, we separated vegetated pixels into needle and non-needle classes to generate a classification map based on the spectral differences between these leaf types ( Figure S10). We trained a deep learning model with custom architecture, detailed in Supporting Information 3, using the site level VSWIR data, combined with additional hand delineated areas containing non-conifer land cover (Chollet et al., 2015;Abadi et al., 2015). The model performed with 0.998 true positive rate and 0.982 true negative rate, with 'positives' being non-needle identification. We assessed false positives relative to the number of conifers (e.g. number of falsely identified conifers divided by the number of conifer training pixels), due to the class imbalance, and accepted a relatively high false-positive rate of 0.2, though we note that some of these false-positive points were shaded pixels, inflating this estimate. This model was applied to all flightlines and uploaded to Google Earth Engine (GEE) for mosaicing.
We then utilize a partial least squares regression (PLSR; Feilhauer, Asner, Martin, & Schmidtlein, 2010;Wold, Sjöström, & Eriksson, 2001) approach to map foliar traits across the study area by leaf type (Figure 4b). PLSR models have been used to map canopy foliar traits in a variety of systems (Asner, Martin, Anderson, & Knapp, 2015;Dahlin, Asner, & Field, 2013;Singh et al., 2015); however, additional methods have been employed for this type of modelling, including radiative transfer to estimate a limited number of canopy level traits from leaf-level spectra, ensembles of statistical methods, and Bayesian approaches (Feilhauer, Asner, & Martin, 2015;Ferreira et al., 2018;Wang et al., 2019). Detailed information on the PLSR model metaparameters used to develop the maps presented here are provided in Supporting Information 3.
To map uncertainty in the model predictions, we generated 10 different models for needle and non-needle leaf species using different testing holdout sets of discrete sites (Singh et al., 2015;Wang et al., 2019). Each of these models was developed with a 100-fold cross validation procedure that utilized a 70% training set and 30% validation set with each fold, and then assessed based on the 10% of testing sites that were not included in that model's development (more detail in Supporting Information 3 and code reference in Supporting Information 4).
Fits of the PLSR testing set predictions from the iterative model development are shown in Figure 5. These fits are determined with the final combination of the needle and non-needle leaf models. The foliar N, LMA and leaf water content (LWC) models perform well, which is consistent with models from many other studies (Serbin, Singh, McNeil, Kingdon, & Townsend, 2016;Singh et al., 2015). Foliar C is not well predicted across this study region, which is likely due to the variety of structural and non-structural forms that C takes in leaves across the wide range of vegetation types considered here, which gives the potential for widely varying effects on the reflectance spectrum.
Models of δ 13 C also do not perform well, likely for similar reasons. It is probable that adaptations that alter δ 13 C, often associated with water use efficiency (Farquhar, Hubick, Condon, & Richards, 1989), would have inconsistent impacts on the reflectance spectrum.

| Trait mapping and model uncertainty
We applied each of the 10 models across the VSWIR mosaic and used the mean value as the ensemble predicted value for each pixel.
Standard deviation was also mapped to identify areas that have unstable predictions across models. We did not shade mask these model applications because our developed shade mask was conservative and might remove otherwise suitable pixels. By generating maps of uncertainty between models, pixels that have unstable predictions can be removed using those criteria or the shade mask at the user's discretion. We uploaded these maps to GEE to integrate the needle and non-needle leaf models, mosaicing, and masking of non-vegetated pixels (Supporting Information 3).
Maps of foliar traits and the standard deviation of their predictions utilizing these methods are available publicly at 1-m resolution on GEE (Supporting Information 4; . Modelling methods and structures utilized for generating these models and their mapping products are constantly being upgraded and improved. We did not completely exhaust efforts to optimize these trait models. These maps, with a clear understanding of the methodologies and associated error estimates, can be used for assessments across the study area both across and within vegetation types ( Figure 6).
Examples of situations that can result in high variability in foliar trait estimates between models are shown in Figure 7. An area that was partially in cloud shade at the time of overflight results in an area having generally high variation in N estimates, and can be used to identify an area that is unsuitable for use in further analysis and interpretation While calibrating spectrometer acquisitions from raw data to radiance values, the AOP team develops a stochastic model to approximate radiance at a small number of bad detector elements. These error maps provide a method for identifying and removing pixels that have unstable trait predictions and are likely unsuitable for trait modelling and interpretation. F I G U R E 5 Comparison scatter plots of field-measured foliar trait values and PLSR-estimated values at the (averaged) top of canopy level for each of the 10 PLSR models developed with linear fits (solid black for study region fits across leaf types) and the 1:1 (dashed) lines displayed. Coloured points and fits indicate needle (blue) and non-needle (orange) leaf types. Foliar N and C are in weight % and C:N is a molar ratio

F I G U R E 6
Predicted trait values and model uncertainties. (a-c) Foliar trait maps of N (weight %), leaf water content (%) and leaf mass per area (LMA, g/m 2 ). Mapped quantities are the mean values from the 10 PLSR models (one per crownbased holdout set) predicted across vegetated pixels within the study area. The site shade mask is not applied. (d-f) The corresponding standard deviation in the predicted values across all 10 PLSR models for the study area. Non-vegetated areas are shown in white

| D ISCUSS I ON
Although tremendous information can be gained from imaging spectroscopy data, they are not a panacea. The scientific utility of the information derived from these data are a function of the quality of ground reference data, care and attention to detail in data processing, and iterations that locate and resolve issues that arise from site-specific conditions. We developed a methodology for integrating airborne and ground sampling across the Upper East River catchments, consisting of vegetation with different leaf types and life strategies, to develop a remote sensing data and physical sample archive that will allow us to leverage these data and resource investments for a wide range of ecological and Earth system studies. In addition, we present a case study of developing foliar trait maps from this archive, with associated processing code, to demonstrate a mapping application. The variable topographic conditions, daily atmospheric conditions, vegetation types and goals of future studies were accommodated by advanced planning and by sampling widely across the study domain. The methods presented here provide a blueprint for other interdisciplinary and team science studies, which can be customized as required based on the questions being posed by a particular community. Some general themes arise that are universal including the need for (a) thorough advanced planning and coordination, (b) considering seasonal rates of change F I G U R E 7 Areas with high variability in predicted traits. A small area that had cloud shadow in true colour (a), mean foliar N predictions (b) and the standard deviation of those predictions (c). Flightline where errors in radiance calibrations lead to highly variable predictions (e, f) that cannot be seen in the true colour imagery (d). Areas of manicured grass (g) have high levels of variability in predictions (i), likely because no calibration samples were collected in this type of vegetation of characteristics of interest, (c) careful processing and screening of reflectance data to ensure high input data quality and (d) iteration between model development and evaluation at the level of extracted spectra, and application of models across the full study area to identify potential issues that arise in the mapping process.
Sample labelling and tracking along with a data management system are critical for integration of such diverse datasets and must be planned before data acquisition. This will allow us to continue to integrate airborne, biogeochemical and geophysical datasets as results from more time-intensive chemical analyses become available.
It was essential to conduct substantial planning and coordination for sample handling, processing and organization well in advance of the field surveys. While we acknowledge that the scope of this project is more extensive than may be possible for many groups, due to the large number of researchers working in the Crested Butte area, the collections presented here could be scaled down or altered based on the science questions being undertaken. When determining a sampling plan, we encourage researchers to consider (a) what land surface characteristics they need to quantify to extrapolate site-level insights and ensure that they will achieve a large enough sample size for their intended statistical methods and (b) sample across the range of variation that occurs in the study area and across flightlines or other external characteristics that may result in variation in spectral reflectance.
The timing of these collections was particularly important due to the rapid phenological change throughout the season, and well-calibrated models were therefore dependent on both spatial and temporal alignment of field and airborne data. In addition, datasets such as microbial composition, soil moisture and soil nutrient status change at different and sometimes unknown rates within and across seasons (Schmidt et al., 2007;Sorensen et al., 2020), so we planned to sample contemporaneously with the AOP survey to avoid temporal mismatch. This too may present a challenge to other groups depending on the characteristics of interest, but whenever possible, we recommend collecting and stabilizing co-aligned samples of characteristics that may change rapidly (i.e. microbial communities, foliar samples) even if there is not yet funding to analyse them so that they are available for future projects. On the other hand, characteristics such as soil texture and mineralogy which are stable for long periods, or species identities of long-lived individuals, could be put off until later and only collected if desired.
Processing VSWIR data also requires consideration of local conditions and there is not, as of this writing, a one-size-fits-all methodology. Here, data were collected in a high montane, topographically complex system. As a result, we did not face many challenges associated with high levels of atmospheric aerosols and water vapour, which can cause challenges in polluted or humid lowland environments. However, we did have to contend with the variable path length issues that arise from the complex terrain. Adapting mosaicking criteria also proved critical to minimize the changes in view angle between flightlines due to the highly variable topography.
These were issues that we were able to address through careful consideration of the system and its particular requirements.
Finally, a further important consideration is the potential for introducing artefacts during watershed-scale applications of PLSR or similar models across flightlines that have inherently variable conditions. Model performance statistics, such as R 2 and root mean square error from testing set ground data, are often insufficient criteria for selection. Some models can show high performance, yet also induce substantial artefacts when applied across flightlines due to variable flight conditions. As an example, despite roughly equivalent model performance statistics, we chose to utilize brightness-normalized spectra during PLSR generation. This ensured that the model fit spectral features rather than simply overall pixel brightness (which can be very characteristic in conifers, for example), and substantially reduced inter-flightline artefacts. This was something that we could only determine with certainty by applying different models across the mosaic to observe what issues arose across the VSWIR dataset. These efforts were assisted by mapping the standard deviation of the model predictions to identify areas that had poor quality spectral data.
While there are significant logistical and financial challenges to implementing this approach, the methods presented above provide a framework for doing so in a collaborative and effective manner. We also note that many ancillary benefits were realized from this intensive field survey: scientists at a range of career stages were able to spend time together working towards a common goal which led to many productive conversations and strengthened relationships; we conducted trainings that were attended by scientists and students, many of whom were only able to participate in one day of field work, but were still able to learn about our project and methodology; and we were able to include a large number of undergraduates, many of whom continue to be engaged in the project.

| CON CLUS IONS
The challenges of integrating imaging spectroscopy data with detailed field collections are worth working through for many applications. Our methodology transformed the use of airborne datasets from single investigator-based research in plant ecology to interdisciplinary research on ecosystem ecology and CZ science, supporting a number of research questions. In particular, coordinated plant, soil and geophysical datasets can open a door to understanding subsurface controls on plant species distribution and dynamics as well as plant controls on soil biogeochemistry. The scope of information gained from these studies can be enormous, enabling a quantitative understanding of landscape distributions and variability in traits at scales and resolutions not previously possible. Here, we provide a set of tools on which others can build at this site or elsewhere to gain landscape scale perspective on questions of ecology and Earth system science.

ACK N OWLED G EM ENTS
We would like to acknowledge the many people with whom we coordinated with to make this work possible. Many people who assisted for a limited time in the field, and we deeply appreciate their time and enthusiasm: Rosemary Carroll, Marco Castaneda, Christian Dewey, Jacob Gerber, Hilary Henry, Lara Kueppers, Taylor Maavara, Hailee McOmber, Alex Neild, Thomas Powell, Chris Still and Chelsea Wilmer.
We would like to thank the full NEON AOP team and Assignable Assets programme. We thank members of the Watershed Function SFA Data Management and ESS-DIVE teams for their assistance in publishing data and registering IGSN samples: Madison Burrus, Valerie Hendrix, Zarine Kakalia and Roelof Versteeg. In addition, we thank Maya Franklin for graphical assistance. Finally, we thank two anonymous reviewers and the associate editor who gave valuable and constructive feedback which improved this manuscript. K.D.C. and portions of physical

PEER R E V I E W
The peer review history for this article is available at https://publo ns. com/publo n/10.1111/2041-210X.13463.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data associated with this manuscript are published on ESS-DIVE, the Department of Energy's data portal, and have been issued searchable DOIs as cited in the manuscript. In addition to this data portal, spatial data products are available via Google Earth Engine to support processing and visualization by a wide range of users without requiring data download and local manipulation.
Furthermore, code repositories used in the development of these data are available on GitHub and versioned through the Zenodo integration. All specific links and DOIs are documented in Supporting Information file 4 and referenced in the Data Sources section.