Mapping Quaking Aspen (Populus tremuloides Michx.) Using Seasonal Sentinel-1 and Sentinel-2 Composite Imagery across the Southern Rockies, USA

1. Introduction

Across the interior western United States, quaking aspen (Populus tremuloides Michx.) is the dominant deciduous tree species in primarily mixed-conifer forests [1,2]. This important forest type provides ecosystem services, including biodiversity, wildlife habitat, recreational value, soil carbon sequestration, and vegetation recovery [1]. Aspen respond readily to canopy opening events (e.g., wildfire) and its persistence on the landscape is closely linked to disturbance type and periodicity [3,4,5,6]. Increasing synchronous disturbances (e.g., drought, bark beetle, and wildfire) may favor future dominance of quaking aspen in some regions [5,7,8]. However, a drier climate and altered disturbance regimes are likely to have varied impact on this foundational species across its range [9,10,11,12]. With growing focus on aspen ecology and management in the context of climate change and forest resiliency, production of high-resolution, accurate maps of its distribution at scales relevant to forest management activities is necessary.

Current nationwide maps of quaking aspen forest cover are developed using Landsat imagery at a 30-m spatial resolution (e.g., [13]) and often fail to identify small stands, which offer prolific root suckering and expansion when managed or following disturbance [14,15]. There are three relevant 30-m maps of aspen distribution widely available for land managers in the Southern Rockies and US: 1 – the LANDFIRE Existing Vegetation Type (Landfire EVT) [13], 2 – the United States Forest Service (USFS) National Individual Tree Species Parameter maps (USFS ITSP) [16], and 3 – the USFS TreeMap [17]. Of these, only the Landfire EVT is provided in more than one time stamp. The remaining databases provide a snapshot in time as valuable baselines for managers. Collectively, these databases are limited in accuracy, their ability to detect small patches of aspen, and in their temporal range. For example, one of the more commonly used land cover maps in the US, the Landfire EVT, has an accuracy of 60-63% within aspen classes across the southwest US region for the 2016 product (https://www.landfire.gov/remapevt_assessment.php, accessed March 31, 2023). The USFS ITSP (c. 2014) and TreeMap (c. 2016) uniquely provide additional information from the USFS Forest Inventory and Analysis (FIA) plot data such as basal area and stem density but provide only a snapshot in time and have similar moderate accuracy for the aspen class. Developing new, higher-resolution maps will improve on these existing products at scales relevant to quaking aspen management.

Remote sensing science has benefited in recent years from increasingly available free and open imagery and improved cloud computing resources, facilitating the application of large scale analysis of the Earth’s surface [18,19]. In particular, the Sentinel missions, developed to support the Copernicus Programme and administered by the European Space Agency, have improved LULC mapping efforts across broad geographic regions [20]. The Sentinel missions consist of a constellation of satellites carrying a range of technologies, including microwave (radar) and multispectral imaging. These data are made free and open to the public through the Copernicus Programme free, full, and open data policy. Importantly, Sentinel products are also hosted in the Google Earth Engine (GEE) platform, which provides access to petabytes of earth systems data co-located with powerful cloud-computing resources [19].

The combination of spatial, spectral, and temporal resolutions of the Sentinel-2 (S2) Multispectral Instrument (MSI) is particular useful for mapping vegetation with distinct phenological patterns [20]. Since 2015, S2 has provided global multispectral imagery at finer spatial (10-60 m), spectral (13-band), or temporal (5-day) resolutions compared to similar satellite missions such as Landsat (30-90 m, 4-8 bands, and 16-day) or Satellite Pour L’Observation de la Terre (1.5-6 m, 4-band, and 1-3 day). The four narrow red-edge bands captured by S2 improve species-level mapping, estimation of canopy chlorophyll content, and derivation of metrics such as gross primary productivity [21,22,23,24]. The 5-day temporal resolution increases availability of cloud-free images and facilitates analysis of land surface phenology and generation of seasonal cloud-free image composites [25]. Seasonal imagery improves species discrimination in forested regions, especially for deciduous types with distinct phenology [26,27,28,29].

However, as with any optical imagery, persistent cloud cover, haze, and snow cover are common issues affecting data acquisition and quality. In the case of vegetation mapping with multispectral imagery, multi-temporal image gap-filling can help overcome this challenge while retaining seasonal characteristics of the land surface [25]. Additionally, combining multispectral images with sensors less affected by atmospheric and environmental conditions, such as radar and microwave, can improve availability and quality of data.

The Sentinel-1 (S1) provides cloud-penetrating synthetic aperture radar (SAR) backscatter imagery at the same spatial and temporal resolution as S2, providing an important complementary data source for applications in land surface monitoring. As an active remote sensing system, SAR measures the amplitude and phase of the backscatter signal, which corresponds to the physical properties (e.g., roughness) of the land surface [30]. In forested systems, the wavelength of radar signal emitted from the sensor determines the sensitivity and depth of penetration into the canopy [31]. The S1 collects dual-polarized C-band SAR backscatter at a nominal frequency range, from 4-8 GHz (3.75-7.5 cm wavelength) in the microwave region of the electromagnetic spectrum. In this region, there is minimal penetration into the forest canopy, making it useful for characterizing tree canopy roughness and texture [31,32].

Recent studies have shown potential for multi-temporal optical and radar (e.g., S1 and S2) imagery to improve vegetation mapping, especially in areas with persistent cloud cover and heterogeneous vegetation types [33,34,35]. The common spatial and temporal resolutions of the S1 and S2 facilitate the use of these two data sources in tandem for classifying vegetation type [36,37]. Radar imagery is particularly useful in characterizing deciduous forest canopy structure during winter months when there is often persistent snow cover and/or cloud cover [33,34]. In deciduous forests specifically, the canopy structure changes at key phenological stages (e.g., onset of greenness increase, greenness maximum, onset of senescence, and greenness minimum), which may be characterized by multi-temporal S1textural data, improving forest species classifications [25]. Furthermore, because SAR is sensitive to changes in the canopy roughness, textural features derived from backscatter imagery have improved classification of forest type [32,38]. Given the advances in classification accuracy using the Sentinel suite in a variety of recent studies, the efficacy in different ecosystems should be explored.

In this study, we have three primary objectives: 1) develop methods for an open-source, accurate, and high-resolution (10 m) map of aspen cover across the Southern Rockies ecoregion (U.S. Environmental Protection Agency Level III) using combined S1 and S2 seasonal composite imagery with spectral indices and textural features; 2) assess the agreement of this Sentinel-based map with existing 30-m products; and 3) analyze and compare aspen patch size and total area across the Southern Rockies and within a targeted Colorado geography where aspen is currently a management focus. Finally, we provide access to results and input data through a public interactive GEE application and make all code and methodology available in public repositories. We anticipate that a higher resolution and more accurate map of aspen cover will facilitate targeted and effective management of aspen forests to increase the valued ecosystem services provided by this forest type.

2. Materials and Methods

2.1. Study Area

The Southern Rockies includes portions of southern Wyoming, central and western Colorado, and northern New Mexico. Extending nearly 130k km²;, the elevation of the region ranges from approximately 1,128-4,268 m above sea level and is characterized by two major mountain belts and intermontane valleys and parks. There are a diverse array of ecosystem types including lower montane/foothills, wet meadows, upper montane, subalpine, and alpine [39]. Forested area accounts for an estimated 72% (93.6k km²) of the landscape (Landfire EVT, c. 2016), with dominant tree species including Engelmann spruce and subalpine fir (18.1%), ponderosa pine and Douglas-fir (17.2%), quaking aspen (10.1%), lodgepole pine (9.4%), and pinon and juniper woodland (7.8%) [13]. The Southern Rockies is characterized by a large proportion of federally managed lands (87.5%) including ten National Forests, two National Parks, and 57 federally listed wilderness areas [40]. Quaking aspen is the dominant deciduous forest type and persists on the landscape in both pure (stable) and mixed (seral) stands and their geographic distribution is highly linked with disturbance history [4,41].

We evaluate our Sentinel-based aspen map across the Southern Rockies blocks and for a targeted geography in Colorado, the White River National Forest (NF), where quaking aspen is actively managed for wildlife habitat and wildfire risk reduction. The White River NF covers over 10,000 km²; in central and western Colorado and encompasses many of the ecosystem types of the Southern Rockies. The region is also recreationally important with 11 ski resorts, ten mountain peaks over 4,267 m, and eight wilderness areas encompassing nearly a third of the total area (https://www.fs.usda.gov/whiteriver; accessed March 31, 2023). The USFS manages 1,966 km²; of both stable and seral stands of quaking aspen. Current management activities aim to promote aspen regeneration in existing aspen stands, improve wildlife habitat by expanding aspen cover, and increase forest disturbance severity and recovery by promoting aspen in conifer stands. Managers are using mechanical harvests (coppice silviculture) and prescribed fire to achieve these outcomes.

Large area image classification presents unique challenges, particularly in the generation of representative and accurate training and validation data, collection of satellite imagery, and assessment of model performance [18]. To overcome these challenges, we divided the Southern Rockies into 129 equal area (50 km²) spatial blocks (Figure 1). These blocks provide an analytical unit to generate training and validation data (Section 2.2), imagery composites (Section 2.3.1 and 2.3.2) and perform spatial block cross-validation during model training (Section 2.4).

2.2. Reference Data

To train and validate the classification models, we generated presence and background learning data for the study region. In binary classification algorithms, presence (or the positive class) represents the target land cover type or species, and background (or the negative class) captures a representative sample of all other cover types. This approach is particularly useful for applications where the desired result is a single land cover type or species and is commonly used in applications of habitat suitability modeling [42].

2.2.1. Quaking Aspen Reference Data

Availability of field-based measurements of quaking aspen forest cover in the region is limited, particularly data aligned with Sentinel pixels. To overcome this challenge, we created a database of aspen presence using photo interpretation of high-resolution (1 m) aerial imagery from the United States Department of Agriculture (USDA) National Agricultural Imagery Program (NAIP, c. 2017) [43]. We created a binary mask of aspen forest cover using the USFS TreeMap (c. 2016) basal area estimates derived from the USFS FIA. We identified pixels where live basal area was greater than 10 stems/ha, a potential minimum unit for established aspen stands in the Southern Rockies [2]. For each spatial block, we placed 10 random points inside the mask and applied a 5 km buffer which served as a spatial support to aid interpretation. An interpreter was trained to identify aspen from NAIP imagery and attempted to label a minimum of 10 pixels of aspen within the spatial support unit. If there were no obvious stands of aspen, a new random point was generated in the block. A random sample of these presence points was taken with a minimum distance of 100 m to limit effects of pseudoreplication. The resulting database includes 12,609 points across 95 spatial blocks in the Southern Rockies.

2.2.2. Background Reference Data

Background reference data are representative of other major land cover types across the Southern Rockies and were derived from the Landfire EVT Sub-class property (c. 2016). To account for regional variation in cover types, proportional stratified samples were generated for each spatial block. The total number of samples in each block was set to ten times the number of presence points and divided amongst classes proportional to their area. A minimum sample size of 10 was set for each class in the block to ensure representation of minority classes, resulting in 73,031 samples across Landfire EVT sub-classes (Table 1).

2.3. Satellite Imagery

We collected freely available Sentinel-1 and Sentinel-2 imagery for one year (2019) in the GEE data catalog [19]. To develop seasonal imagery composites which are consistent across a large and diverse geographic region, we use phenological characteristics of aspen forests to identify temporal windows of key stages of development. Collection and preparation of median seasonal radar composites is described in Section 2.2.1 and optical imagery collection and compositing is described in Section 2.2.2 and 2.2.4, respectively. In addition, we calculated spectral vegetation indices and radar textural features which are described in Section 2.2.3.

2.3.1. Sentinel-1 Data

We collected winter (December-February) and summer (June-August) imagery from the Sentinel-1 (S1) C-Band Ground Range Detected (GRD), which collects data in dual polarization mode (VV and VH). These seasonal windows capture two distinct development stages of the forest canopy: peak (summer) and dormancy (winter). The S1 GRD collection was processed to analysis-ready data involving border noise correction and speckle filtering [44]. We combined both orbital passes (ascending and descending) for each polarization mode and calculated the median backscatter coefficient for the seasonal windows. In total, 556 individual S1 images were used to generate the median composites (324 for summer and 232 for winter).

2.3.2. Sentinel-2 Data

We collected all Sentinel-2 (S2) Multispectral Instrument (MSI) Level-2A surface reflectance scenes for the entire year across the Southern Rockies (5,137 tiles). Each surface reflectance tile includes four spectral bands at 10-m spatial resolution, six bands at 20-m, and three bands at 60-m. The 60-m bands (coastal aerosol, water vapor, SWIR-cirrus), which are uninformative for vegetation analysis, were excluded from the collection and remaining bands were resampled to a common 10-m spatial resolution (Table 1).

2.3.3. Additional Spectral and Textural Feature

In addition to the S1 backscatter and S2 spectral bands, we calculated a suite of additional textural features and vegetation indices (Table 2). Given that radar backscatter imagery is sensitive to canopy surface roughness, textural features of S1 data can provide valuable information for discriminating forest types. Specifically, the Gray Level Co-Occurrence Matrix (GLCM) calculates textural features from image data by assessing the patterns of pixel intensities and spatial arrangement [38]. The derivatives of GLCM have been shown to improve land cover classifications from radar data [32,38,45]. We derived GLCM entropy, contrast, variance, and correlation using a 7x7 kernel for each of the winter and summer VH and VV polarization modes resulting in 16 additional features for classification.

Six vegetation indices targeting forest canopy productivity, moisture, and chlorophyl content were derived from the S2 bands. These included the Chlorophyll Index Red-edge (CIRE), Modified Chlorophyll Absorption in Reflectance Index (MCARI), Inverted Red-edge Chlorophyll Index (IRECI), Specific Leaf Area Vegetation Index (SLAVI), Modified Normalised Difference Water Index (MNDWI), and the Red-edge Normalised Difference Vegetation Index (NDVI705).

2.3.4. Seasonal Sentinel-2 Composites

To generate median spectral image composites which are representative of key phenological stages of aspen, we extracted 1) the spectral response of aspen reference data (Section 2.3.1) for each image in the annual collection, and 2) the phenological characteristics within the USFS TreeMap mask (Section 2.3.1) from the Visible Infrared Imaging Radiometer Suite (VIIRS) Land Cover Dynamics data product (VNP22Q2), which provides global land surface phenology (GLSP) metrics at yearly interval [51]. For the S2 time-series, we excluded any observations which were occluded by clouds or shadows using the Cloud Score Plus [52] with a threshold value for occluded pixels of 0.60. We then calculated the median bi-weekly surface reflectance for each reference point in the time-series for the S2 bands and vegetation indices. The median and range of phenology metrics (in day-of-year) were calculated at the spatial block unit for the VIIRS time-series (2013-2022). The timing of phenological transition varies across spatial block units and years and is heavily influenced by elevation (Figure S1). A more detailed description of the phenology assessment across the Southern Rockies is provided in Appendix S1. Using these data, we defined two ecologically relevant temporal periods for image composites: summer (median onset of greenness increase to the 10^th percentile onset of greenness decrease, May 25-August 13, 2019), and autumn (minimum mid-senescence phase to the 90^th percentile onset of greenness minimum, September 2-November 14).

For each seasonal window, we filtered out images with greater than 80% cloud cover and 10% snow cover. We joined each image in the collection to its corresponding Cloud Score Plus image. The Cloud Score Plus provides a single pixel occlusion score based, enabling flexible and accurate masking of cloud and cloud shadows [52]. A binary cloud mask for each image was created with a cloud score threshold of 0.60 and used to remove occluded pixels. Seasonal image composites were created by calculating the per-pixel median surface reflectance value. A total of 3,956 S2 tiles were used with per-spatial block averages of 146 (+/- 71) tiles for summer and 73 (+/- 30) tiles for autumn.

2.3.5. Topographic Data

Inclusion of topographic information is shown to improve the accuracy of forest type classifications, especially in regions with complex elevation patterns [33]. We include a digital elevation model (DEM) from the United States Geological Service (USGS) 3-dimensional Elevation Program (3DEP) 1/3 arc-second (10 m) elevation product. In addition, we calculated aspect and slope from the DEM, both of which may influence the distribution of aspen forest cover across the landscape [41].

2.4. Image Classification

We developed a binary random forest (RF) classification of aspen cover. For image classification tasks, RF is a robust method which is less sensitive to overfitting while also minimizing computational expense with large input data [42,53]. In order to account for spatial autocorrelation in the training and testing data, we implemented 10-fold spatial block cross-validation during model development and calculated the average model performance metrics and confidence intervals [54]. For each fold, a 70:30 ratio was used to split valid spatial blocks (i.e., blocks which contain reference data) into training and testing sets respectively. Within each subset, the presence and background learning data were extracted, resulting in different subsets of reference data for each fold. Models were trained with 1001 trees and the number of random variables used at each split was set to the square root of the number of input variables. We assessed model performance across different combinations of classification scenarios and selected the best performing combination of input features to use for the final classification.

2.4.1. Model Selection and Accuracy Assessment

We assessed model performance using a set of common metrics based on the confusion matrix for each k-fold model and classification scenario. From the confusion matrix, we calculated the precision, recall and F1-score. To account for class imbalance in our training data, the F1-score was adopted to assess model performance because it typically performs better than overall accuracy in imbalanced classifications [55]. We selected the best performing classification scenario based on the average F1-score across folds. For the best performing scenario, we identified and removed multicollinear bands and performed feature selection to identify the most parsimonious model for classification in GEE. While the predictive ability of RF is likely unaffected by multicollinearity, the assessment of feature importance is sensitive to its presence [53]. We tested for multicollinearity using the ‘multi.collinear’ function in the rfUtilities [56] package in the R Statistical Programming Language [57] with a permutated (N=1001, p<=0.05) leave-one-out method. We removed features from best performing classification scenario based on the frequency that they were identified in the permutation tests (frequency > 0). We then used the ‘rf.modelSel’ function to identify the most parsimonious set of features based on out-of-bag error rate to use in the final classification model in GEE, thereby improving computation time and minimizing noise [58].

We assessed the accuracy of the final model based on the binary classification of the probability surface generated by the RF model. Typically, predicted classes are defined based on a single probability threshold between 0-1 (e.g., probability > 0.5 belongs to class 1). We tested the sensitivity of model performance using a moving threshold where the confusion matrix is calculated for 100 values between 0-1 for each of the 10 folds. At each threshold in the sequence, we classified the testing data into aspen (1) or non-aspen (0) and calculated the precision, recall, and F1-score. We identified the optimal threshold for classification based on the average maximum F1-score across folds.

2.4.2. Feature Importance

We retrieved feature importance from the final model, calculated in GEE as the sum of decrease in Gini impurity index over all trees in the forest [59]. We stored the feature importance values for each fold and present the median values with confidence intervals.

2.5. Agreement with Existing Products

We assessed the agreement between our Sentinel-based aspen map and three existing 30-m products: the Landfire EVT (c. 2016), USFS ITSP (c. 2014), and USFS TreeMap (c. 2016). We compare both the accuracy of our testing data (i.e., aspen presence data held out for model testing) as well as the pixel-based agreement across the Southern Rockies and the White River NF for each product. Each of these reference products were reclassified into aspen and non-aspen classes and resampled to a common 10-m spatial resolution to match the Sentinel-based map. For the USFS ITSP and TreeMap, we used the basal area metric (> 10 stems/ha) to reclassify the grids into aspen and non-aspen. A confusion matrix was created using the testing data across the region and for each spatial block across the Southern Rockies and for the White River NF for the pixel-based assessment. From the confusion matrices, we calculated the precision, recall, and F1-score and present the distribution of these metrics across the study regions.

2.6. Case Study: Landscape Patch Dynamics

To highlight the utility of the Sentinel-based map for management activities, we calculated patch- and class-level landscape metrics across the Southern Rockies and within the White River NF. We compare landscape metrics between the Sentinel-based map and the three reference datasets. These metrics included landscape-level (total area, patch density, and proportion of the landscape) and patch-level (patch size and perimeter to area ratio) statistics. All metrics were calculated using the pylandstats Python package [60].

3. Results

3.1. Annual Spectral Response of Quaking Aspen Forests

The annual surface reflectance of quaking aspen forests highlights the patterns of phenology, canopy development and environmental conditions in the study region (Figure 2A). During winter (December-March), surface reflectance is highly variable as the canopy is in dormancy and stands are often covered in snow, especially at higher elevations. As the canopy develops, a peak in the near-infrared (NIR) and red-edge wavelength occurs around the first week or two of July and corresponds with a mature forest canopy. As the canopy transitions to early autumn, there is a drop in NIR reflectance and increase in both the visible and shortwave-infrared (SWIR) which corresponds to a decrease in green vegetation. A similar pattern is observed for the S2-MSI spectral indices (Figure 2B). The largest difference in reflectance between summer and autumn windows occurs in the NIR and red-edge regions of the spectrum, although there are significant spectral changes across all S2 bands (Figure 2C) and spectral indices (Figure 2D).

3.2. Model Selection and Accuracy Assessment

We selected the final classification model based on the results from the scenario testing. The best performing classification scenario included all features from both S1 and S2 and topographic data with an average F1-score of 0.91 (+/- 0.01, Figure 3). The worst classification scenarios were those in which only S1 data was used (F1-score 0.46-0.62). There was a significant increase in classification accuracy for single-season S2 data (i.e., summer or autumn alone). Combining both the seasonal S2 composites along with spectral indices increased the average F1-score by 0.3 compared with just the summer S2 and indices, whereas the addition of S1 only improved the F1-score by an additional 0.1.

Thus, our final scenario included all S1 and S2 features for each seasonal window along with topographic features. For this best performing classification scenario, we removed 21 collinear bands and identified a set of 17 which minimized the out-of-bag error rate for the most parsimonious classification. We achieved an average F1-score of 0.931 (+/- 0.008) for the final model.

Feature Importance

We extracted the feature importance from the final model based on the sum of decrease in Gini Impurity Index for each model fold. While there were 17 features in the final model, we present only the top 10 most important (Figure 4). The Modified Normalised Difference Water Index (MNDWI) from the autumn S2 composite was the top predictor across all folds. Elevation was the second most important feature followed by the Chlorophyll Index Red-edge (CIRE) and Band 3 (green) from the summer S2 composite and the VV polarization band from the summer S1 composite. Just outside of the top 5 features is the Specific Leaf Area Vegetation Index (SLAVI) from the summer S2 composite. The remaining features in the top 10 are the red bands for both summer and autumn S2, Band 5 (red-edge) for the summer S2, and autumn CIRE. Of the 17 bands in the final model, the radar textural features were the least important in classifying aspen forests.

3.3. Quaking Aspen Forest Map

The resulting high-resolution (10 m) Sentinel-based map of quaking aspen forest cover represents the average probability across the 10 folds for the final model (Figure 5A, 5C). The optimum probability threshold for classification of 0.42 was used to classify the final map into aspen and non-aspen pixels based on the maximum F1-score achieved across the range of thresholds tested for each fold (Figure S2B). A more detailed description of the threshold testing and model performance can be found in Appendix S2. Across the Southern Rockies, we estimate 9,482 km² (7.4% of total land area) as quaking aspen forests. Within the White River NF, we mapped 1,658 km² (16.3% of total land area).

3.4. Agreement with Reference Datasets

We tested the agreement of the Sentinel-based map with three existing products using both the testing data from model training and a pixel-based approach across the entire study region and within the White River NF. Based on testing data (i.e., aspen presence points held out of model training), we see high overall agreement between all products with average F1-score ranging from 0.829-0.865 compared to the 0.931 for the Sentinel-based map (Table 4). However, a notable pattern emerges regarding precision and recall. The Sentinel-based map demonstrates higher precision (0.9516) relative to recall (0.9116). This suggests that while the map is correctly predicting aspen presence 95% of the time, it may be more restrictive, missing about 9% of true aspen presence. Conversely, all three reference products exhibit lower precision than recall. For example, the USFS TreeMap correctly predicts aspen presence just 80% of the time but identifies 93% of all true aspen presence. This indicates that while existing products are relatively accurate, they may be more likely to include non-aspen areas in their predictions based on the testing data.

In the pixel-based agreement assessment between the three existing products and the Sentinel-based map, we calculated the precision, recall, and F1-score across spatial blocks in the Southern Rockies and for the White River NF. There was significant spatial variability in all metrics across spatial blocks with some regions performing significantly better than others (Figure 6). Interestingly, we observe an opposite relationship between precision and recall when comparing the maps directly suggesting that the Sentinel-based map remains restrictive in classifying aspen presence compared to the existing products. Given the higher accuracy of the Sentinel-based map based on testing data (Table 4), this supports the finding that existing products likely over-predict areas of aspen presence. Additionally, the White River NF exhibits significantly higher precision, recall, and F1-score for all existing products than the median across spatial blocks (Figure 6). This region is characterized by extensive aspen forests, suggesting that the agreement may be higher in areas with more consistent aspen cover. Still, the higher precision than recall in the pixel-based assessment indicates that there is likely overestimation occurring for the three reference products. Along with assessment of the test data, the Sentinel-based map shows more accurate identification of actual aspen forests on the landscape.

Landscape and Patch Dynamics

We assessed the landscape patch dynamics for the Sentinel-based map and the three existing products by spatial block across the Southern Rockies and for the White River NF. The Sentinel-based map estimates 43-118% less total aspen area across the Southern Rockies and 21-64% less in the White River NF compared to existing products (Figure 7A). The average patch size varies greatly across the Southern Rockies, highlighting the regional difference in stand characteristics. For the Sentinel-based map, the average patch size is 0.53 ha (+- 23) in the Southern Rockies and 0.74 ha (+- 22) in the White River NF, significantly smaller than existing products (Table 5, Figure7B). Some regions are dominated by large stands of pure aspen (> 22 ha), whereas others are defined by many small, more dispersed stands (Figure 7B and 7C). Notably, the Sentinel-based map identifies 28-93% more individual patches across the Southern Rockies and White River NF. Given the higher accuracy of the Sentinel-based map and the trend in over-prediction of aspen areas from existing products, these results suggest that the higher-resolution map is better able to characterize the actual landscape configuration of aspen forests.

4. Discussion

This overall effort provides a high-resolution quaking aspen forest cover map from 2019 for the Southern Rockies. We map 9,482 km² (7.4% of total area) of aspen forest cover within the SRME and 1,658 km2 (16.3% of total area) within the White River NF region. Leveraging the GEE cloud-compute platform, we were able to map quaking aspen forests across a large geographic region. We used multi-temporal composite satellite imagery, presence and background learning data and a spatial block cross-validation approach to train a random forest classification model which achieved high accuracy (average F1-score of 0.93). We demonstrate that combining seasonal spectral composites (summer and autumn) significantly improves classification accuracy of quaking aspen forests, while the addition of S1 imagery only marginally improves accuracy. In particular, S2 spectral indices targeting canopy moisture and chlorophyl content were important in the classification of aspen forests. Furthermore, we assessed the agreement of this new high-resolution map with existing products, identifying spatial consistencies across the region while achieving a higher level of accuracy when compared to training data. The higher-resolution map facilitates improved understanding of landscape dynamics such as the average patch size, patch density, and overall area. These metrics provide key insights into the management of individual stands across diverse areas.

Similar to previous studies using Sentinel data to classify forest species, we achieve high classification accuracy (average F1-score of 0.93), especially when leveraging multiple seasonal image composites. Previous studies have shown the effectiveness of multi-temporal S2 imagery for vegetation mapping in heterogeneous environments [26,27,28,29,33]. However, these studies often use single S2 tiles to retrieve cloud-free images across time for a given location. Large-area image classification is impacted by persistent cloud-cover and variations in the timing and availability of cloud-free images on a single day [18]. We show that multi-temporal gap-filling and compositing methods to generate cloud-free images is an effective method for overcoming these challenges. Additionally, existing data on land surface phenology is helpful in identifying ecologically relevant temporal windows to create seasonal imagery composites. This method can be applied to other species with unique phenological characteristics.

In western US mixed conifer forests, quaking aspen canopies contain higher levels of canopy moisture and chlorophyl content than other species. We showed that spectral indices from S2 seasonal composites which leverage the red-edge and SWIR bands to estimate canopy moisture and chlorophyl content are important features for classification. For example, two of the most important features in the final model were the Modified Normalized Water Index (MNDWI) and the Chlorophyl Red-edge Index (CIRE) for the summer image composite. Leveraging unique canopy characteristics to created targeted spectral indices improves discrimination of forest species, particularly deciduous species.

In comparison with three primary existing products, we achieve higher precision in the prediction of aspen presence. Existing maps are developed at a coarser spatial resolution (30-m) and tend to over-predict areas of aspen cover, as highlighted by the tradeoffs between precision and recall when comparing the products using known locations of aspen presence. For the Sentinel-based map, we achieve significantly higher precision than recall which contrasts with existing maps. This suggests that the new map is highly accurate in correctly identifying aspen pixels but may be more restrictive. Conversely, existing products typically do identify aspen forests accurately, but tend to include many non-aspen areas in their classifications. However, these trends vary across the Southern Rockies as shown in the distribution of the pixel-based agreement metrics across spatial blocks. Areas with higher consistency between the Sentinel-based map and existing products based on the pixel-level confusion matrices seem to occur in regions with more consistent or homogenous aspen forest cover. Landscape heterogeneity, aspen stand dynamics (e.g., pure or mixed/seral), or various forest densities may influence the success of predictions using remotely sensed data. This underscores the importance of spatial assessments, as the results provide an opportunity to investigate areas with poor precision or recall for model tuning and classification improvements in the future. In these areas, development of more precise training and validation data may improve the overall ability of the model to distinguish aspen forests from other land cover types.

Furthermore, the landscape dynamics of aspen forests vary widely across the Southern Rockies, with some areas dominated by pure or stable stands and others characterized by mid-late successional seral stands. Using the Sentinel-based map, we demonstrate the effectiveness for higher-resolution data to provide detailed insights into the landscape configuration based on metrics like average patch size, patch density, and perimeter to area ratio. We show that these metrics vary across both the Southern Rockies and White River NF. Importantly, the Sentinel-based map can also detect smaller stands (<0.09 ha) and greater perimeter to area ratio, suggesting it is able to better characterize the landscape configuration of aspen forests. These insights are critical in developing management plans targeting different stages of successional development in aspen forests.

Together, this effort provides a crucial tool for land managers in determining where to prioritize active management of aspen forests. Quaking aspen is a shade-intolerant species capable of vegetation regeneration from root suckers, as such it readily occupies canopy openings [14,61]. These characteristics have been used to regenerate and even expand aspen patches using coppice silviculture techniques. The improved capability to map patch dynamics at fine spatial resolution and high accuracy is useful for forest managers who are planning silvicultural management activities which maximize the potential for expansion of aspen patches into new areas. These smaller patches can offer prolific root suckering when managed or following disturbances (Long and Mock 2012; Shepperd et al. 2015). Combined with other datasets including habitat suitability, wildfire or other disturbance risk, and others, these new maps have the potential to assist forest management at both temporal and spatial scales relevant to the ecology of the species.

While we achieve high classification accuracy in this new map, there are improvements that could be made. Addition of field-based measurements of aspen forests across the region would improve accuracy assessment and model refinement. Collecting data from existing field studies, citizen science databases (e.g., iNaturalist), or species occurrence databases (e.g., the Global Biodiversity Information Facility) may provide both additional accuracy assessment opportunities and model improvements while limiting the need for time-intensive photo interpretation efforts. Future efforts may include more state-of-the-art classification techniques such as neural networks and object-oriented classification, which would take advantage of the spatial characteristics of aspen stands. Moreover, the use of upcoming hyperspectral imaging campaigns, very high spatial resolution imagery, or other active remote sensing such as LiDAR will undoubtedly offer continued improvements. The methods presented in this research are fully open and reproducible, enabling others to contribute as the availability of remotely sensed products and cloud-compute resources continues to grow.

5. Conclusion

As human- and climate-related impacts to forested areas continue to grow, the development of new high-resolution and accurate maps of forest species is critical to effective management activities which improve resilience, reduce hazard for communities, and maintain important ecosystem service benefits. We demonstrate an effective method for large-area forest species classification and provide an improved map of quaking aspen forests for the Southern Rockies. Given the uncertainty around aspen growth dynamics in a changing climate, this product is an important source of information for forest managers and will benefit management of quaking aspen forests.