Multiple UAV Flights across the Growing Season Can Characterize Fine Scale Phenological Heterogeneity within and among Vegetation Functional Groups

Abstract

Grasslands and shrublands exhibit pronounced spatial and temporal variability in structure and function with differences in phenology that can be difficult to observe. Unpiloted aerial vehicles (UAVs) can measure vegetation spectral patterns relatively cheaply and repeatably at fine spatial resolution. We tested the ability of UAVs to measure phenological variability within vegetation functional groups and to improve classification accuracy at two sites in Montana, U.S.A. We tested four flight frequencies during the growing season. Classification accuracy based on reference data increased by 5–10% between a single flight and scenarios including all conducted flights. Accuracy increased from 50.6% to 61.4% at the drier site, while at the more mesic/densely vegetated site, we found an increase of 59.0% to 64.4% between a single and multiple flights over the growing season. Peak green-up varied by 2–4 weeks within the scenes, and sparse vegetation classes had only a short detectable window of active phtosynthesis; therefore, a single flight could not capture all vegetation that was active across the growing season. The multi-temporal analyses identified differences in the seasonal timing of green-up and senescence within herbaceous and sagebrush classes. Multiple UAV measurements can identify the fine-scale phenological variability in complex mixed grass/shrub vegetation.

1. Introduction

Rangelands are widely distributed globally and contribute significant ecosystem services [1]. The grasslands, shrublands, and other ecosystems that make up rangelands are comprised of diverse plant species with divergent life-history strategies (e.g., annual vs. perennial), structural differences (e.g., grass vs. shrub), and photosynthetic pathways (e.g., C3 vs. C4), which creates pronounced spatial variability in plant function and high biodiversity [2,3,4]. Interannual variability in climate, coupled with community and abiotic differences, can lead to large interannual variability in rangeland vegetation function [5,6,7,8]. Due to the collective variability and complexity of vegetation communities, understanding of the drivers of rangeland vegetation processes (e.g., carbon accumulation, reproduction, productivity, nutrient availability, etc.) across scales is needed to develop effective management strategies and predictive models of changes in a dynamic world [9,10]. Ground-based measurements provide crucial information for understanding these processes but are limited in scope. Linking in situ measurements with remotely sensed data can expand the geographic extent of observation, furthering the ability to prioritize management efforts and aid in answering key ecological questions [11,12,13].

Data collected at high temporal frequencies and fine spatial extents may be necessary to detect complex patterns and ecologically relevant processes [14,15]. Multi-temporal analyses based on remote sensing are increasingly used to monitor ecosystems and their dynamics over multiple time steps and improve the thematic resolution of remote sensing products. For example, analyses of multiple remotely sensed images from within the same year are used to improve classification accuracy [16,17,18], measure species composition and coverage in shrublands [19], classify crop types [20], examine yields and performance [21,22], differentiate tree species in an urban environment [23], identify plants with C3 versus C4 photosynthetic pathways [24], detect invasive species [25], and examine the number of annual green-up cycles across ecosystems [26]. Importantly, examining continuous patterns can reveal new dynamics of a system across seasons and over multiple years [26,27,28].

Understanding the timing, magnitude, and duration of life-history events, phenology, forms a key aspect of the assessment of the function of vegetation communities [29]. Changes in phenology can lead to ecological disruptions, plant-herbivore mismatch, alterations to competitive interactions, and disrupted nutrient fluxes, e.g., [30,31,32,33]. Phenology and production of rangelands assessed at the community level has identified important spatial and temporal variability and trends [34,35,36,37]. However, differences in growing season length and timing exist among and within vegetation functional groups: aggregations of species with similar characteristics and ecosystem roles. Measuring heterogeneity at smaller spatial extents provides the ability to observe differences in species and functional groups such as temporal patterns in green-up, peak production, and senescence [38]. Measurement at the functional group level is needed for quantifying phenological mismatches, timing management and restoration actions, understanding the species-specific impacts of climate change, and connecting phenological processes between scales [32,38,39,40]. However, the spatial and temporal resolution needed to identify phenological events at smaller spatial scales, such as within species and functional groups, precludes many satellite-based remote sensing approaches [41,42]. In situ approaches can be time-consuming and logistically difficult across large spatial areas, which compromises our ability to observe ecosystem dynamics at all relevant scales in space and time [29].

The use of an unpiloted aerial vehicle (UAV) can examine such questions at appropriate spatial resolutions, link to aerial- and satellite-based systems [23,43,44,45], and build upon ground-based surveys [46,47,48,49]. Multi-flight UAV approaches can measure seasonal changes, e.g., [22,50,51], variability between individual plants and/or species [52,53], and can increase accuracy of classifications [54]. However, UAV approaches are not a panacea; there are multiple challenges and questions about approaches to be resolved. Phenology can be difficult to measure in dryland systems [55]. The sensors on UAVs, normally utilizing visible color and near-infrared (NIR) bands, work best for classifying dominant species, whereas rare species and some herbaceous species can be hard to identify [44], requiring additional specific data collection [51]. While using multi-temporal images can improve classification, matching scenes between UAV flights can be challenging, e.g., [42]. Furthermore, flight times and sensor payloads are limited on lightweight UAVs [43,56], constraining the size of sampling areas. Aerial systems (planes and helicopters) are used to map large extents of vegetation at fine spatial scales [45,57,58], measure canopy structure [59,60,61], and survey for wildlife [62]. While aerial platforms can carry larger payloads for longer time periods, they are more limited in takeoff and landing locations and come with increased costs [56]. Appropriate application of UAVs should consider survey area, the need for repeat flights, and research questions [43]. For example, UAVs are demonstrated to be particularly well suited to sampling ecosystem indicators in rangeland settings, e.g., [44,46,51].

Understanding phenology implies measuring vegetation at multiple intervals, which of course, incurs a cost. It is important to understand what inference can be gained (or lost) by measuring more (or fewer) times within a growing season. Due to onboard GPS accuracy limitations, approaches such as manual co-registration to known surface models [52], ground control points (GCPs) [22,51,63,64], photogrammetric image orientation [65,66], and identifying objects post-flight [42] are used to successfully align imagery from multiple UAV flights. Object-based and/or supervised image classification techniques demonstrate utility for high-resolution imagery, e.g., [54,67,68]. However, the lack of prior knowledge of classes within a scene can impact accuracy [69] and complicate classification decisions. Furthermore, some rangeland species (e.g., sagebrush, Artemisia tridentata) show intra-canopy variability from different growing season lengths for new versus prior season perennial leaves and semi-deciduous phenology [70,71]. However, while supervised classification techniques are often used for increased accuracy and a better understanding of classes, e.g., [72], unsupervised approaches identify spectral variation within a dataset [73] and may be better suited for detecting the phenological heterogeneity that exists among and within vegetation functional groups.

Therefore, a deeper investigation into the ability to use UAVs to identify within functional group phenological variability is a needed step in linking information across scales, approaches, and components of rangeland systems. Specifically, identification within-group phenological variability at fine spatial scales is a key component for monitoring and downscaling vegetation relationships to climate, providing scale-appropriate data on management, and studying ecological consequences of phenological change, e.g., [29,38,52]. The goal of this study is to examine tradeoffs in classification depth (e.g., identification of phenological differences) and accuracy between collecting data from individual versus multiple UAV flights during the growing season in rangelands. Multiple UAV flights were used to: (1) compare vegetation classification based on a single flight versus multiple flights, (2) examine the timing variations of green-up, peak productivity, and senescence within vegetation functional groups, (3) determine the increased information gained from various, multiple flight scenarios over a growing season, and (4) quantify the increased logistical and analytical effort of various flight scenarios. The prediction was that multiple flights over the growing season would increase classification accuracy due to detection of differences in phenology among vegetation functional groups and by accounting for spatial variability within functional groups, thus enabling better spectral differentiation of classes. In addition, accuracy improvement could be achieved with flights over the first half of the growing season, which would capture green-up variation and initial differences in senescence that would improve the differentiation of functional groups.

2. Materials and Methods

2.1. Study Areas

We selected two rangeland sites in southwestern Montana (Mont.), U.S.A., to capture two levels of precipitation of shrub-grassland sites in the region (Figure 1). Argenta, a drier and lower-elevation site, is located 20 km west of Dillon, Mont., in upland sagebrush steppe with primarily alluvial alfisols, sandy-loam soils, at approximately 1900 m elevation [74]. Mean monthly temperatures range from −7.3 to 16.4 °C with annual average precipitation of 283 mm (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created July 2012, accessed on 12 December 2019) [75]. The wetter, cooler, and higher elevation Virginia City site is located 5 km east of Virginia City, Mont., at an elevation of approximately 2240 m. Mean monthly temperatures range from −6.9 to 14.6 °C with annual average precipitation of 588 mm (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created July 2012, accessed on 12 December 2019) [75]. Soils are primarily mollisols with a mix of gravelly loam and stony loam on the steeper slopes [74]. Both sites are composed of sagebrush steppe; however, Virginia City generally has taller sagebrush, denser grasses, and less bare ground than Argenta.

2.2. Data Collection and Processing

We conducted numerous UAV flights (

Table 1. Information for the unpiloted aerial vehicle (UAV) flight dates and for which flights were used in the four classification scenarios at Argenta and Virginia City sites in 2018. Flights with dots indicate good data quality, and the number and color of dots indicate for which classification scenario(s) these data were used. Light green dots –all flight classification; orange dots –limited classification; purple dots –spring growing season classification; dark green dots—peak growing season, single flight UAV classification; grey bar—flights conducted but poor data quality; black bar—no flights conducted.

Flight	1	2	3	4	5	6	7	8	9
Site/Date	May-2	May-30	June-12	June-27	July-19	Aug-4	Aug-19	Sep-10	Oct-1
Argenta
Virginia City

), each consisting of multiple missions per site (individual takeoff and landing to replace batteries, change camera settings, etc.), over the 2018 growing season from early May until October at each of the study sites. We flew a solo quadcopter (3DR, Berkeley, CA, USA) with a RedEdge-M (MicaSense, Seattle, WA, USA) multispectral camera, which captures five-band images (one band per lens corresponding to the following wavelengths, blue: 475 nm with 20 nm width, green: 560 nm with 20 nm width, red: 668 nm with 10 nm width, NIR: 840 nm with 40 nm width, and red edge 717 nm with 10 nm width). We programmed our UAV missions in Mission Planner (1.3.52, Ardupilot.org, accessed on 1 December 2021) to cover an approximately 300 × 300 m scene at each site. We recorded images of a calibrated reflectance panel (provided with the camera from MicaSense) before and after each mission to correct for illumination during processing. We conducted flights using a flight altitude of 43 m and airspeed of 5 m/s with the camera programmed to capture pictures every 2 s. This produced images with 65% forelap and sidelap and a roughly 3 cm ground sample distance (pixels). In addition, the RedEdge-M camera included an upward-facing light sensor that collected corresponding illumination data with each set of five-band images. We distributed GCPs at five marked locations across our scene to provide spatial referencing between flight dates.

At Argenta, we collected data over nine flights, covering the entire growing season (Table 1). Due to weather and snow cover, we were unable to collect early season data at Virginia City for flights 1 and 2 but were able to collect data for flights 3 through 9. However, a camera malfunction during flight 4 at Argenta and flight 7 at Virginia City resulted in missing data for portions of these scenes. Therefore, all multi-flight UAV classifications described hereafter are based on up to eight flights at Argenta and six at Virginia City (Table 1).

For each flight, we used MicaSense Python tools (https://github.com/micasense/imageprocessing, accessed on 1 December 2021) to convert images to aligned, corrected five-band stacks. This process included us using the calibration panel images to convert radiance values to reflectance values, correct for camera lens effects, and align each individual band together into a single image. We then inputted these images into Metashape (v1.5.2, Agisoft, St. Petersburg, Russia) for photogrammetric processing. Although the flights at each site followed the same flight plan, minor differences in extent and resolutions existed between the orthomosaics; therefore, we exported all orthomosaics for each site to the smallest extent and largest resolution identified across all flights to facilitate additional raster calculations. We transferred Metashape-produced products (e.g., orthomosaics, Digital Elevation Models (DEMs)) into ArcGIS 10.6 (ESRI, Redlands, CA) for additional processing steps as detailed below.

We employed a four-dimensional structure-from-motion approach (4D SfM) for our multi-date imagery [65,66] and used ground targets to assess the performance of image and flight alignment. The 4D SfM technique works by including survey photographs from all flights (i.e., dates) in the alignment, camera calibration, and georeferencing of ground control steps of the photogrammetry process. The process located objects (e.g., plants, rocks, litter, etc.) that remained in fixed locations through the entire time series. Effectively, 4D SfM identified thousands of ground control points that are consistent through the temporal period of the imagery to reference all photographs prior to building orthomosaics and other output products [65,66]. Following the 4D SfM alignment, we separated images by date and site to complete the remaining photogrammetry steps (generation of point cloud, surface models, orthomosaics, etc.). This process resulted in single five-band orthomosaics for each flight date: eight for Argenta and six for Virginia City. We calculated the normalized difference vegetation index (NDVI) for each UAV orthomosaic in ArcGIS using Image Analyst with associated red and NIR bands.

Active remote sensing sensors and structure-from-motion techniques utilized in aerial image processing produce three-dimensional models of the target scene, which allows for structural representations of vegetation [59,61,76]. The use of UAVs is expanding for mapping forest canopy structure, e.g., [60,76], and more recently, these techniques were applied to other vegetation types, e.g., [51,68,77,78]. In rangeland environments, shrub and grass structures may help in differentiating vegetation types [68]; therefore, we produced a canopy height model for our study sites. We estimated the vegetation height at each site from the UAV flight at the peak of the growing season (Table 1). We used the automatic Classify Ground Points function in Metashape, which classifies points within the dense point cloud as ground or other (i.e., vegetation) using three user-defined inputs. We explored multiple combinations of inputs (max angle, max distance, and cell size) to balance the ability of the classification algorithm to correctly assign the full canopy of shrubs while not missing a significant portion of ground points (e.g., being misclassified as vegetation). We created a digital surface model (DSM) from all points in the dense cloud and a digital terrain model (DTM), which uses only ground points and interpolates the elevation beneath vegetation. We calculated vegetation height using the raster calculator in ArcGIS by using the difference in DEMs technique, effectively subtracting the DTM from the DSM.

To provide additional information for classification, we created a texture raster layer based on the peak growing season flight for each site. We first created a single band image using the grayscale function from Image Analyst in ArcGIS. We equally weighted the red, green, and NIR bands from the peak growing season flight, effectively creating a grayscale version of a color infrared image. Then, using the focal statistics tool from Spatial Analyst, we calculated the standard deviation of a circular five-pixel region of the color infrared grayscale (single band) image. Window size and shape were chosen after testing other buffer distances for edge and scale effects [79,80]. Vegetation height and texture were normalized to match the range of NDVI values.

2.3. Image Classification

We classified vegetation for four sets of UAV flights (Table 1) to compare phenological information and classification accuracy differences between four scenarios of flight timing and frequency. These scenarios include (1) a single flight at the presumed peak of the growing season (hereafter “single”), (2) three flights over the growing season consisting of an early, peak, and end of season flight (hereafter “limited”), (3) flights during the early and core part of the growing season to cover initial green-up and the beginning of senescence (hereafter “spring”), and (4) using all available data covering the full growing season (hereafter “all flights”). We performed an iterative self-organized (ISO) unsupervised data analysis classification using NDVI from the selected flight(s) and vegetation height and texture data derived from the flight at the peak of the growing season. We allowed for a maximum of 30 classes. We repeated this approach for the multiple flight scenarios but included additional NDVI values with the same normalized vegetation height and texture. For Argenta, the all flights classification was based on 10 inputs (8 flights), the spring classification on 7 inputs (5 flights), and the limited classification on 5 inputs (3 flights). For Virginia City, the all flights classification was based on eight inputs (six flights), the spring classification on six inputs (four flights), and the limited classification on five inputs (three flights).

For each scenario, we examined the mean height and NDVI of the resulting classes along with visual interpretation to assign classes and collapsed classes together that could not be categorized into separate vegetation types, densities, or phenological patterns. For the single flight scenario, we used a single mean NDVI value; however, we used the mean NDVI values across included flight dates (i.e., a partial phenological curve) for each included flight for multiple flight scenarios. We first classified the image to identify a base set of six functional groups: bare ground, litter, sparse, medium, and dense herbaceous, and sagebrush. Differentiation between sparse, medium, and dense herbaceous categories was based on mean NDVI of the class (NDVI or peak NDVI for sparse between 0.2 and 0.3, medium between 0.3 and 0.4, and dense > 0.4), supplemented by visual inspection of pixels and areas of the class. Sagebrush classes also had an NDVI or peak NDVI of about 0.45, but vegetation height was higher (mean about 0.3 m versus < 0.1 m for herbaceous). Soil and litter classes were identified based on NDVI patterns that showed only minor deviations over the growing season (due to changes in moisture content and possible pixel mixing), height (for standing litter), and visual inspection of color and/or location in the scene. In addition, if possible for a given scenario, we identified subcategories of each of these classes to reflect differences in green-up and senescence timing and height.

We chose to use NDVI as our primary input into our classification algorithms as it is a surrogate for plant photosynthetic activity and was an indicator that matched our objective of examining phenological differences. In addition, even after post-processing, we had illumination differences due to clouds and shadows among flights which caused discontinuities on the five original bands produced by the MicaSense RedEdge-M camera that can impact the classification algorithms. Illumination differences were greatly reduced or eliminated by employing a ratio of bands such as NDVI, which removed brightness issues within and across flights.

2.4. Accuracy Assessment

To assess the accuracy of our classification scenarios, we performed analyst-based classification to derive reference data following recommendations for using the same orthomosaics as used for our classification, see [81,82]. At each study site, we created a stratified random sample of 1000 points, stratified on the classes derived from the all flights scenario. We used zonal statistics to calculate the mean NDVI for all flight classes for each flight date to create NDVI curves representing the phenological pattern for each class. Then, for each sample point, we extracted NDVI values at each flight date and vegetation height as derived from the peak growing season flight. Then, at each sample point, we deduced the reference classification for the point by comparing NDVI curves, height ranges, and RGB images from available flights. These reference data were then used to create a confusion matrix for each classification scenario at each site. Although reference data were based on all flight classifications (with subcategories), we collapsed subcategories into the main class to allow comparison across the full set of scenarios. In addition, we used our base set of six functional groups (bare ground, litter, sparse, medium, and dense herbaceous, and sagebrush) as a standard comparison of scenarios, as the single flight classification approach contains no subcategories. For all confusion matrices, we also calculated the kappa statistic (the classification agreement relative to random allocation) as a measure of accuracy [[73], but see [81] for limitations].

3. Results

3.1. Flight Classifications

The single flight classifications based on the 12 June 2018, flight at Argenta and the 27 June 2018, flight at Virginia City resulted in our base set of six functional groups. Most multi-flight scenarios could identify different phenological patterns within functional groups, i.e., subcategories (Figure 2, see also Appendix A). Subcategories were defined in several ways because we found that areas within the same functional group were actively growing (or peaking) at different time periods. Subcategories were defined by differences in the length of the growing season (i.e., how long pixels had NDVI values above baseline), the timing of green-up (NDVI increase), peak values (day of peak NDVI), and senescence (when NDVI begins to quickly decline). Furthermore, at the Virginia City site, sagebrush heights covered a larger range, and we could differentiate sagebrush into two height classes that also coincided with the magnitude of NDVI values. Including subcategories, we differentiated 9 vegetation classes at Argenta (Figure 2, Figure 3 and Figure A1) and 10 classes at Virginia City (Figure 2, Figure 4 and Figure A2), using the all flights scenario for each site. In the sparse flights scenario, we were unable to identify subcategories at Argenta, which resulted in the 6 main classes (Figure A3); however, we identified the same 10 classes as in the all flight scenario at Virginia City (Figure A4). For the spring scenario, we could partially differentiate the 6 classes at Argenta, with differences observed in the herbaceous class but not in the sagebrush class (Figure A5), and fully differentiate all 10 classes at Virginia City (Figure A6).

Subcategories were identified within the herbaceous and sagebrush classes and represent different phenological responses within these functional groups (Figure 2, Figure 3 and Figure 4). At Argenta, NDVI slowly increased between flights 1 and 2, rapidly increased to a peak at flight 3, and then declined, with differences within herbaceous functional groups related to NDVI decline. The subclass ‘short season’ quickly declined after the peak at flight 3, whereas the moderate season remained close to peak values at flight 4, then declined quickly after flight 5. Flights 3 and 4 were 15 days apart, so these phenological differences represent at least two additional weeks of increased photosynthesis. Sagebrush short and moderate seasons were similar to herbaceous classes but shifted one additional week into the growing season (short season declines after flight 4, moderate after flight 5).

At Virginia City, senescence dates across classes were similar, and herbaceous subcategories were differentiated by early or late green-up dates. Early season herbaceous had a peak value at flight 3 (the first available at Virginia City) and slowly declined through flight 7, whereas late-season continued to increase through flights 4 and 5 before declining. Flights 3 and 5 were 37 days apart, representing a peak green-up value about a month later for the late-season phenological subcategories. Differences in sagebrush subcategories were related to height, with a short class (mean height about 0.1 m, primarily in the upslope portions of the scene) and a tall subcategory (mean height about 0.25 m). In addition, there was one class where pixels were mixed between sagebrush and herbaceous vegetation with an early peak green-up at flight 3 (from herbaceous vegetation) and a long season with plateaued NDVI values through flight 5 and a drop after, similar to sagebrush, and with height values (mean height of 0.24 m), similar to the tall class.

3.2. Accuracy, Class Differentiations, and Comparisons between Scenarios

To compare the four scenarios, we initially assessed the accuracy of each scenario in identifying the six base land cover classes (

Table 2. Summary of accuracy results from confusion matrices derived from vegetation classification scenarios at sites near Virginia City and Argenta, Montana. Base categories include bare ground, litter, sparse, medium, and dense herbaceous, and sagebrush. Subcategories include differences between green-up and senescence timing and/or height within a base class.

	Scenario
Site	Single		Limited		Spring		All
	Overall	Kappa	Overall	Kappa	Overall	Kappa	Overall	Kappa
Argenta
Base Categories	50.6%	0.50	51.6%	0.51	52.5%	0.52	61.4%	0.61
Subcategories	---	---	---	---	45.6%	0.46	51.8%	0.52
Virginia City
Base Categories	59.0%	0.59	60.2%	0.60	61.6%	0.62	64.4%	0.64
Subcategories	---	---	44.9%	0.39	46.6%	0.40	53.2%	0.53

,

Table A1. Confusion matrix of vegetation classification accuracy from a single unpiloted aerial vehicle (UAV) flight on June 12th, 2018, at the Argenta site. Kappa = 0.50.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	105	39	10	5	1	2	162	64.8%
	Litter	3	14	3	0	0	7	27	51.9%
	Sparse Herb	29	51	92	14	4	12	202	45.5%
	Medium Herb	1	2	52	100	20	50	225	44.4%
	Dense Herb	0	2	10	62	18	76	168	10.7%
	Sagebrush	0	1	9	26	3	176	215	81.9%
	Total	138	109	176	207	46	323	999
	Producers Accuracy	76.1%	12.8%	52.3%	48.3%	39.1%	54.5%		50.6%

,

Table A2. Confusion matrix of vegetation classification accuracy from three unpiloted aerial vehicle (UAV) flights (the limited flights scenario) in 2018 at the Argenta site. Kappa = 0.52.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	71	21	7	3	0	1	103	68.9%
	Litter	29	37	30	12	2	8	118	31.4%
	Sparse Herb	36	46	97	44	8	45	276	35.1%
	Medium Herb	1	1	32	99	24	43	200	49.5%
	Dense Herb	0	1	1	24	8	19	53	15.1%
	Sagebrush	1	3	9	25	4	207	249	83.3%
	Total	138	109	176	207	46	323	999
	Producers Accuracy	51.5%	33.9%	55.1%	47.8%	17.4%	64.1%		51.9%

,

Table A3. Confusion matrix of vegetation classification accuracy from 4 unpiloted aerial vehicle (UAV) flights (spring flights scenario) in 2018 at the Argenta site. Kappa = 0.52.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	106	28	33	6	1	2	176	60.2%
	Litter	19	38	25	13	0	6	101	37.6%
	Sparse Herb	10	33	90	65	2	9	209	43.1%
	Medium Herb	2	8	25	113	27	132	307	36.8%
	Dense Herb	0	0	0	0	6	6	12	50.0%
	Sagebrush	1	2	3	10	10	168	194	86.6%
	Total	138	109	176	207	46	323	999
	Producers Accuracy	76.8%	34.9%	51.1%	54.6%	13.0%	52.0%		52.2%

,

Table A4. Confusion matrix of vegetation classification accuracy from 8 unpiloted aerial vehicle (UAV) flights in 2018 (all flights scenario) at the Argenta site. Kappa = 0.61.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	83	15	8	3	0	2	111	74.8%
	Litter	24	41	21	10	1	14	111	36.9%
	Sparse Herb	29	47	110	28	1	7	222	49.6%
	Medium Herb	0	4	34	135	11	38	222	60.8%
	Dense Herb	1	1	2	29	30	48	111	27.0%
	Sagebrush	1	1	1	2	3	214	222	96.4%
	Total	138	109	176	207	46	323	999
	Producers Accuracy	60.1%	37.6%	62.5%	65.2%	65.2%	66.3%		61.4%

,

Table A5. Confusion matrix of vegetation classification accuracy from a single unpiloted aerial vehicle (UAV) flight on June 27th, 2018, at the Virginia City site. Kappa = 0.59.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	66	22	37	2	0	0	127	51.9%
	Litter	21	45	51	25	0	24	166	27.1%
	Sparse Herb	0	2	14	25	5	15	61	22.9%
	Medium Herb	0	0	0	43	36	6	85	50.6%
	Dense Herb	1	0	0	3	141	38	183	77.1%
	Sagebrush	2	16	6	52	21	281	378	74.3%
	Total	90	85	108	150	203	364	1000
	Producers Accuracy	73.3%	52.9%	12.9%	28.7%	69.5%	77.2%		59.0%

,

Table A6. Confusion matrix of vegetation classification accuracy from three unpiloted aerial vehicle (UAV) flights in 2018 (limited flights scenario) at the Virginia City site. Kappa = 0.60.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	78	45	25	2	0	7	157	49.7%
	Litter	4	9	7	11	0	26	57	15.8%
	Sparse Herb	5	19	74	26	0	9	133	55.6%
	Medium Herb	0	3	1	65	58	37	164	39.6%
	Dense Herb	0	0	0	2	104	13	119	87.4%
	Sagebrush	3	9	1	44	41	272	370	73.5%
	Total	90	85	108	150	203	364	1000
	Producers Accuracy	86.7%	10.6%	68.5%	43.3%	51.2%	74.7%		60.2%

,

Table A7. Confusion matrix of vegetation classification accuracy from five unpiloted aerial vehicle (UAV) flights in 2018 (spring flights scenario) at the Virginia City site. Kappa = 0.615.

		Ground Truth
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	78	34	33	2	0	2	149	52.4%
	Litter	7	25	23	8	0	19	82	30.5%
	Sparse Herb	1	9	45	21	0	11	87	51.7%
	Medium Herb	0	2	0	54	26	43	125	43.2%
	Dense Herb	0	0	0	7	139	14	160	86.9%
	Sagebrush	4	15	7	58	38	275	397	69.3%
	Total	90	85	108	150	203	364	1000
	Producers Accuracy	86.7%	29.4%	41.7%	36.0%	68.5%	75.6%		61.6%

and

Table A8. Confusion matrix of vegetation classification accuracy from six unpiloted aerial vehicle (UAV) flights in 2018 (all flights scenario) at the Virginia City site. Kappa = 0.64.

		Reference
	Class Names	Bare Ground	Litter	Sparse Herb	Medium Herb	Dense Herb	Sagebrush	Total	User’s Accuracy
Classification	Bare Ground	73	19	6	1	0	1	100	73.0%
	Litter	7	38	33	6	0	16	100	38.0%
	Sparse Herb	6	16	57	18	0	3	100	57.0%
	Medium Herb	2	5	12	93	17	71	200	46.5%
	Dense Herb	0	0	0	16	147	37	200	73.5%
	Sagebrush	2	7	0	16	39	236	300	78.7%
	Total	90	85	108	150	203	364	1000
	Producers Accuracy	81.1%	44.7%	52.8%	62.0%	72.4%	64.8%		64.40%

). The accuracy from the single flight classifications at the Argenta site (50.6%) was marginally improved under the limited flights scenario (51.9%) and spring scenario (52.2%). However, the all flight scenario improved by almost 10% (61.4%). The accuracy of the all flight scenario, when including subcategories (the full nine identified including phenological differences), was similar to the single flight and other multiple flight scenarios (

Table A9. Confusion matrix of vegetation classification accuracy from eight unpiloted aerial vehicle (UAV) flights (all flights scenario) in 2018 at the Argenta site. Kappa = 0.52. S—short duration growing season, and M—moderate duration growing season.

		Reference
	Class Value	Bare Ground	Litter	Sparse Herb (S)	Sparse Herb (M)	Medium Herb (S)	Medium Herb (M)	Dense Herb	Sagebrush (S)	Sagebrush (M)	Total	User’s Accuracy
Classification	Bare Ground	83	15	3	5	0	3	0	2	0	111	74.8%
	Litter	24	41	8	13	2	8	1	7	7	111	36.9%
	Sparse Herb (S)	11	15	50	14	12	6	0	3	0	111	45.0%
	Sparse Herb (M)	18	32	18	28	3	7	1	1	3	111	25.2%
	Medium Herb (S)	0	2	7	3	72	7	4	13	3	111	64.9%
	Medium Herb (M)	0	2	7	17	9	47	7	9	13	111	42.3%
	Dense Herb	1	1	1	1	11	18	30	26	22	111	27.0%
	Sagebrush (S)	0	1	0	0	1	1	2	85	21	111	76.6%
	Sagebrush (M)	1	0	0	1	0	0	1	27	81	111	73.0%
	Total	138	109	94	82	110	97	46	173	150	999
	Producers Accuracy	60.1%	37.6%	53.2%	34.1%	65.5%	48.5%	65.2%	49.1%	54.0%		51.8%

, overall accuracy = 51.8%). We were unable to identify any phenological differences (i.e., no subcategories) under the limited flights scenario. However, we were able to identify seasonal differences within the herbaceous categories with the spring scenario but with limited accuracy (45.6%).

At the Virginia City site, a similar class comparison (with six classes) resulted in an overall accuracy of 59.0% for the single flight scenario, 60.2% for the limited scenario, 61.6% for the spring scenario, and 64.4% for the all flight scenario (Table 2). When considering phenological differences, we were able to identify the same 10 classes in each of the multiple flight scenarios, but the overall accuracy of the limited (38.5%) and spring (40.0%) were poor compared to the all flight scenario (53.2%,

Table A10. Confusion matrix of vegetation classification accuracy from six unpiloted aerial vehicle (UAV) flights (all flights scenario) in 2018 at the Virginia City site. Kappa = 0.53. Herb—herbaceous, E—early green-up, L—later green-up, U—short (generally upslope), T—tall, Mix—mixed herbaceous and sagebrush.

		Reference
	Class Value	Bare Ground	Litter	Sparse Herb (E)	Medium Herb (S)	Medium Herb (L)	Dense Herb (S)	Dense Herb (L)	Sagebrush (U)	Sagebrush (Mix)	Sagebrush (T)	Total	User’s Accuracy (%)
Classification	Bare Ground	73	19	6	0	1	0	0	1	0	0	100	73.0
	Litter	7	38	33	4	2	0	0	12	3	1	100	38.0
	Sparse Herb	6	16	57	16	2	0	0	1	2	0	100	57.0
	Medium Herb (E)	1	1	9	43	11	8	1	13	10	3	100	43.0
	Medium Herb (L)	1	4	3	6	33	3	5	34	0	11	100	33.0
	Dense Herb (E)	0	0	0	7	0	75	16	1	0	1	100	75.0
	Dense Herb (L)	0	0	0	0	9	4	52	9	2	24	100	52.0
	Sagebrush (U)	0	1	0	1	10	4	2	42	4	36	100	42.0
	Sagebrush (Mix)	0	1	0	0	1	11	8	3	57	19	100	57.0
	Sagebrush (T)	2	5	0	0	4	5	9	9	4	62	100	62.0
	Total	90	85	108	77	73	110	93	125	82	157	1000
	Producers Accuracy (%)	81.1	45.7	52.8	55.8	45.2	68.2	55.9	33.6	69.5	39.5		53.2

).

Typically, the lowest class-specific accuracy rates at both sites were for the litter class across all scenarios (Table A1, Table A2, Table A3, Table A4, Table A5, Table A6, Table A7, Table A8, Table A9 and Table A10), with improvement in most situations between the single and multi-flight approaches (Table 2). The misclassification of litter was primarily with bare ground and sparse herbaceous classes. Herbaceous misclassification errors tended to be between densities (sparse, medium, and dense), in addition to some errors between dense herbaceous and sagebrush misclassifications at Argenta and to a lesser extent at Virginia City. The identification of phenological subcategories resulted in increased classification errors within the multi-flight classifications. These errors tended to be between the phenological classes (e.g., short versus moderate season length herbaceous); however, at Argenta, we also found classification errors between herbaceous and sagebrush classes. More flights (i.e., moving from single, to limited, to spring, to the all flights scenario) improved the classification of subcategories within functional groups (Table 2, Table A9 and Table A10).

4. Discussion

This study examined the ability of UAVs to identify fine-scale phenological heterogeneity, can be used to inform logistical decisions of future UAV studies (timing, number of flights, etc.), and illustrated options to effectively process data from multiple UAV flights. This approach successfully identified plant level phenological differences (i.e., differences in growing season length and timing within vegetation functional groups) and demonstrated improved classification accuracy by utilizing multi-flight UAV classification approaches. At both study locations, multi-flight scenarios improved vegetation classification accuracy over a single flight, mainly due to asynchronies in the timing and duration of herbaceous cover. Two main findings emerge from this study: (1) phenological heterogeneity at fine spatial scales can be identified by UAV flights, and (2) heterogeneity in phenological timing causes accuracy issues between class designations. The implications of these results include the ability to use these data to make finer scale ecological comparisons than satellite-based land surface phenology measures. In addition, while land surface phenology measures from satellite remote sensing can be influenced by the proportion of vegetation functional group, which may exhibit different phenological patterns, within a pixel [83], this study’s results suggest that within functional group heterogeneity needs to be considered as well. Furthermore, UAV-based studies provide phenological information at the plant or sub-plant level to compliment the broader spatial coverage available from satellite-based systems [22,42,52].

4.1. Phenological Heterogeneity within Functional Groups

Remote sensing is increasingly used to measure land surface phenology, which incorporates the mixed response of vegetation and background materials (e.g., soil and litter) within a sensor’s spatial resolution [83,84,85]. This study demonstrated the ability to examine patterns within functional groups at a fine resolution. Many of our multi-flight classification scenarios based on NDVI at multiple points during the growing season identified spatial variation within functional groups. These differences highlighted the role of topographic variations in phenology. For example, small valleys at Argenta have a longer growing season than upland locations, and there are seasonal asynchronies in herbaceous species between south-facing slopes and other areas at Virginia City. A study in marsh vegetation also found phenological differences in timing and duration within close spatial proximity, with implications for restoration and management [86]. Additional asynchronies in plant responses identified from remote sensing include differences in photosynthetic pathways [24], water use efficiency [28], growing conditions within fields [22], and variable intra- and inter-specific phenology [52]. Changing or variable vegetation phenology has consequences to other species and processes within the ecosystem [30,31,32,33,87]. Fine-scale measurement of variable vegetation functional group phenology is needed to inform restoration planning [86] and to better quantify the degree of phenological mismatch between members of an ecosystem under a changing climate. Understanding plant-level phenological differences as demonstrated here helps advance understanding of these processes.

At the drier Argenta site, there were functional group differences when plants were senescing and differences in the timing of green-up at the more mesic Virginia City site. These differences in limiting factors between sites may be difficult to identify a priori. In other sagebrush systems, communities in meadow (more mesic) locations had longer growing seasons compared to upland communities, with the difference due to earlier start of season dates [70]. In addition, new growth on sagebrush has a longer growing season [70], matching some of the sub-canopy variation identified herein, possibly also tied to the semi-deciduous pattern of sagebrush [71]. Phenology differences between and within species can be highly variable [23,52], asynchronies can be in the spring [70] or fall [42], and can be driven by variable growing conditions [22,70]. Building on phenology, temporal rangeland assessments can aid in ecological and land use studies. For example, high spatial resolution imagery can be used for within-season forage utilization [51], monitoring vegetation trends over time [44], and identifying invasive species [25].

4.2. Accuracy and Tradeoffs of Single versus Multi-Flight Approaches

Utilizing multiple flights over the growing season increased the accuracy of our classifications by up to 5–10% when comparing the different multi-flight scenarios to the single-flight scenarios (Table 2). In a single-flight classification, shadows are difficult to classify accurately [88], and spectral differences through time cannot be utilized. Multi-temporal UAV imagery improved accuracy in other cases, resolving issues with flowering and shadows, as well as utilizing phenology to identify different spectral signatures over time [25,54]. While our overall accuracy was higher at the more densely vegetated site (mesic), the greatest increase in accuracy occurred at the sparser, drier site. The improved accuracy at both sites was from increasing the correct classifications of short duration herbaceous, bare ground, and litter, with these classes covering more of the land surface at Argenta. Specifically, the spectral patterns of soil and litter pixels separate out over the growing season. Sparse, short-duration herbaceous material may be missed with one flight, even at the presumed peak of the growing season, as it is spectrally similar to litter or soil except for a short time period, which may differ temporally across a scene. This study only included two sites; therefore, to further test the hypothesis of increased accuracy at drier sites but better overall accuracy at wetter sites, a study examining accuracy across a precipitation gradient would be advantageous.

Even with multiple flights, classification challenges remain between spectrally similar classes, such as bare ground and sparse vegetation [17]. Likewise, the phenological cycle is very similar between dense herbaceous and sagebrush groups (Figure 3 and Figure 4), which also had overlapping height distributions, leading to some misclassification errors between these categories. Overall, classification errors were more likely between similar classes or subcategories (e.g., bare ground vs. litter, the density of herbaceous classes, and length of the growing season) than distinct classes (e.g., herbaceous vs. bare ground or sagebrush vs. litter). However, due to scene heterogeneity, fuzzy class boundaries, and possible mixed pixels, these errors exist between distinct classes in the classification. Classification of our rangeland sites, even with high spatial resolution, still had challenges (Figure 5), and we note several specific issues for our study areas.

First, some of the hardest pixels to identify are edges (e.g., sagebrush exterior pixels), as these were often mixed vegetation types, despite our small pixel size (~3 cm). An example of mixed classes within a pixel is grass growing up through sagebrush or through sagebrush skeletons. The mixing of desired classes within a pixel results in classification challenges and complicates the creation of reference data. Second, shadows move through time-series data and complicate the classification even when using NDVI to limit these effects. Third, the shapes of the phenology curves are highly variable; very few actually match up with the calculated class mean curves (Figure 3 and Figure 4). Fourth, relatively small plants (or short sagebrush) can be hard to visually classify or separate between classes. The advantage of the multiple flight approach is while accuracy drops when classifying based on within functional group heterogeneity, this accuracy is the same as the single flight but provides more information about the timing and duration of when herbaceous and sagebrush classes are photosynthesizing.

We tested the limited and spring flight scenarios to assess whether reducing the frequency of flights or concentrating flights during the primary growing season would have similar accuracy and ability to identify within-group phenological heterogeneity as compared to the all flights scenario. While we predicted that we could achieve similar accuracy for the classification of both vegetation categories and subcategories containing phenological differences, we instead only found increased land-cover (at the class level) accuracy with the spring and limited scenarios. We found limited accuracy and/or ability to identify subcategory phenological differences, as the timing of the flights was key for identifying these within functional groups phenological differences. In our approach, late-season flights helped separate sagebrush (higher NDVI in the early and late season as most leaves stay on) from herbaceous classes (fully brown/senesced). Sagebrush and herbaceous (medium and dense classes) had similar green-up timing and midseason patterns, as well as overlapping heights in some cases, so differences between these classes occurred in late summer. The limited scenario worked better at the higher elevation Virginia City site indicating the timing of flights should be concentrated during the green-up period (near the peak rather than the early season when there is little activity). For a study in the nearby U.S. Great Basin, accurately separating invasive species and other classes from fine spatial resolution multi-flight UAV imagery was found to be tied to having the correct timing of flights that captured phenological differences, more than increased spectral information from any given point in time [25]. In addition, the growing season at Virginia City (with the exception of the NW part of the scene) was longer than at Argenta, increasing the window to capture flights that are diagnostic in separate classes. In differentiating between herbaceous, especially sparse, limited duration classes, and bare ground/litter classes, there was often only a short window (i.e., short green/growing season) that was asynchronous across the scene. Classification accuracy can be improved with multiple flights, and these should be focused on a combination of spring and late season flights. However, accurately differentiating subcategories required flights across the whole growing season in our study and was more successful in areas of the scene without herbaceous species mixed within the shrub canopy.

Tradeoffs between the single and multiple flight approaches beyond the type of information gained are primarily logistical and processing-time based. We reused flight plans and marked ground target locations with rebar for easy relocation. However, despite increased efficiency in operations through the growing season, multiple flights take up a significant amount of time and should be weighed against specific objectives. Specifically, in our case, we conducted 16 total flights, and with ideal flight conditions, theoretically could have mapped vegetation distribution at 16 sites instead of two with a similar effort but with lower accuracy. It took us twofour missions to cover each site based on battery life of 9–10 min, wind speed, cloud cover, and other operational factors (e.g., software issues, clear air space, etc.). Including setup/takedown, the time to swap batteries, and reacquire GPS position, our total operation took between 40 and 60 min per site, per round. Processing time also increased with multiple flights, although most of this effort was in increased background processing. After initial loading and the creation of project workspaces in Python and Metashape, we utilized batch processing. Therefore, manual input time was marginally increased, while total processing time was increased by approximately a factor of six–eight (representing the number of flights at each site). However, in many situations, these increased data collection and processing times can greatly increase the accuracy of classification approaches [17,23,54]. Specific applications and study questions will dictate required information, accuracy, precision, and sample size. Individual applications must weigh these factors and make a decision that balances economy with accuracy needs.

Multi-temporal UAV analyses can complement and build on techniques from multi-temporal satellite-based remote sensing, e.g., [14,89]. Furthermore, high-resolution imagery was used to create training data for satellite-based image classification, e.g., invasive species in [90], relatively pure signals for spectral unmixing, e.g., [91], and to validate continuous cover models, e.g., [92]. However, additional research is needed to improve accuracy and application, as well as to capture spatial and temporal variability to answer questions such as measuring the success of management actions, identifying consequences of a changing climate, and quantifying impacts to ecosystem functions and services. Other approaches such as object-based classification [54], are used in multi-temporal classification; however, the variety of scales of vegetation sizes, and desire to identify possible within canopy phenological differences in sagebrush [70], precluded use in this study. Additional improvements to classification approaches such as classification based on vegetation functional groups with continuous values for phenology (e.g., start of spring vs. early green-up), as well as object-based approaches that can capture fine scale phenological differences are needed. UAVs are a tool to sample the high spatial variability of rangeland ecosystems, and emerging areas of research and application of multi-flight UAV studies show great promise for improving monitoring and assessment of these systems.