Mapping Progression and Severity of a Southern Colorado Spruce Beetle Outbreak Using Calibrated Image Composites

An ongoing spruce beetle (Dendroctonus rufipennis Kirby.) epidemic in southern Colorado has resulted in the death of thousands of acres of forests primarily dominated by Engelmann spruce (Picea engelmannii Parry.). To evaluate the ecological and economic impacts of this massive mortality event, researchers and land managers need to efficiently track its progression, spread, and severity across large spatial extents. In this study, mortality severity (0–100% dead) was successfully mapped at the Landsat pixel scale (30 × 30 m) across a large (5000 km2), persistently cloud-covered study area using multi-sensor (Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI)) harmonized tasseled cap image composites as spectral predictors of gray stage spruce beetle mortality. Our maps display the distribution and severity of this landscape-scale mortality event in 2011 (R2 = 0.48, root mean squared error (RMSE) = 7.7) and 2015 (R2 = 0.55, RMSE = 11.6). Potential applications of this study include efficient landscape-scale forest health monitoring, targeted forest and timber management, and assessment of ecological impacts of bark beetle outbreaks.


Introduction
Timely monitoring of forest mortality across large spatial extents is a traditionally complex and resource-intensive endeavor [1][2][3].In the wake of a changing climate, multi-year drought conditions, and severe insect and disease activity, the complexity of monitoring and understanding change across vast expanses of public and privately owned forestlands in the western United States has only increased [4].A decade-long spruce beetle (Dendroctonus rufipennis Kirby.)outbreak in southern Colorado has killed Engelmann spruce (Picea engelmannii Parry.)trees across thousands of acres of spruce and fir dominated forests [5].With spruce/fir forest types representing approximately 20% of statewide forest cover [6] and encompassing the largest number of forested acres under public ownership within the state [7], monitoring and management of these forests are top priorities for both public and private land managers and owners.
The spruce beetle is one of the most damaging agents in mature spruce stands in Colorado.A native insect, the spruce beetle primarily attacks Engelmann spruce but can infest any spruce species found within the Colorado subalpine zone [8].Spruce beetles generally have a two-year life cycle; however, one-to three-year life cycles have been recorded [8,9].Adult female spruce beetles Forests 2018, 9, 336 2 of 14 bore into a tree and create galleries for their eggs in the tree's phloem tissue.Once the eggs hatch, the larvae overwinter in the galleries and eventually tunnel out of the tree, feeding on phloem tissue as they create additional galleries.These feeding galleries inhibit the flow of nutrients, weakening and eventually killing a tree [6].The beetles can be present at both endemic and outbreak (or "epidemic") population levels [8], typically beginning in areas that have experienced disturbances, often at sites affected by blowdown events or where woody debris have accumulated [10].At endemic levels, adult beetles attack downed woody material and debris.At outbreak levels, beetles will attack trees of all sizes; although large (>40 cm DBH) spruce trees are usually attacked first, trees of any size, including saplings can serve as suitable hosts as an outbreak progresses [6].Following infestation, the tree's needles slowly fade from green to yellowish-green until entering the gray phase where they ultimately drop.The entirety of this process typically takes less than two years [9,10].As beetle populations increase, the majority of suitable host trees within a stand can be killed [8,11].
In Colorado and the Intermountain West, resource managers currently rely on annual forest monitoring programs, such as the United States Forest Service (USFS) Aerial Detection Survey (ADS), to evaluate the impacts of bark beetle-induced tree mortality, and in turn, to plan and implement forest management projects [12], monitor forest carbon dynamics [13] and to help keep the public informed on the status and health of spruce/fir forest resources [14].While this program has been an important component of forest health monitoring in the United States, it is costly and measurements of insect-induced mortality, particularly the spatial extent, are not highly accurate and report mortality intensity at coarse spatial scales [15].As such, researchers have supplemented information provided by the USFS ADS program by remotely sensing bark beetle-induced tree mortality using moderate (30 × 30 m 2 ) and high resolution (<5 × 5 m 2 ) satellite imagery in combination with modelling [3].
Remote sensing has long been shown to be an effective method to detect mortality in coniferous forests [16].Many studies have employed moderate resolution satellite imagery, such as Landsat (30 × 30 m 2 ), to map the presence and absence of insect-induced canopy mortality [17][18][19].This is a relatively efficient and cost-effective method when compared to collecting similar data via aerial or field survey [3], but presence/absence maps of tree mortality do not convey the severity of an outbreak-information which is needed to accurately quantify spread, intensity, distribution, or ecosystem impacts of an infestation.
Recent research has improved methodologies for detecting mortality severity (or the percentage of an area with dead canopy present) at the stand and even single tree level [20] using remotely sensed imagery.In one such example, Long and Lawrence [21] focused on detecting mountain pine beetle (Dendroctonus ponderosae) induced tree mortality in Montana and demonstrated that it is possible to accurately detect mortality severity (0-100% dead) at the Landsat pixel scale.This method offers a promising opportunity for researchers to supplement information provided by aerial detection surveys with remotely sensed maps of canopy mortality severity.Applying these methods to larger study areas and at multiple time steps would allow for an understanding of the progression of outbreaks and gives a more complete representation of landscape dynamics.However, mapping mortality severity at larger scales and across time presents additional challenges.
When working with remotely sensed data to detect insect induced mortality, researchers must weigh the advantages and disadvantages of each available remotely sensed data product's spectral, temporal, and spatial resolution [20,22].Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) and Operational Land Imager's (OLI) spatial resolution of 30 × 30 m 2 has been shown to be of sufficient resolution for resource managers to explore spatial patterns and trends of bark beetle outbreaks across landscapes [19].Although Landsat's spectral resolution has been shown in many cases to be sufficient for used in predictive detection models of bark beetle attack, changing sensor properties through time (TM to ETM+ to OLI within our study period) complicates image selection and differencing and can result in spectral mismatches between neighboring image collections.
Landsat's temporal resolution of 16 days can limit image availability, especially in persistently cloudy areas such as Colorado's subalpine zone, and the Landsat 7 ETM+ scan line corrector failure [23] Forests 2018, 9, 336 3 of 14 further limits image availability.In some seasons it is possible that little or no cloud or scan line free Landsat imagery is available for a study area of interest [24].In many scenarios, this can be easily and acceptably rectified by selecting imagery that contain clouds only partially obscuring a study area of interest and excluding affected areas with a cloud mask.While this process is necessary in some applications, doing so leaves some areas unmonitored and limits the applicability of methods to other studies [3].These challenges with cloud cover and spectral inconsistency across sensors and space can be overcome through spectral harmonization and compositing of multiple images.This multi-step process is carried out by collecting all available Landsat imagery (including regularly discarded cloudy image collections) for a time period and modelling spectral relationships of images across sensors to ensure spectral alignment and spectral harmonization across time and space [25][26][27].The tasseled cap transformation is a particularly powerful transformation that can be employed in spectral harmonization and compositing procedures because of its reported spectral consistency through time and its cross-sensor application [28].Once an index of interest, like tasseled cap, is applied to available satellite collections and a spectral harmonization between sensors is applied, all available pixels from multiple dates are composited by taking the mean or median of all cloud and shadow-free observations for a given pixel [26,27].The resulting product is a spectrally aligned composited image that reduces cloudy patches to create the largest study area extent possible.Since gray stage spruce trees are relatively spectrally static, pixel-by-pixel compositing is a relevant way to increase the area sampled while concurrently maintaining the spectral information required to detect gray stage spruce canopy mortality.
The long-lived and spatially extensive spruce beetle outbreak in southern Colorado provides a unique opportunity to explore methods for detecting bark beetle mortality severity.The primary objective of this study was to map the severity and progression (2011-2015) of the spruce beetle outbreak occurring in southern Colorado spruce/fir forests.To achieve this objective, we utilized harmonized and composited Landsat indices to enhance our ability to monitor large, cloud-covered study areas.In addition, we tested the effectiveness of using composited Landsat indices to detect spruce beetle-induced outbreak severity at multiple time steps and of using multiple Landsat sensors.Existing methods to detect bark beetle canopy mortality severity are thought to be untested in these particular scenarios.

Study Area
We selected a c. 5000 km 2 study region (Figure 1) composed principally of Engelmann spruce (Picea engelmannii Kirby.) and subalpine fir (Abies lasiocarpa Nuttall.)located within and around the Rio Grande and San Juan National Forests in the southern Colorado Rocky Mountains.The study area was restricted to spruce/fir forest types using the publicly available LANDFIRE existing vegetation type layer [29].Burned areas were excluded from sampling using Monitoring Trends in Burn Severity [30] fire history records.Elevations in the study area range between ~1800 m and ~4000 m.While average temperatures and rainfall are quite variable in the study region, Elliot & Baker [31] averaged conditions reported at three nearby weather stations and found mean annual temperatures to be between −5.9 • C and 12.4 • C, with mean annual precipitation of 50.8 cm.Most of the study region is managed by the USFS for multiple-use objectives, including conservation of public lands, recreational activities, timber and resource extraction, as well as cattle grazing.
Aerial detection surveys show that spruce/fir forest types in the study region experienced varying levels of spruce beetle caused tree mortality over the past decade, which was reported at low levels beginning in 2004, later intensifying to outbreak proportions across the landscape between 2009 and 2015.The study area is an area of ecological interest because of the recent and intense nature of spruce beetle-induced tree mortality.This study region also emulates challenges encountered when researchers use satellite data to model ecological phenomena at the landscape level: (1) cloud cover obstructing areas of scientific interest, and spectral inconsistencies between imagery caused by (2) sensor mismatches and (3) working across multiple satellite collection paths.As such, the study area and period covered portions of two Landsat Worldwide Reference System 2 (WRS-2) path/rows (P/R), including P034, R034 and P035 R034, spanned multiple sensor periods (TM/ETM+ to OLI), and was distributed across high elevation spruce/fir forests that are often cloud covered.Aerial detection surveys show that spruce/fir forest types in the study region experienced varying levels of spruce beetle caused tree mortality over the past decade, which was reported at low levels beginning in 2004, later intensifying to outbreak proportions across the landscape between 2009 and 2015.The study area is an area of ecological interest because of the recent and intense nature of spruce beetle-induced tree mortality.This study region also emulates challenges encountered when researchers use satellite data to model ecological phenomena at the landscape level: (1) cloud cover obstructing areas of scientific interest, and spectral inconsistencies between imagery caused by ( 2) sensor mismatches and (3) working across multiple satellite collection paths.As such, the study area and period covered portions of two Landsat Worldwide Reference System 2 (WRS-2) path/rows (P/R), including P034, R034 and P035 R034, spanned multiple sensor periods (TM/ETM+ to OLI), and was distributed across high elevation spruce/fir forests that are often cloud covered.

Data Collection
We randomly distributed 410 sampling locations within the study area to facilitate ocular estimation of canopy mortality using high-resolution (1 m 2 ) National Agricultural Imagery Program (NAIP) imagery for 2011 and 2015.We delineated 30 × 30 m Landsat pixel boundaries at the locations where the 410 sampling points were located and overlaid a 100-square grid to aid in ocular estimation of four categories: percent gray stage canopy mortality, live canopy, other live vegetation, and "other" within each of the plots.The sampling strategy was designed so that estimates of percent canopy mortality would spatially coincide and be at the same spatial scale with the pixel-level spectral values extracted from Landsat imagery [21,32].While neither NAIP nor Landsat is geolocated with perfect precision, their geolocation accuracies have been shown to be sufficient to have a minimal effect on sampling accuracy, which was supported through ocular comparisons of the image sets.All ocular estimation of canopy mortality using NAIP imagery was conducted within Google Earth Engine's API [33] using a scripted interface that allowed for near instant mosaicking and display of NAIP orthorectified quarter quad tiles for all years of interest.

Data Collection
We randomly distributed 410 sampling locations within the study area to facilitate ocular estimation of canopy mortality using high-resolution (1 m 2 ) National Agricultural Imagery Program (NAIP) imagery for 2011 and 2015.We delineated 30 × 30 m Landsat pixel boundaries at the locations where the 410 sampling points were located and overlaid a 100-square grid to aid in ocular estimation of four categories: percent gray stage canopy mortality, live canopy, other live vegetation, and "other" within each of the plots.The sampling strategy was designed so that estimates of percent canopy mortality would spatially coincide and be at the same spatial scale with the pixel-level spectral values extracted from Landsat imagery [21,32].While neither NAIP nor Landsat is geolocated with perfect precision, their geolocation accuracies have been shown to be sufficient to have a minimal effect on sampling accuracy, which was supported through ocular comparisons of the image sets.All ocular estimation of canopy mortality using NAIP imagery was conducted within Google Earth Engine's API [33] using a scripted interface that allowed for near instant mosaicking and display of NAIP orthorectified quarter quad tiles for all years of interest.Ocular estimation was carried out by two calibrated image interpreters with multiple years of experience working together to conduct image interpretation.We classified tree mortality as "mortality" only if the tree was characteristic of spruce beetle-induced mortality in the gray stage.Other types of disturbance that were seen in the study area (i.e., windthrow, management activities) were not included in estimates of canopy mortality to ensure our estimates were characteristic of those resulting from spruce beetle attack.

Remotely Sensed Data
To derive predictors of spruce beetle outbreak severity, we obtained all available Landsat imagery for three sensors (TM, ETM+ and OLI) across two primary study periods: 2011 (July-August) and 2015 (July-August).We also obtained all available imagery for the same months in 2000 and 2007 to characterize pre-and mid-outbreak forest conditions, respectively.All products were obtained pre-processed to surface reflectance through the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On Demand Interface [34].Imagery was selected without consideration of cloud cover, but the study period was restricted to these two peak summer months to ensure the high elevation study area was free of snow and ice and deciduous vegetation was leaf-on so as to not be confused with gray stage spruce.All imagery was obtained in the Albers equal-area conic projection [35].
The primary predictive indices of focus in this study were Landsat-derived tasseled cap transformations [28,36] and derivations thereof.The tasseled cap transformation transforms spectral information from the six reflective Landsat bands into three interpretable bands directly associated with landscape and vegetative characteristics, including brightness (which represents soil and image brightness), greenness (which represents vegetative greenness) and wetness (which represents soil and vegetation wetness) [28].We selected tasseled cap transformations as the primary set of predictive indices because of their reported spectral consistency through time and their cross-sensor application [28], both of which are important when compositing multiple images.In addition, tasseled cap transformations have strong experimental precedence as robust predictors of forest disturbances and changes in canopy characteristics [37][38][39].
We processed all imagery using LandsatLinkr [26,27], an R package, to obtain near cloud-free tasseled cap brightness, greenness, and wetness (BGW) indices that were used as predictors of spruce beetle induced tree mortality for both 2011 and 2015 (Figure 2, Table 1).LandsatLinkr automatically conducts basic image preparation (decompression, stacking, reprojection), masks all cloud covered pixels within a Landsat scene using Fmask [40], and applies standard reflectance based tasseled cap transformations [36] to TM/ETM+ sensors.Landsat 8 indices are then created by integrating ETM+ tasseled cap indices with near-date OLI surface reflectance imagery to model spectral relationships and generate OLI tasseled cap indices, helping to ensure spectral consistency across sensors [25,27].The tool then composites unobstructed (cloud and shadow free) portions of all harmonized images using a mean of overlapping pixel values for all images available in the year, including those areas that overlap across Landsat WRS-2 scene boundaries.The resulting products are single growing season composites of tasseled cap BGW indices for each year (hereafter referred to as compositeTCAP products).A detailed description and workflow for the LandsatLinkr tool is available in Vogeler et al. [27].
We differenced the 2011 and 2015 BGW compositeTCAP with previous years BGW compositeTCAP products (2015-2011, 2011-2007 and 2015/2011-2000).These differenced layers captured changes in TCAP values between pre-, mid-, and post-spruce beetle outbreak conditions, a practice commonly employed in land and forest change studies [41].Finally, to capture the potential influence of topography on spruce beetle mortality characteristics, we derived elevation, slope, and aspect from the 1 arc-second Shuttle Radar Topography Mission v2.0 digital elevation model product.After all processing was completed, 12 individual predictive data layers were available for both 2011 and 2015 (Table 1).
Finally, to test the robustness and application of compositeTCAP, which used 16 individual image collections versus single image date tasseled cap derived indices, we created one additional set of TCAP BGW indices for 2015 that used only one image date per scene.The two scenes with the lowest cloud cover for the study period were downloaded pre-processed to surface reflectance, transformed to TCAP indices using coefficients specifically created for Landsat 8 OLI [42], and using the top layer values, were mosaicked for analysis and comparison to the compositeTCAP product (hereafter, this product is referred to as singleTCAP).

Mapping Spruce Mortality Severity in 2011 and 2015
We used a combination of the mortality severity data collected at plots in 2011 and 2015 and extracted associated spectral values from the twelve predictor variables (Table 1) to build two (2011 and 2015) random forest models of spruce mortality severity.Random forests are a robust decision tree based prediction algorithm that can be used in both classification and regression based problems [43].Random forests are a particularly powerful tool for remote sensing based analyses [44] because they are non-parametric and difficult to overfit [43,45].Random forests also facilitate simple evaluation of model performance without the use of separate testing data because the model is built with a randomly selected subset of the predictors and training data, resulting in "built in" cross validation in each model run [46].
Both models were tested and built in the R statistical software using the randomForest package [45] using 2000 trees and the remainder of parameters left at default.We conducted model selection using the rfUtilities package model selection function [47] to achieve a balance of model parsimony and predictive power.We evaluated the performance of each model using out-of-bag estimates of variance explained, root mean squared error (RMSE) and percent RMSE of the maximum observation.After selecting the best performing models for 2011 and 2015, we applied the models to generate spatial predictions of mortality severity across the entire study area.

Effectiveness of Composited TCAP Indices
To perform an initial test to determine whether compositeTCAP BGW indices were an effective substitute for traditional singleTCAP BGW in modelling spruce beetle induced tree mortality, we compared them in two ways.First, we compared the proportion of the total study area that remained following cloud and cloud shadow masking.This was completed by simply tabulating within a GIS the number of pixels remaining following completion of both processing types.Next, we built a simplistic random forest model for each image type and compared evaluation metrics.The models used the same plots and the same response variable.The only difference was the predictor variable used: compositeTCAP or singleTCAP (Table 2).We extracted values of the BGW predictors for the 342 plots (out of 410 original plots) that were still present after cloud masking in both the singleTCAP and compositeTCAP products and built a random forest model for each dataset.

Comparison of Composited TCAP Indices with Single TCAP Indices
The use of composited TCAP indices, which combined 16 tasseled cap transformed predictor layers of BGW from all available 2015 July and August Landsat imagery, increased the spruce/fir area that was clear of clouds and available for analysis.The use of the lowest cloud cover single image dates in the same area and time period reduced the study area size by 16% (Table 2).This does, however, represent a snapshot in time and serves only as an example of the power that compositing imagery can have to expand a study area in a situation where cloud-free imagery is simply unavailable.While many studies are forced to select suboptimal image dates or to significantly reduce their study area size because of cloud coverage, compositing allowed us to expand the area of analysis to cover nearly all (99.5%) spruce/fir forests within the study region, while maintaining the optimal time period of interest in the peak summer months.
Each of the random forest models performed similarly, with singleTCAP very slightly outperforming compositeTCAP.Since both models used the same set of observations and are based upon the same response variable, we can directly compare the two model's evaluation metrics.Root mean square error (RMSE), a measure of the differences between the observed and predicted values of the model, was similar for both at 14.9 (18%) for compositeTCAP and 13.9 (17%) for singleTCAP.This initial test gave us confidence in the use of compositeTCAP products in the more refined modelling techniques employed in the final models of spruce mortality, which used an expanded set of predictors and training dataset.We were satisfied by the increase in study area size (16%, ~900,000 pixels) produced through compositing and decided it a worthwhile tradeoff for the very minor reduction in model performance, thoughmore comparisons are needed in additional modelling scenarios.Gray stage spruce trees are likely spectrally stagnant within a summer season, so if composited imagery was to be used in an analysis more sensitive to phenological variations, we recommend additional testing of the suitability of composited imagery or indices.

Final Models of Spruce/Fir Mortality Severity in 2011/2015
The final models for both 2011 and 2015 performed quite well in predicting spruce mortality across the expansive and environmentally diverse study area.The 2011 and 2015 models had an RMSE of 7.7 (14%) and 11.6 (14%) and explained 48% and 55% of variance, respectively.The twelve predictors for each year were pared down to six in both models after variable selection (Figure 3).The predictor variables retained and their importance only varied slightly between 2011 and 2015.Brightness and differenced wetness were the top predictor variables in both the 2011 and 2015 models.The inclusion of topographic indices provided no additional predictive power to the model in 2011 but, interestingly, elevation contributed in 2015.
cover single image dates in the same area and time period reduced the study area size by 16% (Table 2).
This does, however, represent a snapshot in time and serves only as an example of the power that compositing imagery can have to expand a study area in a situation where cloud-free imagery is simply unavailable.While many studies are forced to select suboptimal image dates or to significantly reduce their study area size because of cloud coverage, compositing allowed us to expand the area of analysis to cover nearly all (99.5%) spruce/fir forests within the study region, while maintaining the optimal time period of interest in the peak summer months.
Each of the random forest models performed similarly, with singleTCAP very slightly outperforming compositeTCAP.Since both models used the same set of observations and are based upon the same response variable, we can directly compare the two model's evaluation metrics.Root mean square error (RMSE), a measure of the differences between the observed and predicted values of the model, was similar for both at 14.9 (18%) for compositeTCAP and 13.9 (17%) for singleTCAP.This initial test gave us confidence in the use of compositeTCAP products in the more refined modelling techniques employed in the final models of spruce mortality, which used an expanded set of predictors and training dataset.We were satisfied by the increase in study area size (16%, ~900,000 pixels) produced through compositing and decided it a worthwhile tradeoff for the very minor reduction in model performance, thoughmore comparisons are needed in additional modelling scenarios.Gray stage spruce trees are likely spectrally stagnant within a summer season, so if composited imagery was to be used in an analysis more sensitive to phenological variations, we recommend additional testing of the suitability of composited imagery or indices.

Final Models of Spruce/Fir Mortality Severity in 2011/2015
The final models for both 2011 and 2015 performed quite well in predicting spruce mortality across the expansive and environmentally diverse study area.The 2011 and 2015 models had an RMSE of 7.7 (14%) and 11.6 (14%) and explained 48% and 55% of variance, respectively.The twelve predictors for each year were pared down to six in both models after variable selection (Figure 3).The predictor variables retained and their importance only varied slightly between 2011 and 2015.Brightness and differenced wetness were the top predictor variables in both the 2011 and 2015 models.The inclusion of topographic indices provided no additional predictive power to the model in 2011 but, interestingly, elevation contributed in 2015.The 2011 spruce mortality map shows the highest severity outbreaks were located in the northeastern corner of the study area, between San Luis Peak and Wolf Creek Pass and just west of this location (Figure 4).Other low severity outbreaks were detected across fairly small spatial extents, but the vast majority of the study region appeared to have been unaffected or affected at low levels of severity (10-20% dead) in 2011 (Figure 5).Between the 2011 and 2015 study periods, spruce beetle    It is important to note that the study period included only the height of the spruce beetle outbreak (2011-2015) occurring in the region.To understand the full history and progression of spruce beetle activity within this outbreak, additional sampling from previous years would be required, although this comes with an additional caveat.This method is successful because we have a wide range of mortality severity occurring across the landscape.Modelling mortality severity at the beginning of the outbreak could be less successful because a randomly distributed dataset would contain mostly absence values, which may be better suited for a classification modelling approach.
Finally, while we took care to exclude disturbances that were not caused by spruce beetle in sampling and initial pre-processing steps (such as excluding areas affected by fire and by restricting the study area to the spruce/fir vegetation class) we are almost certainly capturing disturbances in the final maps that are not induced by spruce beetle.Forest canopies in the area were observed to be almost entirely green in scans of pre-outbreak imagery, but other insects, tree diseases, and disturbances are present in the area, which may be classified by the model as gray stage spruce mortality.Because of the massive extent and high severity of the ongoing outbreak, we believe it is reasonable to assume that the vast majority of predicted mortality was the result of spruce beetle attack.By 2015, the spruce beetle attack seemed to expand in all directions from the relatively spatially isolated mortality events seen in 2011, with areas along the northern border of the Rio Grande National Forest being attacked very severely (50-70% dead).The outbreak occurring near Wolf Creek Pass increased in severity, and an additional high severity mortality pocket appeared in the southeastern portion of the study area.The outbreak dynamics conveyed in these maps are consistent with findings in other studies showing that spruce beetle outbreaks spread from multiple population centers instead of from a single "epicenter" [48].By 2015, outbreaks converged and very little of the study region had been left unaffected by spruce beetle attack (Figures 4 and 5).These maps are consistent with the broad patterns of new attacks conveyed through the USFS ADS maps (Figure 6).This lends credence to the concept that maps such as those produced in this study could be combined with ADS data to provide an enhanced product that more clearly and precisely conveys information on the severity and distribution of insect and disease activity.
This study shows that multiyear and multi-sensor Landsat imagery is effective for monitoring large, persistently cloudy areas.While many studies have focused upon a single year of monitoring, we have mapped the progression and severity of a bark beetle outbreak across a five-year period using multi-sensor tasseled cap composites.Updates to this progression can be efficiently reproduced with little to no field work whenever high-resolution aerial imagery is available.We further show that the methods are effective across spatially expansive (5000 km 2 ) study regions, a larger size than explored in previous studies using similar methods [21,49].In addition, the results show that composited tasseled cap indices are effective in detecting spectrally subtle, slow-moving disturbances spanning low to high canopy mortality.While similar predictive indices have been used in studies exploring deforestation [50], bark beetle induced disturbances are often considered more difficult to detect than more severe and fast moving disturbances [18], such as fire or forest harvest, when using moderate resolution imagery like Landsat.Composited tasseled cap indices and their derivatives were effective predictors of gray stage spruce mortality when applied within a regression tree based modelling framework.The variable selection procedure highlighted the important role that differencing indices from different stages of an outbreak can play in detecting changes in vegetation through time, a finding that is consistent with similar studies [41].While the 2011 and 2015 models both slightly overpredict low levels of mortality and underpredict more severe levels of mortality (Figure S1), we believe the RMSE is within an acceptable range for the models to be valuable tools to be applied in forest management and in evaluating ecosystem impacts of spruce beetle outbreaks.The resulting maps represent a significant improvement over products that display only presence/absence of canopy mortality.
It is important to note that the study period included only the height of the spruce beetle outbreak (2011-2015) occurring in the region.To understand the full history and progression of spruce beetle activity within this outbreak, additional sampling from previous years would be required, although this comes with an additional caveat.This method is successful because we have a wide range of mortality severity occurring across the landscape.Modelling mortality severity at the beginning of the outbreak could be less successful because a randomly distributed dataset would contain mostly absence values, which may be better suited for a classification modelling approach.
Finally, while we took care to exclude disturbances that were not caused by spruce beetle in sampling and initial pre-processing steps (such as excluding areas affected by fire and by restricting the study area to the spruce/fir vegetation class) we are almost certainly capturing disturbances in the final maps that are not induced by spruce beetle.Forest canopies in the area were observed to be almost entirely green in scans of pre-outbreak imagery, but other insects, tree diseases, and disturbances are present in the area, which may be classified by the model as gray stage spruce mortality.Because of the massive extent and high severity of the ongoing outbreak, we believe it is reasonable to assume that the vast majority of predicted mortality was the result of spruce beetle attack.

Conclusions
This study displayed an effective methodology for detecting gray stage spruce mortality in a large, cloudy study region spanning multiple Landsat scenes.We have shown that (1) multi-image date and multi-sensor tasseled cap composites are a powerful tool when working in a cloudy study area and serve as an effective predictor of gray stage spruce mortality severity across time and in large study regions; (2) differencing images from multiple time steps is important when attempting to detect long, slow moving disturbances like spruce beetle outbreaks, and (3) data produced through models like those in this study can serve as a potential supplement to existing forest management and monitoring programs, such as the USFS Aerial Detection Survey.Finally, we created a map of an ongoing outbreak event that supports a better understanding of outbreak progression, spread, and intensity for this particular spruce beetle outbreak.

Forests 2018, 9 , 15 Figure 1 .
Figure 1.The study area located in and around San Juan and Rio Grande National Forests.

Figure 1 .
Figure 1.The study area located in and around San Juan and Rio Grande National Forests.

Figure 3 .
Figure 3. Radar charts displaying the variable importance of predictors retained in the final models for (A) 2011 and (B) 2015 models.

Figure 3 .
Figure 3. Radar charts displaying the variable importance of predictors retained in the final models for (A) 2011 and (B) 2015 models.

Forests
around the San Juan and Rio Grande National Forests became more severe, with a much wider spatial distribution of outbreak events.

Forests 2018, 9 , 15 Figure 4 .
Figure 4. Final predictions of percent gray stage canopy mortality in and around Rio Grande and San Juan National Forests in 2011 and 2015.Figure 4. Final predictions of percent gray stage canopy mortality in and around Rio Grande and San Juan National Forests in 2011 and 2015.

Figure 4 .
Figure 4. Final predictions of percent gray stage canopy mortality in and around Rio Grande and San Juan National Forests in 2011 and 2015.Figure 4. Final predictions of percent gray stage canopy mortality in and around Rio Grande and San Juan National Forests in 2011 and 2015.

ForestsFigure 5 .
Figure 5. Relative area of percent mortality as the outbreak progressed between 2011 and 2015.Mortality intensity increased as the outbreak progressed.

Figure 5 .
Figure 5. Relative area of percent mortality as the outbreak progressed between 2011 and 2015.Mortality intensity increased as the outbreak progressed.

Figure 6 .
Figure 6.Spruce beetle-induced mortality as documented by the USFS Aerial Detection Survey (ADS) between 2004 and 2015.There is a close match between the general location and timing of mortality events as documented by ADS and our models.Data provided by Pete Barry, Colorado State Forest Service.

:
Plots of predicted values versus observed values for the 2011 and 2015 final models.Author Contributions: B.D.W., P.H.E. and A.G.V. conceived, designed, and coordinated the research; B.D.W. performed the research; B.W. and A.G.V. analyzed the data; P.H.E.contributed to data collection and organization.B.D.W. led the writing of the paper, with assistance from all co-authors.

Table 2 .
Characteristics of the two predictive layer stacks that were used in the initial comparison of composited predictors versus single image date predictors and the number of cloud free pixels remaining in each product following the two processing schemes.