A Multi-Sensor Unoccupied Aerial System Improves Characterization of Vegetation Composition and Canopy Properties in the Arctic Tundra

Changes in vegetation distribution, structure, and function can modify the canopy properties of terrestrial ecosystems, with potential consequences for regional and global climate feedbacks. In the Arctic, climate is warming twice as fast as compared to the global average (known as ‘Arctic amplification’), likely having stronger impacts on arctic tundra vegetation. In order to quantify these changes and assess their impacts on ecosystem structure and function, methods are needed to accurately characterize the canopy properties of tundra vegetation types. However, commonly used ground-based measurements are limited in spatial and temporal coverage, and differentiating low-lying tundra plant species is challenging with coarse-resolution satellite remote sensing. The collection and processing of multi-sensor data from unoccupied aerial systems (UASs) has the potential to fill the gap between ground-based and satellite observations. To address the critical need for such data in the Arctic, we developed a cost-effective multi-sensor UAS (the ‘Osprey’) using off-the-shelf instrumentation. The Osprey simultaneously produces high-resolution optical, thermal, and structural images, as well as collecting point-based hyperspectral measurements, over vegetation canopies. In this paper, we describe the setup and deployment of the Osprey system in the Arctic to a tundra study site located in the Seward Peninsula, Alaska. We present a case study demonstrating the processing and application of Osprey data products for characterizing the key biophysical properties of tundra vegetation canopies. In this study, plant functional types (PFTs) representative of arctic tundra ecosystems were mapped with an overall accuracy of 87.4%. The Osprey image products identified significant differences in canopy-scale greenness, canopy height, and surface temperature among PFTs, with deciduous low to tall shrubs having the lowest canopy temperatures while non-vascular lichens had the warmest. The analysis of our hyperspectral data showed that variation in the fractional cover of deciduous low to tall shrubs was effectively characterized by Osprey reflectance measurements across the range of visible to near-infrared wavelengths. Therefore, the development and deployment of the Osprey UAS, as a state-of-the-art methodology, has the potential to be widely used for characterizing tundra vegetation composition and canopy properties to improve our understanding of ecosystem dynamics in the Arctic, and to address scale issues between ground-based and airborne/satellite observations. Remote Sens. 2020, 12, 2638; doi:10.3390/rs12162638 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2638 2 of 24

. Osprey unoccupied aerial systems (UAS) airframe, sensor suite, sensor footprints, and example data. In (a) and (b), the white, yellow, and red dots, respectively, indicate the three types of sensors applied on the Osprey, i.e., Canon EOS M6 camera, ICI 9640 P-series thermal camera, and Ocean Optics FLAME spectrometer. The green dot indicates the WiFi connector, which enables real-time data logging from the Osprey. The light blue dot is the onboard global positioning system (GPS). The field of view (FOV) of the down-facing FLAME spectrometer, and the resulting image size from a Canon camera and ICI camera are shown in (c). An example of the optical RGB photo and thermal photo (overlapped with RGB) are shown in (d), and (f) is an example spectral reflectance data from the FLAME spectrometer.

Sensors
The Osprey used three types of sensors: (1) A digital single-lens reflex (DSLR) camera (Canon EOS M6), (2) a thermal infrared camera (ICI 9640 P-series), and (3) two point-sampling spectrometers (Ocean Optics FLAME) (Figure 1a,b; technical details in Table 1). The Canon camera can collect up to seven RGB photos per second with an image size of 6000 × 4000 pixels (Figure 1c,d). The ICI thermal camera collects TIR images with a size of 640 × 480 pixels (Figure 1c,e) and a thermal accuracy of +/− 1 °C. The FLAME spectrometers collect spectra in the visible to near-infrared (350~1000 nm) wavelengths with a nominal spectral resolution of 1.5 nm (Figure 1c,f). To minimize the impacts of variable illumination during each flight on our retrieval of surface spectral reflectance, we used a dual-FLAME Figure 1. Osprey unoccupied aerial systems (UAS) airframe, sensor suite, sensor footprints, and example data. In (a) and (b), the white, yellow, and red dots, respectively, indicate the three types of sensors applied on the Osprey, i.e., Canon EOS M6 camera, ICI 9640 P-series thermal camera, and Ocean Optics FLAME spectrometer. The green dot indicates the WiFi connector, which enables real-time data logging from the Osprey. The light blue dot is the onboard global positioning system (GPS). The field of view (FOV) of the down-facing FLAME spectrometer, and the resulting image size from a Canon camera and ICI camera are shown in (c). An example of the optical RGB photo and thermal photo (overlapped with RGB) are shown in (d), and (f) is an example spectral reflectance data from the FLAME spectrometer.

Sensors
The Osprey used three types of sensors: (1) A digital single-lens reflex (DSLR) camera (Canon EOS M6), (2) a thermal infrared camera (ICI 9640 P-series), and (3) two point-sampling spectrometers (Ocean Optics FLAME) (Figure 1a,b; technical details in Table 1). The Canon camera can collect up to seven RGB photos per second with an image size of 6000 × 4000 pixels (Figure 1c,d). The ICI Table 1. Technical details of the three types of sensor applied on the Osprey. The total cost of the platform was about $18,000, including the aircraft, sensors, and other accessories (estimated based on the price in 2017). We connected the sensor suite to an onboard Raspberry Pi microcomputer (Raspberry Pi Foundation, UK) to facilitate synchronized control of the multiple sensors and data logging. Instrument management during flight was handled by our open-source Python software called the Modular Data Collection System (MoDaCS, https://github.com/TESTgroup-BNL/MoDaCS). MoDaCS provided centralized triggering of instruments, data logging, record management, and error handling, as well as remote dashboard functionality for real-time data collection monitoring from a ground station. In addition to the Canon, ICI, and FLAMEs, MoDaCS was also synced with the PixHawk for recording telemetry and trigger signals as included in flight mission plans (from Mission Planner). All flight and instrument metadata were recorded with every trigger event in JSON format to link with different data outputs and the real-time location, attitude, and heading of the Osprey.

Study Site
In July 2018, we deployed the Osprey at a Kougarok Hillslope tundra field site located in western Alaska on the Seward Peninsula that was established as part of the Next-Generation Ecosystem Experiments-Arctic (NGEE-Arctic, https://ngee-arctic.ornl.gov/; [64]). The field site is located in the interior of the Seward Peninsula and is about 103 km northeast from the town of Nome (Figure 2a). The hillslope is comprised of an exposed rocky outcrop with alpine vegetation at its summit surrounded Remote Sens. 2020, 12, 2638 6 of 24 by steep well-drained slopes that transition to gently sloping tussock tundra and alder savanna, and then to lowland wet tundra. The hillslope spans a roughly 100-m change in elevation and a variety of plant communities are present across the varying topography [45]. Alder shrublands are found along the well-drained slopes below the crest of the hill and are interspersed with patches of willow-birch and dwarf shrub lichen tundra. The weather at the site is characterized by warm dry summers and long cold winters, with a short growing season (June~September). The mean annual air temperature calculated from the nearby historic ATLAS site (  (b) Species of which the canopy spectra were collected using a Spectra Vista Corporation (SVC) HR-1024i spectrometer (except for the dark fruticose lichen, which was mapped in the case study) and an example photo of ground sampling. The colored rectangles indicate the plant function types (PFTs) of the surveyed species.

Osprey Flights
We performed 12 Osprey flights (bottom panel in Figure 2a) at the Kougarok Hillslope field site in 2018, one of which (the bold yellow rectangle) was examined later in a case study to demonstrate the data products and present the performance statistics. For each flight, an automated flight mission was first designed in the Mission Planner from a ground control station (Getac B300, Getac Technology Corp., 15,495 Sand Canyon Rd., Suite 350 Irvine, CA, USA). Each flight was executed in routines of parallel rows. Data from the multi-sensor suite were simultaneously collected and stored in the field with MoDaCS. The forward-overlap and side-overlap (i.e., the overlap between photos on the same flight line and on parallel flight lines, respectively) for the optical DSLR camera were set to 85%, which resulted in an image overlap of about 75% for the ICI thermal camera. The flying height was maintained at ~40 m above the home location with a flight speed of 5 m/s determined based on the wind speed, battery capacity, and size of the flight region.
For each trigger of the optical DSLR camera, three photos with different exposures (0: auto, −1: underexposure, and +1: overexposure) were taken. We later visually selected the exposure that provided the best contrast between the vegetation canopy and the background (non-vegetation: shadow, bare soil, rock, and water) for mapping tundra vegetation. For spectral measurements, we used an FOV of 14° for the downward-facing FLAME spectrometer, which optimized the spectral footprint size (radius ~4 m) for integrating sufficient radiation and generating high-quality spectral signals based on the specified (b) Species of which the canopy spectra were collected using a Spectra Vista Corporation (SVC) HR-1024i spectrometer (except for the dark fruticose lichen, which was mapped in the case study) and an example photo of ground sampling. The colored rectangles indicate the plant function types (PFTs) of the surveyed species.

Osprey Flights
We performed 12 Osprey flights (bottom panel in Figure 2a) at the Kougarok Hillslope field site in 2018, one of which (the bold yellow rectangle) was examined later in a case study to demonstrate the data products and present the performance statistics. For each flight, an automated flight mission was first designed in the Mission Planner from a ground control station (Getac B300, Getac Technology Corp., 15,495 Sand Canyon Rd., Suite 350 Irvine, CA, USA). Each flight was executed in routines of parallel rows. Data from the multi-sensor suite were simultaneously collected and stored in the field with MoDaCS. The forward-overlap and side-overlap (i.e., the overlap between photos on the same flight line and on parallel flight lines, respectively) for the optical DSLR camera were set to 85%, which resulted in an image overlap of about 75% for the ICI thermal camera. The flying height was maintained at~40 m above the home location with a flight speed of 5 m/s determined based on the wind speed, battery capacity, and size of the flight region.
For each trigger of the optical DSLR camera, three photos with different exposures (0: auto, −1: underexposure, and +1: overexposure) were taken. We later visually selected the exposure that provided the best contrast between the vegetation canopy and the background (non-vegetation: shadow, bare soil, rock, and water) for mapping tundra vegetation. For spectral measurements, we used an FOV of 14 • for the downward-facing FLAME spectrometer, which optimized the spectral footprint size (radius~4 m) for integrating sufficient radiation and generating high-quality spectral signals based on the specified flying height and speed. A dark current measurement (i.e., the electric current from the spectrometers) was recorded before each flight with the lens caps placed on both FLAMEs. We then collected the white-reference spectra of a 10 inch × 10 inch Spectralon reflectance standard plate (LABSPHERE, INC. 231 Shaker Street, North Sutton, NH, USA) to convert the radiance of each triggered spectrum to surface reflectance.

Ground Measurements
To assist with georeferencing Osprey flights, we deployed eight ground control points (GCPs) within the survey region. We recorded the coordinates of each GCP using a hand-held differential GPS (dGPS, Trimble Geo7x, 935 Stewart Drive, Sunnyvale, CA, USA). In addition, coordinates for the centers of 25 individual tall shrubs or distinct low shrub patches were also collected. We processed the dGPS data using Pathfinder (Trimble Inc., Sunnyvale, CA, USA), and the final registration error was less than 10 cm. To validate the spectral reflectance collected by the Osprey, in situ canopy spectra of seven representative species (Figure 2b; Table 2) were collected using a Spectra Vista Corporation (SVC) HR-1024i spectrometer. These species were further classified into five PFTs and are as follows: ). Specifically, we pointed the SVC over relatively homogeneous parts of selected targets and collected multiple spectra for each species (at least three depending on the size of the target). Using the same method, we also collected the spectra of two non-vegetated surfaces (metasedimentary schist rock and water). The mean spectrum of each species was then calculated for usage in later Osprey spectral validation, described in Section 2.4.2 below.

Osprey Data Processing
Our workflow to process the multi-sensor data streams provided by the Osprey used a range of software. The workflow included two main parts ( Figure 3). Part I: The Osprey raw data collections were processed to generate primary products, including ortho-mosaicked RGB and TIR images, a digital surface model (DSM), and the canopy reflectance spectra samples. Part II: We produced a digital elevation model (DEM) based on geometric interpolation and vegetation mapping, from which we generated a canopy height model (CHM) by subtracting the DEM from the DSM.

Processing the Osprey Raw Data Collections
To generate ortho-mosaicked RGB and TIR images, we used the SfM algorithm in the Metashape software (previously called PhotoScan; Agisoft LLC, 11 Degtyarniy per., St. Petersburg, Russia, 191144). SfM generates a three-dimensional "point cloud" of surveying targets based on dense image matching and stereoscopic photogrammetry principles. It has been widely used in UAS applications for characterizing land surface structure [20,58,66]. As the same method applied for both RGB and TIR, here we demonstrate the image processing in Metashape using just the RGB as an example. First, we extracted the center location (i.e., latitude, longitude, and altitude) of each DLSR photo and the

Processing the Osprey Raw Data Collections
To generate ortho-mosaicked RGB and TIR images, we used the SfM algorithm in the Metashape software (previously called PhotoScan; Agisoft LLC, 11 Degtyarniy per., St. Petersburg, Russia, Remote Sens. 2020, 12, 2638 9 of 24 191144). SfM generates a three-dimensional "point cloud" of surveying targets based on dense image matching and stereoscopic photogrammetry principles. It has been widely used in UAS applications for characterizing land surface structure [20,58,66]. As the same method applied for both RGB and TIR, here we demonstrate the image processing in Metashape using just the RGB as an example. First, we extracted the center location (i.e., latitude, longitude, and altitude) of each DLSR photo and the corresponding aircraft attitude (i.e., pitch, roll, yaw) from the MoDaCS software (JSON files in Figure 3I). In Metashape, we then generated geo-referenced dense point clouds by linking the image center location and aircraft attitude to the DLSR photos. Next, by following a step-by-step workflow (Agisoft, 2018), we derived the ortho-mosaicked RGB image. In this process, a DSM image was also generated. After that, we further corrected the registration error by manually georeferencing the RGB and DSM using in situ GCPs. In this study, we applied 'high' for generating dense cloud points. The output spatial resolution of RGB, TIR, and DSM (from RGB processing) was~1, 5, and 2 cm, respectively, using default export settings.
The reflectance spectra samples collected by the dual-FLAME system were extracted via the MoDaCS software. Using a locally estimated scatterplot smoothing (LOESS) algorithm (Cleveland, 1979), we resampled the data to a consistent 1 nm spectral resolution. To visualize the spectral measurements, we further extracted the footprint (i.e., the corresponding content of each spectrum on the surface, Figure 1c) of each spectral measurement from the ortho-mosaicked images using the ENVI+IDL software (Harris Geospatial Solution Inc.). For this, we first located the center (i.e., latitude and longitude) of each spectral measurement on all the previously mentioned images. The radius of the ith footprint (R i ) was then calculated using Equation (1): where H i,UAS is the altitude of the Osprey and E i,ground is the altitude of the ground/canopy derived from the DSM. The FOV is 14 • in this study. Once the radius was obtained, we then clipped out the range of the footprints from all images.

Deriving the Canopy Height Model (CHM)
The Metashape DSM contains both ground and top-of-vegetation data points and requires additional processing to be used for investigating vegetation structures. Here, we present our approach to derive CHM from the DSM with three steps ( Figure 3II):

1.
Identify 'on-ground' pixels from the DSM: We first performed an object-based classification on the four-band layer-stack of RGB and TIR (details for classification can be found in mapping tundra vegetation below). The 'on-ground' pixels were then identified as non-vegetation classes, e.g., schist rock and water, or non-vascular vegetation, e.g., lichen. By doing this, we avoided including any vegetation canopies that have a discernible height into the later reconstruction of the DEM, which serves as a base layer for deriving CHM.

2.
Interpolate the DEM using the 'on-ground' pixels: To produce a spatially continuous smooth DEM (i.e., without surface structure), we used a thin plate spline (TPS) algorithm [67] in ENVI+IDL to interpolate the 'on-ground' pixels. TPS provides smooth interpolation of given scattered data in two or more dimensions by considering geometric properties [68,69] and has been widely used for creating DEMs from LiDAR cloud points [70,71].

3.
Calculate the CHM: Once the DEM was obtained, we calculated the CHM by subtracting the DEM from the DSM. The resulting CHM has the same spatial resolution as the DSM.
Using this workflow, we processed the Osprey flight data collected from the Kougarok study site. The results included four image products (i.e., ortho-mosaicked RGB and TIR, DSM, and CHM), point-sampled surface spectral reflectance, and the footprint for each spectral measurement. In this study, we validated the obtained DEM against the dGPS measurements (see ground measurements) using linear regression. R 2 and root mean square error (RMSE) were calculated to assess their agreement. However, we did not validate the CHM since previous studies have shown that SfM-generated CHM from digital photography can be as accurate as LiDAR products in the Arctic [21,49] and other ecosystems [58,72,73].

Methods for the Case Study
To illustrate the utility of our Osprey data products for characterizing vegetation canopy properties, we chose a single flight (bold yellow rectangle in Figure 2a, flight date: 25 July 2018) from our data collection at the Kougarok Hillslope and performed a case study. The flight covered a total area of 1.8 ha that includes a variety of plant species. Using the data processing methods described above, both image (RGB, TIR, DSM, and CHM) and spectral products (167 spectral measurements) were acquired and processed. We first mapped the distribution of select dominant individual species and PFTs using an object-based image classification. Based on the species and PFT maps, we then (1) evaluated the quality of Osprey spectral reflectance; (2) characterized the canopy-scale height, surface temperature (i.e., skin temperature), and greenness of different PFTs; and (3) analyzed the spectral separability among of deciduous low to tall shrub PFT fractional covers.

Mapping Tundra Vegetation
Tundra species and PFTs were mapped using the four-band layer-stack of RGB and TIR. CHM was not used because deriving CHM requires image classification as a model input. Prior to the classification, we first resized and resampled the RGB and TIR to the same spatial extent (i.e., the TIR extent as it is smaller than RGB) and resolution (1cm, i.e., the RGB resolution as it is higher than TIR). To reduce 'salt-and-pepper' noise (Boncelet, 2009) in the high-resolution Osprey images, we employed a multi-resolution segmentation algorithm (scale = 20, shape = 0.2) in eCognition 8.9 (Trimble Inc.) to segment the layer stack of RGB and TIR to generate relatively homogeneous "objects". Eleven features related to the objects were generated, including the object mean and variance in blue, green, red, and TIR, respectively; object area; object border length; and object-averaged Green Chromatic Coordinate (GCC, Equation (2)) which is an effective indicator of vegetation greenness (Leduc and Knudby 2018): where B, G, and R are the values of blue, green, and red channels of the optical DSLR camera.
To perform classification, we used a supervised maximum likelihood classifier (MLC) in ENVI+IDL for its demonstrated good performance [74,75]. The 11 object features were used as input for the MLC classifier. Classification training samples were visually selected from the layer stack of RGB and TIR with the aid of in situ survey information and our expert knowledge about the study site. We generated a sample dataset with~3600 eCongition-provided objects, including eight species classes (Siberian alder, Arctic dwarf birch, bog blueberry, black crowberry, Alaskan mountain-avens, tussock cotton-grass, reindeer lichen, and dark fruticose lichen) and three non-vegetation classes (shadow, schist rock, and water). Other plant species [76] are also present in the case flight region; however, they are either not dominant or cannot easily be distinguished from the Osprey images. For example, Bigelow's sedge (Carex bigelowii) is known to occur at the site and was lumped with the tussock cotton-grass species class for this first case study, along with other sedges, grasses, and rushes that were collectively merged into graminoid in the later PFT-level classification. The sample dataset was randomly split into two chunks: 50% for training and 50% for accuracy assessment. The training set was then used as input to the supervised MLC classifier to map the 11 classes. We then merged our seven plant species classes into PFTs (Table 2) based on the PFT descriptions of Walker et al. (2000) [77]. The confusion matrix and overall accuracy [78,79] were calculated for both the species-and PFT-level classifications, using the 50% objects.

Evaluating Osprey Spectral Reflectance
To evaluate the Osprey-collected spectral reflectance measurements, we employed a spectral simulation approach (similar to Yang et al. (2017) [80], Figure S1). The simulation assumes that the spectrum of a footprint (i.e., the extent of a spectral measurement) is a linear mixture of its containing components [48], as shown in Equation (3): where R mix is the simulated reflectance of a footprint; f i and R i are the fractional cover and reflectance of the ith component, respectively. λ is the wavelength and n is the total number of components identified in the footprint. The simulated footprint spectra were used as a 'reference' to compare with the Osprey spectral data.
In the case study Osprey flight, a total of 167 spectral measurements were acquired. We extracted the footprint of each spectral measurement from the species classification described above and calculated the fractional cover of each species class component. The input reflectance spectrum of each component for the simulation was acquired from the in situ SVC measurements. Multiple spectra for each component were collected and their mean was used as the input of R i . Not all species present in the spectral footprints were included among the in situ SVC measurements; therefore, we conducted the spectral simulation using the seven dominant species components (Siberian alder, Arctic dwarf birch, bog blueberry, black crowberry, tussock cotton-grass, Alaskan mountain-avens, and reindeer lichen) and two non-vegetation components (schist rock and water) where we had associated SVC measurements (as described in ground measurements). Spectral measurements where the above nine components occupied over 95% of the total area of the footprint were selected. This resulted in 105 spectral simulations to be used for our spectral reflectance evaluation. The spectra of these footprints were then simulated using Equation (3). A linear regression was used to compare the Osprey-collected against the simulated spectral reflectance. R 2 and RMSE were calculated to assess the agreement between the two datasets.

Characterizing Canopy Height, Temperature, and Greenness of PFTs
Changes in vegetation composition affect energy and matter exchange through the atmosphere-plant-soil continuum, mainly via the different canopy properties among vegetation types. By aligning the PFT map with the RGB, CHM, and TIR, we investigated the canopy height, temperature, and greenness of different PTFs. The PFT-level classification was used because (1) it is broadly used in Earth system models (ESMs) to represent global vegetation dynamics [77,[81][82][83], and (2) the mapping accuracy of PFTs can be much higher than individual species, thus largely reducing the analysis uncertainty caused by misclassification (e.g., a high degree of confusion between the two deciduous low shrub species). In this study, we merged the species classes into five PFTs, including deciduous low to tall shrub (DLTS), deciduous low shrub (DLS), evergreen shrub (ES), graminoid (GR), and lichen (LI). Non-vegetation classes were excluded for the analysis. The mean and standard deviation (sd) of GCC, CHM, and TIR were then calculated for each PFT. Tukey tests were conducted to examine if there is a significant difference in canopy GCC, CHM, and TIR among PFTs.

Quantifying Spectral Separability among Deciduous Low to Tall Shrub PFT Covers
Climate warming in arctic tundra regions has led to a significant increase in low and tall shrub abundance and this increase affects energy balance by causing changes in surface reflectance [84]. In this study, we examined the spectral separability among deciduous low to tall shrub PFT covers using Osprey spectral products. To this end, we first calculated the fractional cover of low to tall shrub PFT for the 167 Osprey spectral measurements by extracting their footprints from the PFT-level classification. This PFT originated from the species-level class 'Siberian alder' comprised of mainly Siberian alder and other low to tall shrub species that were less common in the area, such as diamond-leaf willow (Salix pulchra) and resin birch (Betula glandulosa), and were lumped with Siberian alder. Based on the fractional cover, we then binned the spectral measurements into seven groups (i.e., group 1: 0-10%; group 2: 10-20%; group 3: 20-30%; group 4: 30-40%; group 5: 40-50%; group 6: 50-60%; group 7: 60-68%). Using an instability index (ISI, Equation (4); [85]), we examined the spectral separability among these groups. ISI measures the ratio of intra-group spectral variability to inter-group spectral variability at each spectral band. ISI < 1.0 indicates high separability, and ISI > 1.0 indicates lower separability: In this equation, R p,λ and R q,λ are the mean reflectance of group p and q, respectively, at wavelength λ; δ p,λ and δ q,λ are the standard deviation of ground p and q, respectively; and n is the number of groups (i.e., seven).

Osprey Image Products, Classification, and Canopy CHM, TIR, and GCC of Different PFTs
The pre-planned Osprey mission for our case study flight yielded a very high resolution (~1 cm) RGB image mosaic (Figure 4a). The TIR mosaic, while having a lower resolution (~5 cm) given the sensor characteristics, still provided very high spatial details for the flight area (Figure 4b). The resulting RGB-based DSM (~2 cm) displayed a generally smooth gradient in topography from 75 to 102 m (Figure 4c). In contrast, the TIR ( Figure 4b) and CHM (Figure 4d) images showed high spatial variation over the landscape, associated with fine-scale variation in the vegetation type. Canopy temperature varied from 18 to 30°C, where cooler canopies (dark blue regions in Figure 4b) were found at middle-to-upper elevations of the case flight area, while warmer canopies (bright yellow regions in Figure 4b) were found on the top and bottom. The maximum canopy height was about 2.5 m, with taller plants (yellow to red colors in Figure 4d) in regions where cooler canopies occurred (Figure 4b). The derived DEM was validated against the field GCPs using linear regression. An R 2 of 0.997 and RMSE of 0.25 m were achieved ( Figure S2).
The dominant tundra species mapped from the combination of RGB and TIR images (Figure 5a) had an overall accuracy of 81.5% (Table 3), although higher errors were found for the Arctic dwarf birch and bog blueberry classes (birch: producer's accuracy 78.2%, user's accuracy 63.6%; blueberry: producer's accuracy 45.5%, user's accuracy 49.7%). Clusters of Siberian alder (dark green color) were found in the middle-to-upper elevations of the study area (DSM: 85~100 m), as shown in Figure 5a. Arctic dwarf birch and bog blueberry shrubs were mixed with tussock cotton-grass and were distributed at lower elevations (DSM: 75~90 m). The mapping accuracy by PFT (87.4%, Table S1) was much higher than species, as expected. Evergreen shrubs (Alaskan mountain-avens and black crowberry) were distributed mostly on the top of the hillslope above the Siberian alder shrubs (DSM: 95~102 m) (Figure 5b). We also observed a high abundance of lichens in the area; however, they showed very similar optical features with the schist rocks found on the top of the Kougarok Hillslope (Figures  4a and 5b). Within the study area, we identified four representative plant community types (PCTs), including 1) Dryas-lichen dwarf shrub tundra, 2) Alder shrubland, 3) Birch-Ericaceous-lichen shrub tundra, and 4) Alder savanna in tussock tundra ( Figure S3) [76]. These PCTs represented several of the dominant vegetation communities at the Kougarok Hillslope site.
95~102 m) (Figure 5b). We also observed a high abundance of lichens in the area; however, they showed very similar optical features with the schist rocks found on the top of the Kougarok Hillslope (Figures 4a, 5b). Within the study area, we identified four representative plant community types (PCTs), including 1) Dryas-lichen dwarf shrub tundra, 2) Alder shrubland, 3) Birch-Ericaceous-lichen shrub tundra, and 4) Alder savanna in tussock tundra ( Figure S3) [76]. These PCTs represented several of the dominant vegetation communities at the Kougarok Hillslope site.

Classified results
Siberian alder  The canopy GCC, CHM, and TIR of different PFTs are summarized as boxplots in Figure 6. As expected, deciduous low to tall shrub had a higher mean (0.98 m) and variance (sd: 0.68 m) in the CHM, and deciduous low shrub showed a lower mean height (mean: 0.31 m) and variance (sd: 0.43 m) (Figure 6a). The Tukey test suggested no significant difference in CHM among other PFTs (p-value > 0.1). However, the canopy TIR exhibited several discernible differences among PFTs (Figure 6b). On average, lichen had the highest canopy temperature (mean: 24.2 ± 1.4°C) and deciduous low to tall shrub had the lowest canopy temperature (mean: 20.1 ± 1.5°C). A general trend of decreasing TIR with CHM existed, except for graminoid (i.e., tussock cotton-grass), which slightly departed from this trend (Figure 6b). The maximum value of GCC (mean: 0.41 ± 0.03) correlated with deciduous low to tall shrub and dropped gradually to a mean of 0.34 ± 0.01 for lichen (Figure 6c). Generally, taller shrub PFTs had higher GCCs than low shrub PFTs, and shrub PFTs had higher GCCs than non-shrub PFTs (graminoid and lichen) in the case study region.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 23 6a). The Tukey test suggested no significant difference in CHM among other PFTs (p-value > 0.1). However, the canopy TIR exhibited several discernible differences among PFTs (Figure 6b). On average, lichen had the highest canopy temperature (mean: 24.2 ± 1.4 ℃) and deciduous low to tall shrub had the lowest canopy temperature (mean: 20.1 ± 1.5 ℃). A general trend of decreasing TIR with CHM existed, except for graminoid (i.e., tussock cotton-grass), which slightly departed from this trend (Figure 6b). The maximum value of GCC (mean: 0.41 ± 0.03) correlated with deciduous low to tall shrub and dropped gradually to a mean of 0.34 ± 0.01 for lichen (Figure 6c). Generally, taller shrub PFTs had higher GCCs than low shrub PFTs, and shrub PFTs had higher GCCs than non-shrub PFTs (graminoid and lichen) in the case study region.

Osprey Spectra Products, Validation, and Spectral Separability among Deciduous Low to Tall Shrubs Covers
The reflectance spectra collected from the Osprey platform are shown in Figure 7. In general, near-infrared reflectance increased while visible reflectance declined with an increase in Siberian alder shrub cover, resulting in large landscape-scale variation in reflectance across the visible to nearinfrared wavelengths. The extracted spectra footprints had well-aligned features between RGB and TIR images, despite the differences in the nominal pixel resolution, and captured fine-scale details of the vegetation composition across the case flight region (Figure 7a-c).
We validated the Osprey reflectance spectra against simulated spectra with a linear mixture of in situ spectral measurements (Figure 8a). An R 2 of 0.95 and RMSE of 0.04 were achieved when assessed using linear regression (Figure 8b), indicating that the dual-FLAME configuration on the Osprey was sufficient to overcome confounding issues, including in-flight aircraft attitudinal changes, sun-sensor geometry, instrument noises, etc., to generate high-quality spectral reflectance products. However, higher agreements between the two datasets were seen in the visible range (higher degree of clustering around 1:1 line where the reflectance <~0.4), while there was a larger degree of scattering in the near-infrared range (i.e., where the reflectance >~0.4).
We binned the reflectance spectra into seven groups using the deciduous low to tall shrub PFT fractional cover. As expected, the group-averaged reflectance decreased with deciduous low to tall shrub PFT fractional cover in the visible to red edge range (i.e., 400-700 nm, Figure 9a), with group 1 (i.e., 0~10% cover) having the highest reflectance and group 7 (i.e., 60~68% cover) having the lowest reflectance (Figure 9a,b). Instead, the opposite pattern was observed in the near-infrared range (700-900 nm). This finding suggests the fractional cover of deciduous low and tall shrubs, which is expected to increase in response to climate change, has a significant impact on canopy spectral reflectance and can impact some wavelengths more than others. The ISI value was generally below 1.0 in the visible range, indicating a high spectral separability among these groups (Figure 9c). However, we did not find the same separability in the near-infrared range where the ISI was above 1.0 (Figure 9c). To explore this, we further examined the δ (i.e., intra-group standard deviation) and |Rp − Rq| (i.e., inter-group difference) terms in Equation (4). The result showed that, on average, the values of δ are almost twice

Osprey Spectra Products, Validation, and Spectral Separability among Deciduous Low to Tall Shrubs Covers
The reflectance spectra collected from the Osprey platform are shown in Figure 7. In general, near-infrared reflectance increased while visible reflectance declined with an increase in Siberian alder shrub cover, resulting in large landscape-scale variation in reflectance across the visible to near-infrared wavelengths. The extracted spectra footprints had well-aligned features between RGB and TIR images, despite the differences in the nominal pixel resolution, and captured fine-scale details of the vegetation composition across the case flight region (Figure 7a-c).
We validated the Osprey reflectance spectra against simulated spectra with a linear mixture of in situ spectral measurements (Figure 8a). An R 2 of 0.95 and RMSE of 0.04 were achieved when assessed using linear regression (Figure 8b), indicating that the dual-FLAME configuration on the Osprey was sufficient to overcome confounding issues, including in-flight aircraft attitudinal changes, sun-sensor geometry, instrument noises, etc., to generate high-quality spectral reflectance products. However, higher agreements between the two datasets were seen in the visible range (higher degree of clustering around 1:1 line where the reflectance <~0.4), while there was a larger degree of scattering in the near-infrared range (i.e., where the reflectance >~0.4).
We binned the reflectance spectra into seven groups using the deciduous low to tall shrub PFT fractional cover. As expected, the group-averaged reflectance decreased with deciduous low to tall shrub PFT fractional cover in the visible to red edge range (i.e., 400-700 nm, Figure 9a), with group 1 (i.e., 0~10% cover) having the highest reflectance and group 7 (i.e., 60~68% cover) having the lowest reflectance (Figure 9a,b). Instead, the opposite pattern was observed in the near-infrared range (700-900 nm). This finding suggests the fractional cover of deciduous low and tall shrubs, which is expected to increase in response to climate change, has a significant impact on canopy spectral reflectance and can impact some wavelengths more than others. The ISI value was generally below 1.0 in the visible range, indicating a high spectral separability among these groups (Figure 9c). However, we did not find the same separability in the near-infrared range where the ISI was above 1.0 (Figure 9c). To explore this, we further examined the δ (i.e., intra-group standard deviation) and |R p − R q | (i.e., inter-group difference) terms in Equation (4). The result showed that, on average, the values of δ are almost twice as large as the |R p − R q | in the near-infrared range, while being almost the same as |R p − R q | in the visible range ( Figure S4; the displayed δ and |R p − R q | are the averaged intra-group standard deviation of the seven groups and the averaged inter-group difference among the seven groups). This disproportionately large intra-group standard deviation (δ) thus leads to high ISI values in the near-infrared range, which corresponds to low spectral separability. In addition, a spike in ISI was found in the red edge range, as a result of the significantly low |R p − R q | (the ribbon highlighted region in Figure S4).
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 as large as the |Rp − Rq| in the near-infrared range, while being almost the same as |Rp − Rq| in the visible range ( Figure S4; the displayed δ and |Rp − Rq| are the averaged intra-group standard deviation of the seven groups and the averaged inter-group difference among the seven groups). This disproportionately large intra-group standard deviation (δ) thus leads to high ISI values in the nearinfrared range, which corresponds to low spectral separability. In addition, a spike in ISI was found in the red edge range, as a result of the significantly low |Rp − Rq| (the ribbon highlighted region in Figure S4).   as large as the |Rp − Rq| in the near-infrared range, while being almost the same as |Rp − Rq| in the visible range ( Figure S4; the displayed δ and |Rp − Rq| are the averaged intra-group standard deviation of the seven groups and the averaged inter-group difference among the seven groups). This disproportionately large intra-group standard deviation (δ) thus leads to high ISI values in the nearinfrared range, which corresponds to low spectral separability. In addition, a spike in ISI was found in the red edge range, as a result of the significantly low |Rp − Rq| (the ribbon highlighted region in Figure S4).

Discussion
Accurately characterizing vegetation composition, structure, and function in sufficient detail has been a key challenge in Arctic tundra regions where these conditions are often highly variable [12,35,64]. To address this challenge, we developed a cost-effective multi-sensor UAS (the 'Osprey') that simultaneously collects hyperspectral, structural, and thermal data of vegetation canopies. We provided a workflow to process UAS collections to generate ortho-mosaicked optical and thermal images, surface elevation, canopy height, and spectral reflectance data. Using a case study from a low Arctic tundra site in Western Alaska, we showed that our platform provides an efficient means for accurate mapping of tundra vegetation as well as for characterizing canopy-scale greenness, vegetation height, and 'skin' temperature. In addition, we demonstrated the potential of the Osprey system for investigating and interpreting important ecological patterns by quantifying the impacts of PFT cover, namely the abundance of deciduous low to tall shrub PFT, on canopy biophysical properties (reflectance and surface temperature).
By using an object-based classification, we illustrated that the Osprey RGB and TIR images are effective for mapping representative tundra species and PFTs. The accuracy of mapped species in this study (81.5%) was similar to that found by Fraser et al. (2016) [49], who mapped low Arctic plant species with an overall accuracy of 82%. We observed that Arctic dwarf birch and bog blueberry shrubs have similar canopy appearance because they were mixed within local communities, leading to higher mapping uncertainties between them. Since the mapped PFTs (overall accuracy: 87.4%) were derived from the individual species maps, they were more accurate than the individual species classification. Moreover, these maps, when combined with the other datasets provided by the Osprey, were effective at capturing landscape-scale variation in canopy height, temperature, and greenness ( Figure 6).
The canopy of the deciduous low to tall shrubs was found to be cooler than other PFTs, likely due to their higher transpiration rate and strong 'shading' effect [5]. On the other hand, the 'drier' non-vascular lichen had higher canopy temperatures on average than other PFTs, because they tend to have a lower water content, regardless of their high albedo [86]. An analysis of the TIR of schist rock showed that, on average, the canopy temperature of lichen (24.2 °C) was very similar to schist rock (24.3 °C), which was one of the hottest surface components in our study area. These hot, dry lichens, together with increased summer temperatures, might play a role in the increased wildfires in the Arctic [87]. We also observed a negative relationship between PFT canopy height and canopy temperature where taller canopies had cooler surface temperatures ( Figure 6), but this pattern did show some departures. For example, taller graminoids commonly mixed with deciduous low shrubs at lower elevations were found to be cooler than low-growing evergreen shrubs encompassed by schist rocks and lichens on the hillslope summit (Figures 5b and 6b). Typically, these departures were driven by strong differences in background temperature from non-vegetated or non-vascular plants. In addition, taller shrub PFTs

Discussion
Accurately characterizing vegetation composition, structure, and function in sufficient detail has been a key challenge in Arctic tundra regions where these conditions are often highly variable [12,35,64]. To address this challenge, we developed a cost-effective multi-sensor UAS (the 'Osprey') that simultaneously collects hyperspectral, structural, and thermal data of vegetation canopies. We provided a workflow to process UAS collections to generate ortho-mosaicked optical and thermal images, surface elevation, canopy height, and spectral reflectance data. Using a case study from a low Arctic tundra site in Western Alaska, we showed that our platform provides an efficient means for accurate mapping of tundra vegetation as well as for characterizing canopy-scale greenness, vegetation height, and 'skin' temperature. In addition, we demonstrated the potential of the Osprey system for investigating and interpreting important ecological patterns by quantifying the impacts of PFT cover, namely the abundance of deciduous low to tall shrub PFT, on canopy biophysical properties (reflectance and surface temperature).
By using an object-based classification, we illustrated that the Osprey RGB and TIR images are effective for mapping representative tundra species and PFTs. The accuracy of mapped species in this study (81.5%) was similar to that found by Fraser et al. (2016) [49], who mapped low Arctic plant species with an overall accuracy of 82%. We observed that Arctic dwarf birch and bog blueberry shrubs have similar canopy appearance because they were mixed within local communities, leading to higher mapping uncertainties between them. Since the mapped PFTs (overall accuracy: 87.4%) were derived from the individual species maps, they were more accurate than the individual species classification. Moreover, these maps, when combined with the other datasets provided by the Osprey, were effective at capturing landscape-scale variation in canopy height, temperature, and greenness ( Figure 6).
The canopy of the deciduous low to tall shrubs was found to be cooler than other PFTs, likely due to their higher transpiration rate and strong 'shading' effect [5]. On the other hand, the 'drier' non-vascular lichen had higher canopy temperatures on average than other PFTs, because they tend to have a lower water content, regardless of their high albedo [86]. An analysis of the TIR of schist rock showed that, on average, the canopy temperature of lichen (24.2 • C) was very similar to schist rock (24.3 • C), which was one of the hottest surface components in our study area. These hot, dry lichens, together with increased summer temperatures, might play a role in the increased wildfires in the Arctic [87]. We also observed a negative relationship between PFT canopy height and canopy temperature where taller canopies had cooler surface temperatures ( Figure 6), but this pattern did show some departures. For example, taller graminoids commonly mixed with deciduous low shrubs at lower elevations were found to be cooler than low-growing evergreen shrubs encompassed by schist rocks and lichens on the hillslope summit (Figures 5b and 6b). Typically, these departures were driven by strong differences in background temperature from non-vegetated or non-vascular plants. In addition, taller shrub PFTs were greener (i.e., higher GCCs), which was expected given their higher leaf area index (LAI) and potential higher chlorophyll content (Mentis and Barbour, 1990). Reindeer lichen (white) and dark fruticose lichen (black and dark brown) both had low GCCs, given their low chlorophyll content. Evergreen shrub species had more variable canopy color (i.e., dark green mature leaves, young red-green leaves, light yellow-green senescing leaves, and black fruits for black crowberry), resulting in lower GCC than deciduous shrubs as well. These findings confirm that the high-resolution Osprey image products are effective for characterizing tundra vegetation canopy properties, improving our understanding of tundra vegetation across the landscape.
The spectral assessment showed that our dual-spectrometer (i.e., Ocean Optics FLAME) setup yielded high-quality spectral reflectance retrievals with good consistency and stability across the flight, and was effective at accounting for the variation in illumination over the flight (Figure 8b). However, we found higher uncertainty in the near-infrared wavelengths. This may have occurred because the linear mixture simulation and ground spectral measurements do not account for the effects of a varying three-dimensional canopy structure within and across PFTs, which is higher in the near-infrared wavelengths [88], leading to a less accurate representation of simulated spectra in the near-infrared [89]. In addition, the FLAME spectrometers have a decreasing signal-to-noise ratio with an increasing wavelength, with significant noise found in wavelengths beyond~850 nm (SNR: 400/1), thus impacting the near-infrared more than the visible wavelengths. Despite these challenges, our spectral analysis showed that, as the fractional cover of deciduous low to tall shrub PFT increases, the visible to red-edge reflectance decreases, while the near-infrared reflectance increases (Figures 7 and 9). This finding aligns with previous studies [32] and suggests that an increase in deciduous low to tall shrub abundance, as predicted with climate warming, will make the land surface darker in the visible to red-edge wavelengths but brighter in near-infrared wavelengths. This does not contradict the previous finding that 'increased shrub abundance reduces the full-range shortwave surface albedo in the Arctic', as suggested by ground-based measurements [90] and radiative transfer model simulation [11,84]. On the other hand, it reveals a likely divergent response of surface albedo at different wavelength ranges to the increased shrub abundance, which needs to be considered when scaling up to the full-range surface albedo. Additionally, our analysis quantified the spectral differences among low to tall shrub PFT fractional cover groups (Figure 9a,b) and found higher separability in the visible to red-edge range and lower separability in the near-infrared range (Figure 9c). This suggests that spectral reflectance at visible to red-edge wavelengths may be valuable to determine the fractional cover of low to tall shrubs, which have dominant impacts on energy balance, resource allocation, and primary production in the arctic tundra biome. However, it is worth pointing out that this pattern may not be the case for other PFTs, such as lichen, graminoid, and deciduous low shrub, or for differentiating tundra plant species or communities where near-infrared wavelengths have shown a high spectral separability [30,91].
In recent studies, the use of canopy spectroscopy, such as hyperspectral images from AVIRIS and the Hyperspectral Infrared Imager (HyspIRI) and Surface Biology and Geology (SBG) missions, to predict vegetation composition, structure, and function has been shown to be effective in a variety of biomes [13][14][15][16][17][18]. However, similar work remains technically challenging in the arctic tundra as a result of low data availability and high spatial heterogeneity. In this paper, we demonstrated that the Osprey system is effective for detailed mapping of tundra vegetation, while providing high-quality canopy spectral measurements. This provides us with a method to link the landscape patterns in canopy spectroscopy with the spatial variation in vegetation composition, as partly suggested by the case study spectral analysis on the low to tall shrub PFT (Figure 9). This approach can help determine the drivers of variation in canopy spectroscopy, as well as addresses the scaling gaps between the ground-based spectroscopy and airborne/satellite observations [26,92]. In addition, the Osprey spectroscopy per se can be a valuable resource to validate spectral measurements provided by airborne hyperspectral imager, which is yet to be broadly validated in the Arctic.
While we demonstrated the effectiveness of the Osprey system for improving the characterization and mapping of vegetation composition, structure, and function in the tundra ecosystems, we identified some limitations. First, surface topography and aircraft attitude (pitch, yaw, and roll) were not considered for deriving the footprints of spectral measurements. This can cause a certain degree of error in the calculation of corresponding fractional cover of different species or PFTs. However, sensors on the Osprey are mounted using a powered three-dimensional gimbal (Gremsy H3) to maintain a nadir-viewing orientation and our flight records show that the Osprey is relatively stable for most of the flight missions, with the pitch, roll, and yaw all less than 1.5 • , which is much lower than other platforms, such as DJI Matrice 600 Pro (5-10 • ; [93]). Second, though we identify significant canopy temperature differences among PFTs, the absolute temperature of PFT canopies can change with time and weather conditions (e.g., illumination, wind speed, precipitation). This might affect the temperature patterns observed in this study. In the future, we will utilize ancillary measurements of incident radiation and atmospheric temperature to account for variation in temperatures during flights or to standardize by calculating the canopy minus air temperature difference as a way to account for changing weather conditions. Third, we derived the CHM by interpolating 'on-ground' pixels to generate a DEM. However, these 'on-ground' pixels may not represent the true 'bare ground' in areas with microtopography, such as tussock tundra, where lichen can occur both on top of a tussock and on the ground between tussocks. This might cause underestimation for the CHM of low-growing PFTs, though this method has been proven effective for characterizing much taller shrub PFTs [21,49]. However, we suggest not using this method, when survey point clouds are sparse or study region is dominated by highly rugged surfaces. Lastly, while the CHM provides additional information to improve the classification, we did not utilize this in the current study. Future work with the Osprey platform will validate the CHM and explore incorporating CHM into the classification more directly over a larger number of flights.

Conclusions
In this paper, we described the development of a cost-effective multi-sensor UAS for acquiring data on the spectral, structural, and thermal properties of arctic tundra vegetation. We deployed the UAS at a field site located on the Seward Peninsula in western Alaska, a location at the forefront of climate change. Based on data from our flights, we then developed a workflow to convert the UAS collections to data useful for characterizing vegetation biophysical properties. We presented our results using a case study to highlight the products we can generate to characterize vegetation composition, structure, and function. Our results showed that representative arctic tundra vegetation species could be mapped with an overall accuracy of 81.5%. In addition, significant differences in canopy-scale greenness, canopy height, and 'skin' temperature were observed among PFTs using the UAS image products. In a spectral analysis, we also showed that the spectral data has the potential to discern different fractional covers of deciduous low to tall shrubs. Therefore, the developed UAS, as a state-of-the-art methodology, has the potential to be widely used to improve our understanding of ecosystem dynamics in the Arctic.