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

Quaking aspen (Populus tremuloides Michx.) is an important deciduous species in western U.S. forests. Existing maps of aspen distribution are based on Landsat imagery and often miss small stands (<0.09 ha or 30 m²), which rapidly regrow when managed or following disturbance. In this study, we present methods for deriving a new regional map of aspen forests using one year of Sentinel-1 (S1) and Sentinel-2 (S2) imagery in Google Earth Engine. We assessed the annual phenology of aspen in the Southern Rockies and leveraged the frequent temporal resolution of S1 and S2 to create ecologically relevant seasonal composites. Additionally, we derived spectral indices and radar textural features targeting the canopy structure, moisture, and chlorophyll content. Using spatial block cross-validation and Random Forests, we assessed the accuracy of different scenarios and selected the best performing set of features for classification. Comparisons were then made with existing landcover products across the study region. The resulting map improves on existing landcover products in both accuracy (average F1-score = 93.1%) and detection of smaller forest patches. These methods enable accurate mapping at spatial and temporal scales relevant to forest management for one of the most widely distributed tree species in North America.

Keywords:

Subject: Environmental and Earth Sciences - Remote Sensing

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.

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org. Figure A1: Quaking aspen phenology in the Southern Rockies; Figure B1: Model performance metrics from the final classification.

Author Contributions

Conceptualization, M.C, T.C., S.H., and J.B.; methodology, M.C. and T.C.; software, M.C.; formal analysis, M.C.; investigation, M.C. and T.C.; data curation, M.C.; writing—original draft preparation, M.C.; writing—review and editing, M.C., T.C., S.H., A.P, and J.B.; visualization, M.C.; supervision, S.H. and J.B.; project administration, S.H. and J.B; funding acquisition, S.H., J.B., and M.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Joint Fire Science Program (JFSP, Project ID 21-2-02-29).

Data Availability Statement

Data associated with this paper are available as public Google Earth Engine assets. An interactive web viewing application provides access to the results and input satellite data (https://maco.users.earthengine.app/view/aspen-viewer). All analyses and code used in this paper are available in a public GitHub repository (https://github.com/maxwellCcook/aspen-fire) and a public Google Earth Engine repository (https://code.earthengine.google.com/?accept_repo=users/maco/aspen-fire).

Acknowledgments

Support from the Earth Lab, Cooperative Institute for Research in Environmental Science (CIRES), Department of Geography University of Colorado Boulder, and the Department of Forest and Rangeland Stewardship Colorado State University. The authors would like to acknowledge Cibele Do Amaral, Catherine Schloegel, and Kyle Rodman for their intellectual contributions and for discussions related to the manuscript and undergraduate research assistant Adam Thomas for his dedicated work towards developing the aspen presence data.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A

Phenology of Quaking Aspen across the Southern Rockies

The Southern Rockies study region is an ecologically diverse landscape with complex topography. In complex systems, exploring the phenological patterns of vegetation development is useful when generating seasonal satellite image composites [25]. For quaking aspen forests in particular, the timing of canopy growth differs by elevation, slope, and climate [62]. These differences across a large and diverse region can influence the generation of satellite image composites which are representative of specific phenological stages for this deciduous forest species.

To examine the seasonal variation of canopy development within aspen forests across the Southern Rockies, we summarized phenological metrics from the Visible Infrared Imaging Radiometer Suite (VIIRS) Land Cover Dynamics data product (VNP22Q2), which provides global land surface phenology metrics at yearly intervals since 2013 [51]. These include day-of-year (DOY) for the onset of greenness increase, mid-greenup phase, onset of greenness maximum, onset of greenness decrease, mid-senescence phase, and the onset of greenness minimum. For each spatial block, we calculated the median day-of-year for each metric within aspen forest area based on the USFS TreeMap product (see Section 2.3.1 of the main text). Analysis was carried out in the Google Earth Engine (GEE) Code Editor [19]. We examined the variation in the average DOY across spatial blocks (Figure A1) for each metric and calculated the average and standard deviation across blocks and years for the Southern Rockies. The variation in quaking aspen phenology is well-explained by elevation as demonstrated by the inset plots in Figure A1. To capture this variation when selecting temporal windows for creating seasonal image composites, we subtracted the standard deviation (in days) from the metric average (day-of-year) and used these values as start and end dates. Specifically, we defined two temporal windows: summer, or the mid-greenup phase to the onset of greenness decrease; and autumn, or the mid-senescence phase to the onset of greenness minimum.

Figure A1. Quaking aspen phenology derived from the VNP22Q1 in the Southern Rockies and relationship with elevation. Day-of-year estimates represent the median across all years (2013-2022). (A) Mid-greenup phase (late May to late June) and (B) the onset of greenness decrease (early August to early September) defines the summer temporal window, and both show a strong positive linear relationship with elevation. (C) Mid-senescence phase (mid-September to early October) and (D) the onset of greenness minimum (late October to mid-November) defines the autumn temporal window. Senescence phase indicates a strong relationship with elevation while the onset of greenness minimum is less influenced by elevation.

Appendix B

Accuracy Assessment and Optimal Threshold for Classification

To assess model performance and to calculate the optimum threshold for classification based on the probability grids from the random forest classifier, we calculated a confusion matrix for 100 different classification threshold values between 0-1 based on the training data in each of the 10 model iterations. From this, we calculated an AUC-ROC curve which represents the true positive rate against the false positive rate (Figure B1A). To calculate the optimum threshold, we found the maximum F1 score of the average across all 10 model iterations (Figure B1B). We used this optimum threshold on the average probability of aspen forest grid to calculate the binary output (aspen and no aspen).

Figure A2. Model performance metrics from the final classification. (A) AUC-ROC or the true-positive rate against the false-positive rate with colors representing the 10 model iterations; (B) The F1 score across threshold values. The black line is the average across model iterations and the black point represents the location of the maximum F1 from the average across folds. At this point is our optimum threshold for classification (0.42).

References

Rogers, P.C.; Pinno, B.D.; Šebesta, J.; Albrectsen, B.R.; Li, G.; Ivanova, N.; Kusbach, A.; Kuuluvainen, T.; Landhäusser, S.M.; Liu, H.; et al. A Global View of Aspen: Conservation Science for Widespread Keystone Systems. Global Ecology and Conservation 2020, 21, e00828. [Google Scholar] [CrossRef]
DeByle, N.V.; Winokur, R.P. Aspen: Ecology and Management in the Western United States. USDA Forest Service General Technical Report RM-119. Rocky Mountain Forest and Range Experiment Station, Fort Collins, Colo. 283 p. 1985, 119. [Google Scholar] [CrossRef]
Bartos, D.L. Landscape Dynamics of Aspen and Conifer Forests. In: Shepperd, Wayne D.; Binkley, Dan; Bartos, Dale L.; Stohlgren, Thomas J.; Eskew, Lane G., comps. Sustaining aspen in western landscapes: Symposium proceedings; 13-15 June 2000; Grand Junction, CO. Proceedings RMRS-P-18. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. p. 5-14. 2001, 18, 5–14. [Google Scholar]
Rogers, P.C.; Landhäusser, S.M.; Pinno, B.D.; Ryel, R.J. A Functional Framework for Improved Management of Western North American Aspen (Populus Tremuloides Michx.). Forest Science; Bethesda 2014, 60, 345–359. [Google Scholar] [CrossRef]
Landhäusser, S.M.; Deshaies, D.; Lieffers, V.J. Disturbance Facilitates Rapid Range Expansion of Aspen into Higher Elevations of the Rocky Mountains under a Warming Climate. Journal of Biogeography 2010, 37, 68–76. [Google Scholar] [CrossRef]
Gill, N.S.; Sangermano, F.; Buma, B.; Kulakowski, D. Populus Tremuloides Seedling Establishment: An Underexplored Vector for Forest Type Conversion after Multiple Disturbances. Forest Ecology and Management 2017, 404, 156–164. [Google Scholar] [CrossRef]
Andrus, R.A.; Hart, S.J.; Tutland, N.; Veblen, T.T. Future Dominance by Quaking Aspen Expected Following Short-Interval, Compounded Disturbance Interaction. Ecosphere 2021, 12, e03345. [Google Scholar] [CrossRef]
Nigro, K.M.; Rocca, M.E.; Battaglia, M.A.; Coop, J.D.; Redmond, M.D. Wildfire Catalyzes Upward Range Expansion of Trembling Aspen in Southern Rocky Mountain Beetle-Killed Forests. Journal of Biogeography 2022, 49, 201–214. [Google Scholar] [CrossRef]
Worrall, J.J.; Rehfeldt, G.E.; Hamann, A.; Hogg, E.H.; Marchetti, S.B.; Michaelian, M.; Gray, L.K. Recent Declines of Populus Tremuloides in North America Linked to Climate. Forest Ecology and Management 2013, 299, 35–51. [Google Scholar] [CrossRef]
Rosenblum, A. Altered Fire Regimes and the Persistence of Quaking Aspen in the Rocky Mountains: A Literature Review. Open Journal of Forestry 2015, 5, 563–567. [Google Scholar] [CrossRef]
Krasnow, K.D.; Stephens, S.L. Evolving Paradigms of Aspen Ecology and Management: Impacts of Stand Condition and Fire Severity on Vegetation Dynamics. Ecosphere; Washington 2015, 6, 1–16. [Google Scholar] [CrossRef]
Shinneman, D.J.; McIlroy, S.K. Climate and Disturbance Influence Self-Sustaining Stand Dynamics of Aspen (Populus Tremuloides) near Its Range Margin. Ecological Applications 2019, 29, e01948. [Google Scholar] [CrossRef] [PubMed]
Picotte, J.J.; Dockter, D.; Long, J.; Tolk, B.; Davidson, A.; Peterson, B. LANDFIRE Remap Prototype Mapping Effort: Developing a New Framework for Mapping Vegetation Classification, Change, and Structure. Fire 2019, 2, 35. [Google Scholar] [CrossRef]
Shepperd, W.D.; Smith, F.W.; Pelz, K.A. Group Clearfell Harvest Can Promote Regeneration of Aspen Forests Affected by Sudden Aspen Decline in Western Colorado. Forest Science 2015, 61, 932–937. [Google Scholar] [CrossRef]
Landhäusser, S.M.; Pinno, B.D.; Mock, K.E. Tamm Review: Seedling-Based Ecology, Management, and Restoration in Aspen (Populus Tremuloides). Forest Ecology and Management 2019, 432, 231–245. [Google Scholar] [CrossRef]
Ellenwood, J.; Jr, F.; Romero, S. National Individual Tree Species Atlas; 2015;
Riley, K.L.; Grenfell, I.C.; Shaw, J.D.; Finney, M.A. TreeMap 2016 Dataset Generates CONUS-Wide Maps of Forest Characteristics Including Live Basal Area, Aboveground Carbon, and Number of Trees per Acre. JOURNAL OF FORESTRY 2022. [Google Scholar] [CrossRef]
Hermosilla, T.; Wulder, M.A.; White, J.C.; Coops, N.C. Land Cover Classification in an Era of Big and Open Data: Optimizing Localized Implementation and Training Data Selection to Improve Mapping Outcomes. Remote Sensing of Environment 2022, 268, 112780. [Google Scholar] [CrossRef]
Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-Scale Geospatial Analysis for Everyone. Remote Sensing of Environment 2017, 202, 18–27. [Google Scholar] [CrossRef]
Phiri, D.; Simwanda, M.; Salekin, S.; Nyirenda, V.R.; Murayama, Y.; Ranagalage, M. Sentinel-2 Data for Land Cover/Use Mapping: A Review. Remote Sensing 2020, 12, 2291. [Google Scholar] [CrossRef]
Bayle, A.; Carlson, B.Z.; Thierion, V.; Isenmann, M.; Choler, P. Improved Mapping of Mountain Shrublands Using the Sentinel-2 Red-Edge Band. Remote Sensing 2019, 11, 2807. [Google Scholar] [CrossRef]
Clevers, J.G.P.W.; Gitelson, A. Using the Red-Edge Bands on Sentinel-2 for Retrieving Canopy Chlorophyll and Nitrogen Content. European Space Agency, (Special Publication) ESA SP 2012, 707. [Google Scholar]
Lin, S.; Li, J.; Liu, Q.; Li, L.; Zhao, J.; Yu, W. Evaluating the Effectiveness of Using Vegetation Indices Based on Red-Edge Reflectance from Sentinel-2 to Estimate Gross Primary Productivity. Remote Sensing 2019, 11, 1303. [Google Scholar] [CrossRef]
Wong, C.Y.S.; D’Odorico, P.; Bhathena, Y.; Arain, M.A.; Ensminger, I. Carotenoid Based Vegetation Indices for Accurate Monitoring of the Phenology of Photosynthesis at the Leaf-Scale in Deciduous and Evergreen Trees. Remote Sensing of Environment 2019, 233, 111407. [Google Scholar] [CrossRef]
Kollert, A.; Bremer, M.; Löw, M.; Rutzinger, M. Exploring the Potential of Land Surface Phenology and Seasonal Cloud Free Composites of One Year of Sentinel-2 Imagery for Tree Species Mapping in a Mountainous Region. International Journal of Applied Earth Observation and Geoinformation 2021, 94, 102208. [Google Scholar] [CrossRef]
Immitzer, M.; Link to external site, this link will open in a new window; Neuwirth, M.; Böck, S.; Brenner, H.; Vuolo, F.; Atzberger, C.; Link to external site, this link will open in a new window Optimal Input Features for Tree Species Classification in Central Europe Based on Multi-Temporal Sentinel-2 Data. Remote Sensing 2019, 11, 2599. [CrossRef]
Li, H.; Shi, Q.; Wan, Y.; Shi, H.; Imin, B. Using Sentinel-2 Images to Map the Populus Euphratica Distribution Based on the Spectral Difference Acquired at the Key Phenological Stage. Forests 2021, 12, 147. [Google Scholar] [CrossRef]
Macintyre, P.; van Niekerk, A.; Mucina, L. Efficacy of Multi-Season Sentinel-2 Imagery for Compositional Vegetation Classification. International Journal of Applied Earth Observation and Geoinformation 2020, 85, 101980. [Google Scholar] [CrossRef]
Persson, M.; Lindberg, E.; Reese, H. Tree Species Classification with Multi-Temporal Sentinel-2 Data. Remote Sensing 2018, 10, 1794. [Google Scholar] [CrossRef]
Barber, B.C. Review Article. Theory of Digital Imaging from Orbital Synthetic-Aperture Radar. International Journal of Remote Sensing 1985, 6, 1009–1057. [Google Scholar] [CrossRef]
Udali, A.; Lingua, E.; Persson, H.J. Assessing Forest Type and Tree Species Classification Using Sentinel-1 C-Band SAR Data in Southern Sweden. Remote Sensing 2021, 13, 3237. [Google Scholar] [CrossRef]
Numbisi, F.N.; Numbisi, F.N.; Coillie, F.V.; Wulf, R.D. MULTI-DATE SENTINEL1 SAR IMAGE TEXTURES DISCRIMINATE PERENNIAL AGROFORESTS IN A TROPICAL FOREST-SAVANNAH TRANSITION LANDSCAPE. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2018; XLII–1, 339–346. [Google Scholar] [CrossRef]
Dobrinić, D.; Gašparović, M.; Medak, D. Sentinel-1 and 2 Time-Series for Vegetation Mapping Using Random Forest Classification: A Case Study of Northern Croatia. Remote Sensing 2021, 13, 2321. [Google Scholar] [CrossRef]
Lechner, M.; Dostálová, A.; Hollaus, M.; Atzberger, C.; Immitzer, M. Combination of Sentinel-1 and Sentinel-2 Data for Tree Species Classification in a Central European Biosphere Reserve. Remote Sensing 2022, 14, 2687. [Google Scholar] [CrossRef]
Liu, Y.; Gong, W.; Hu, X.; Gong, J. Forest Type Identification with Random Forest Using Sentinel-1A, Sentinel-2A, Multi-Temporal Landsat-8 and DEM Data. Remote Sensing 2018, 10, 946. [Google Scholar] [CrossRef]
Steinhausen, M.J.; Wagner, P.D.; Narasimhan, B.; Waske, B. Combining Sentinel-1 and Sentinel-2 Data for Improved Land Use and Land Cover Mapping of Monsoon Regions. International Journal of Applied Earth Observation and Geoinformation 2018, 73, 595–604. [Google Scholar] [CrossRef]
Dusseux, P.; Corpetti, T.; Hubert-Moy, L.; Corgne, S. Combined Use of Multi-Temporal Optical and Radar Satellite Images for Grassland Monitoring. Remote Sensing 2014, 6, 6163–6182. [Google Scholar] [CrossRef]
Haralick, R.M.; Shanmugam, K.; Dinstein, I. Textural Features for Image Classification. IEEE Transactions on Systems, Man, and Cybernetics, 1973; SMC-3, 610–621. [CrossRef]
R. K. Peete Forest and Meadows of the Rocky Mountains. In North American terrestrial vegetation; Cambridge University Press, New York, NY, 1988; pp. 63–102.
Johnson, R.M. Help Build the Protected Areas Database of the United States (PAD-US); Fact Sheet; Reston, VA, 2023;
Shepperd, W. A Classification of Quaking Aspen in the Central Rocky Mountains Based on Growth and Stand Characteristics. Western Journal of Applied Forestry 1990, 5, 69–75. [Google Scholar] [CrossRef]
Deng, X.; Li, W.; Liu, X.; Guo, Q.; Newsam, S. One-Class Remote Sensing Classification: One-Class vs. Binary Classifiers. International Journal of Remote Sensing 2018, 39, 1890–1910. [Google Scholar] [CrossRef]
Davis, D. 2017 NAIP Information Sheet.; US Department of Agriculture, 2011;
Mullissa, A.; Vollrath, A.; Odongo-Braun, C.; Slagter, B.; Balling, J.; Gou, Y.; Gorelick, N.; Reiche, J. Sentinel-1 SAR Backscatter Analysis Ready Data Preparation in Google Earth Engine. Remote Sensing 2021, 13, 1954. [Google Scholar] [CrossRef]
Farwell, L.S.; Gudex-Cross, D.; Anise, I.E.; Bosch, M.J.; Olah, A.M.; Radeloff, V.C.; Razenkova, E.; Rogova, N.; Silveira, E.M.O.; Smith, M.M.; et al. Satellite Image Texture Captures Vegetation Heterogeneity and Explains Patterns of Bird Richness. Remote Sensing of Environment 2021, 253, 112175. [Google Scholar] [CrossRef]
Gitelson, A.A.; Gritz †, Y.; Merzlyak, M.N. Relationships between Leaf Chlorophyll Content and Spectral Reflectance and Algorithms for Non-Destructive Chlorophyll Assessment in Higher Plant Leaves. Journal of Plant Physiology 2003, 160, 271–282. [Google Scholar] [CrossRef] [PubMed]
Frampton, W.J.; Dash, J.; Watmough, G.; Milton, E.J. Evaluating the Capabilities of Sentinel-2 for Quantitative Estimation of Biophysical Variables in Vegetation. ISPRS Journal of Photogrammetry and Remote Sensing 2013, 82, 83–92. [Google Scholar] [CrossRef]
Kobayashi, N.; Tani, H.; Wang, X.; Sonobe, R. Crop Classification Using Spectral Indices Derived from Sentinel-2A Imagery. Journal of Information and Telecommunication 2020, 4, 67–90. [Google Scholar] [CrossRef]
Evangelides, C.; Nobajas, A. Red-Edge Normalised Difference Vegetation Index (NDVI705) from Sentinel-2 Imagery to Assess Post-Fire Regeneration. Remote Sensing Applications: Society and Environment 2020, 17, 100283. [Google Scholar] [CrossRef]
Jiang, W.; Ni, Y.; Pang, Z.; He, G.; Fu, J.; Lu, J.; Yang, K.; Long, T.; Lei, T. A NEW INDEX FOR IDENTIFYING WATER BODY FROM SENTINEL-2 SATELLITE REMOTE SENSING IMAGERY. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2020; V-3–2020, 33–38. [Google Scholar] [CrossRef]
Zhang, X.; Liu, L.; Liu, Y.; Jayavelu, S.; Wang, J.; Moon, M.; Henebry, G.M.; Friedl, M.A.; Schaaf, C.B. Generation and Evaluation of the VIIRS Land Surface Phenology Product. Remote Sensing of Environment 2018, 216, 212–229. [Google Scholar] [CrossRef]
Pasquarella, V.J.; Brown, C.F.; Czerwinski, W.; Rucklidge, W.J. Comprehensive Quality Assessment of Optical Satellite Imagery Using Weakly Supervised Video Learning. In Proceedings of the 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW); June 2023; pp. 2125–2135. [Google Scholar]
Belgiu, M.; Drăguţ, L. Random Forest in Remote Sensing: A Review of Applications and Future Directions. ISPRS Journal of Photogrammetry and Remote Sensing 2016, 114, 24–31. [Google Scholar] [CrossRef]
Karasiak, N.; Dejoux, J.-F.; Monteil, C.; Sheeren, D. Spatial Dependence between Training and Test Sets: Another Pitfall of Classification Accuracy Assessment in Remote Sensing. Mach Learn 2021. [Google Scholar] [CrossRef]
Uhl, J.H.; Leyk, S. A Framework for Scale-Sensitive, Spatially Explicit Accuracy Assessment of Binary Built-up Surface Layers. Remote Sensing of Environment 2022, 279, 113117. [Google Scholar] [CrossRef]
Evans, J.S.; Murphy, M.A. rfUtilities; 2018;
R Core Team R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2022.
Evans, J.S.; Cushman, S.A. Gradient Modeling of Conifer Species Using Random Forests. Landscape Ecol 2009, 24, 673–683. [Google Scholar] [CrossRef]
Nembrini, S.; König, I.R.; Wright, M.N. The Revival of the Gini Importance? Bioinformatics 2018, 34, 3711–3718. [Google Scholar] [CrossRef] [PubMed]
Bosch, M. PyLandStats: An Open-Source Pythonic Library to Compute Landscape Metrics. PLOS ONE 2019, 14, e0225734. [Google Scholar] [CrossRef] [PubMed]
Long, J.N.; Mock, K. Changing Perspectives on Regeneration Ecology and Genetic Diversity in Western Quaking Aspen: Implications for Silviculture. Can. J. For. Res. 2012, 42, 2011–2021. [Google Scholar] [CrossRef]
Meier, G.A.; Brown, J.F.; Evelsizer, R.J.; Vogelmann, J.E. Phenology and Climate Relationships in Aspen (Populus Tremuloides Michx.) Forest and Woodland Communities of Southwestern Colorado. Ecological Indicators 2015, 48, 189–197. [Google Scholar] [CrossRef]

Figure 1. Map of the Southern Rockies study area showing the spatial block grid, case study landscape (White River NF), quaking aspen reference data and forest cover. The inset map highlights an area of aspen photo interpretation points with aerial imagery from the National Agricultural Imagery Program (NAIP, c. 2019). Forested area makes up an estimated 72% of land area in the Southern Rockies (National Land Cover Database Forest Canopy Cover c. 2019). The study area is divided into 129 equal area (50 km²) spatial blocks which serve as the analytical unit for image analysis and spatial cross-validation during model training.

Figure 2. Annual (A-B) and seasonal (C-D) spectral response of aspen presence data (c. 2019). Grey shading defines the temporal windows. (A) Median bi-weekly surface reflectance values for S2 bands; (B) Median weekly response of derived spectral indices; (C) Seasonal differences in S2 bands from image composites (summer and autumn); (D) Seasonal difference in spectral indices.

Figure 3. Results from the classification scenarios. We report the average and standard error F1-score across folds for each scenario. Topographic data (elevation, slope, and aspect) were included in all scenarios.

Figure 4. Feature importance for the top 10 features across folds for the final model measured as the sum of decrease in Gini impurity index. Color shade indicates the seasonality.

Figure 5. Map of quaking aspen forest cover (c. 2019) for the Southern Rockies with reference location (Hahns Peak, Steamboat Lake State Park, Routt County, CO). (A) Quaking aspen forest cover probability (0-100%); (B) True color imagery from NAIP for the reference location showing large stands of quaking aspen forests near Hahns Peak; (C) Quaking aspen forest cover probability (0-100%) for the reference location; (D) Quaking aspen classification based on the optimum threshold of the probability (0.42).

Figure 6. Distribution of pixel-based precision, recall, and F1-score of existing products compared to the Sentinel-based map across spatial blocks within the Southern Rockies. Diamonds indicate the values for the targeted geography (White River NF).

Figure 7. Comparison of patch metrics for the Southern Rockies and the White River NF for the Sentinel-based map and the three existing products. The White River NF is dominated by large, often pure stands of quaking aspen which is highlighted in the much larger average patch density and patch size when compared to the Southern Rockies.

Table 1. Distribution of background samples from the Landfire EVT Sub-class (c. 2016). The number of samples represents the total across all spatial blocks with at least 100 presence points. We removed any classes which included quaking aspen forests.

Landfire EVT Sub-class	Number of Samples
Evergreen closed tree canopy	23,957
Mixed evergreen-deciduous shrubland	13,421
Evergreen open tree canopy	12,991
Perennial graminoid grassland	7,637
Annual Graminoid/Forb	3,537
Evergreen shrubland	3,273
Sparsely vegetated	2,869
Mixed evergreen-deciduous open tree canopy	1,016
Developed	896
Non-vegetated	808
Perennial graminoid	716
Evergreen dwarf-shrubland	669
Evergreen sparse tree canopy	624
Deciduous open tree canopy	617
Total	73,031

Table 2. Sentinel-1 and Sentinel-2 input features. All bands were resampled to a common spatial resolution of 10-m.

Satellite	Abbrev.	Name	Center Wavelength	Seasonal Windows*
Sentinel-1	VV	Vertical-Vertical	5.5cm	Peak Greenness / Dormancy
Sentinel-1	VH	Vertical-Horizontal	5.5cm	Peak Greenness / Dormancy
Sentinel-2	B2	Blue	490 nm	Summer / Autumn
	B3	Green	560 nm
	B4	Red	665 nm
	B5	Red-edge 1	705 nm
	B6	Red-edge 2	740 nm
	B7	Red-edge 3	783 nm
	B8	Near Infrared	842 nm
	B8A	Red-edge 4	865 nm
	B11	Shortwave Infrared 1	1610 nm
	B12	Shortwave Infrared 2	2190 nm
* Seasonal windows defined based on key phenological periods within quaking aspen forests (see Section 2.2.4). Peak greenness (full canopy development), dormancy (complete canopy senescence), summer (mid greenup phase to onset of greenness decrease), and autumn (mid senescence phase to onset of greenness minimum).

Table 3. Spectral vegetation indices and radar textural features. The formula and associated reference to calculate each index is provide where relevant. Indices and textural features were calculated for each seasonal window.

Satellite	Index	Abbreviation	Formula	Seasonal Windows	Reference
Sentinel-1	GLCM Entropy	VV_ent, VH_ent	GLCM 7x7	Peak / Dormancy	[45]
Sentinel-1	GLCM Variance	VV_var, VH_var	GLCM 7x7	Peak / Dormancy
Sentinel-1	GLCM Correlation	VV_corr, VH_corr	GLCM 7x7	Peak / Dormancy
Sentinel-1	GLCM Contrast	VV_contrast, VH_contrast	GLCM 7x7	Peak / Dormancy
Sentinel-2	Chlorophyll Index Red-edge	CIRE	(B8 / B5) - 1	Summer / Autumn	[46]
Sentinel-2	Inverted Red-edge Chlorophyll Index	IRECI	(B8 - B4) / (B5 / B6)	Summer / Autumn	[47]
Sentinel-2	Specific Leaf Area Vegetation Index	SLAVI	B8 / (B4+ B12)	Summer / Autumn	[48]
Sentinel-2	Modified Chlorophyll Absorption in Reflectance Index	MCARI	((B5 - B4) - 0.2 * (B5 - B3)) * (B5 / B4)	Summer / Autumn	[33]
Sentinel-2	Red-edge Normalized Difference Vegetation Index	NDVI705	(B6 - B5) / (B6 + B5)	Summer / Autumn	[49]
Sentinel-2	Modified Normalized Difference Water Index	MNDWI	(B3 - B11) / (B3 + B11)	Summer / Autumn	[50]

Table 4. Precision, recall, and F1-score for aspen test data across the three reference datasets and for the Sentinel-based map.

Data Source	Total Area (Km²)	Number of Patches	Patch Density	Average Patch Size (ha)	Average Perimeter/Area Ratio
Sentinel-based Map	9,384.41	1,760,386	35.73	0.53	2,745.57
LANDFIRE EVT	13,441.59	728,370	14.22	1.85	1,060.98
USFS TreeMap	13,931.85	1,268,131	24.82	1.10	1,145.88
USFS ITSP	20,477.69	266,762	4.82	7.68	941.75

Table 5. Landscape patch dynamics summarized across spatial blocks for the Southern Rockies. We report the total area, number of patches, average patch density, average patch size, and average perimeter to area ratio.

Data Source	Total Area (Km²)	Number of Patches	Patch Density	Average Patch Size (ha)	Average Perimeter/Area Ratio
Sentinel-based Map	9,384.41	1,760,386	35.73	0.53	2,745.57
LANDFIRE EVT	13,441.59	728,370	14.22	1.85	1,060.98
USFS TreeMap	13,931.85	1,268,131	24.82	1.10	1,145.88
USFS ITSP	20,477.69	266,762	4.82	7.68	941.75

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer