Object-Based Classification of Forest Disturbance Types in the Conterminous United States

Forest ecosystems provide critical ecosystem goods and services, and any disturbance-induced changes can have cascading impacts on natural processes and human socioeconomic systems. Forest disturbance frequency, intensity, and spatial and temporal scale can be altered by changes in climate and human activity, but without baseline forest disturbance data, it is impossible to quantify the magnitude and extent of these changes. Methodologies for quantifying forest cover change have been developed at the regional-to-global scale via several approaches that utilize data from high (e.g., IKONOS, Quickbird), moderate (e.g., Landsat) and coarse (e.g., Moderate Resolution Imaging Spectroradiometer (MODIS)) spatial resolution satellite imagery. While detection and quantification of forest cover change is an important first step, attribution of disturbance type is critical missing information for establishing baseline data and effective land management policy. The objective here was to prototype and test a semi-automated methodology for characterizing high-magnitude (>50% forest cover loss) forest disturbance agents (stress, fire, stem removal) across the conterminous United States (CONUS) from 2003–2011 using the existing University of Maryland Landsat-based Global Forest Change Product and Web-Enabled Landsat Data (WELD). The Forest Cover Change maps were segmented into objects based on temporal and spatial adjacency, and object-level spectral metrics were calculated based on WELD reflectance time series. A training set of objects with known disturbance type was developed via high-resolution imagery and expert interpretation, ingested into a Random Forest classifier, which was then used to attribute disturbance type to all 15,179,430 forest loss objects across CONUS. Accuracy assessments of the resulting classification was conducted with an independent dataset consisting of 4156 forest loss objects. Overall accuracy was 88.1%, with the highest omission and commission errors observed for fire (32.8%) and stress (31.9%) disturbances, respectively. Of the total 172,686 km2 of forest loss, 83.75% was attributed to stem removal, 10.92% to fire and 5.33% to stress. The semi-automated approach described in this paper provides a promising framework for the systematic characterization and monitoring of forest disturbance regimes.


Introduction
Forests provide critical ecosystem goods and services and are deeply intertwined with socioeconomic systems [1,2].Consequently, any disturbance-induced change in forest cover can have cascading impacts on surface energy balances, carbon dynamics, wildlife habitat, and human activities [3][4][5][6].In the United States, wildfires, insect outbreaks, and land use (e.g., timber harvest, urbanization) are a few of the major disturbances that can lead to large-scale changes in forest cover [7][8][9][10][11][12][13][14].Disturbance regimes (frequency, intensity, spatial and temporal scale) vary significantly within and between disturbance types, making widespread, automated characterization on year of forest loss and spatial adjacency.Object-level spectral trajectories (including spectral indices and temporal profile metrics) were identified by adapting the approach proposed by Kennedy et al. [39], using web-enabled Landsat Data [40].Visual interpretation of high-resolution Google Earth imagery, in combination with ancillary datasets, was used to create training and validation datasets of objects with known disturbance type.A Random Forest classifier was then used to classify all disturbance objects.The classification was then refined using the Moderate Resolution Imaging Spectroradiometer (MODIS) Active Fire Product [41].The classified forest loss objects were compared to an independent validation dataset and accuracy was assessed using error matrices.

University of Maryland Global Forest Change Product
The Global Forest Change Product, produced by University of Maryland, provides 30 m spatial resolution global tree canopy cover extent and forest cover change (i.e., forest cover loss and forest cover gain) from 2000-2016 [23,42].The forest loss product was designed to capture stand replacing events and may not detect low-to-moderate-intensity disturbances that result in partial forest cover loss.Tree canopy cover and change datasets are distributed in 10 × 10 degree tiles.User's accuracy for the forest loss product in temperate forests was 88.2%, Producer's accuracy was 93.9%, and overall accuracy was 99.8% [20].Accuracy is likely to be higher in CONUS than other regions due to the higher number of available Landsat scenes [43].

Web Enabled Landsat Data (WELD)
The Web Enabled Landsat Data (WELD) dataset consists of temporally composited mosaics of CONUS and Alaska using Landsat 7 Enhanced Thematic Mapper Plus (ETM+) Level 1 Terrain Corrected (L1T) data.All ETM+ pixels with cloud cover ≤ 80% are considered for inclusion in weekly, monthly, or annual composites [40].The Level 1T data used to generate the WELD product is processed using radiometric correction, systematic geometric correction, precision correction with ground control chips, and correction of parallax error aided by digital elevation models [37].The WELD product includes top of atmosphere (TOA) reflectance (bands 1-5, 7) and brightness temperature (band 6), TOA normalized difference vegetation index (NDVI), and the date of pixel acquisition.WELD version 1.5 Products were used in this study.The WELD products are designed to provide consistent data that can be easily used to derive land cover and geo-physical and bio-physical products, such as land cover change [44], and cropland delineation [45,46].Image composites are distributed in 5000×5000 30 m pixel tiles in the Albers equal area projection.In this study, annual CONUS WELD mosaics from 2002-2012 were used.TOA reflectance for ETM+ bands 1 (0.45-0.52 µm), 2 (0.53-0.61 µm), 3 (0.63-0.69 µm), 4 (0.78-0.90 µm), 5 (1.55-1.75µm), and 7 (2.09-2.35µm) from WELD mosaics were used to calculate object-level spectral metrics for each forest cover disturbance (see Section 3.4).

Moderate Resolution Imaging Spectroradiometer (MODIS) Active Fire Product
The MODIS active fire product detects fires that are actively burning at the time of NASA satellite (Terra and Aqua) overpass, with up to four overpasses possible per day at CONUS latitude [41].The Collection 5 Level 3 8-day MODIS Terra active fire product, MOD14A1 was used for this study.The MOD14A1 product contains 8 layers indicating the location of daily active fire detections at a nominal 1 km spatial resolution, includes ancillary detection confident information, and is distributed in the native MODIS sinusoidal projection [47].Annual active fire detection files were created by temporally compositing MOD14A1 tiles for each year from 2003-2012.The compositing procedure is described in detail [48]; to provide the greatest number of possible active fire detections, all confidence levels were used.These composites were used in the post-classification phase, to refine forest cover loss objects classified as fire-caused, by confirming that they were close in space and time to MOD14A1 detections.

Insect Mortality-Area Dataset
The insect mortality-area dataset developed by Meddens et al. [31] provides a 1 km gridded bark beetle-caused tree mortality-area product for British Columbia and the western United States from 1997-2010.This dataset was created by converting aerial detection survey (ADS) polygon data collected by the USDA Forest Health Monitoring Program in United States, to percent crown mortality per unit area.This product includes percent crown mortality resulting from nine major bark beetle species (e.g., mountain pine beetle, pinyon ips, ips engraver beetle, western balsam bark beetle, fir engraver, Douglas-fir beetle, spruce beetle, western pine beetle, and unspecified bark beetles) in the western United States and British Columbia, Canada.In this study, the insect mortality-area dataset was used solely as an ancillary dataset to inform the selection of training and validation dataset locations.

Monitoring Trends in Burned Severity (MTBS) Burned Area Perimeters
The MTBS project produces annual burned area polygons and fire severity maps for CONUS and Alaska from 1984 to the present date.Burned area perimeters are created via on-screen expert interpretation of Landsat data and U.S. Forest Service fire incident information when available [30].Only large fires are mapped (>400 ha in western U.S. and >200 ha in eastern U.S.).Like the insect detection data, MTBS fire perimeters were used solely as ancillary dataset to inform the selection of training and validation data locations.

GoogleEarth High Spatial-Resolution Imagery
GoogleEarth imagery consists of high spatial resolution (typically <5 m) data acquired by U.S. government programs, such as the National Agricultural Imagery Program (NAIP), and commercial companies, such as DigitalGlobe.This imagery allows for the interpretation of tree crown absence/presence and was used in the production of training and validation datasets (see Sections 3.5 and 3.7 for more detail).

Data Preprocessing
Following [48], all datasets were projected to the WELD Albers equal area projection (datum: WGS84).The 1 km datasets (insect mortality-area dataset and MOD14A1) were reprojected and resampled to 30 m spatial resolution.As mentioned previously, the insect detection dataset was used only as a guide to find suitable training and validation dataset locations, and MOD14A1 was used to confirm that objects classified as fire were close in space and time relative to MOD14A1 detections.Consequently, geolocation errors arising during reprojection and resampling are not critical for this analysis.

Image Segmentation
The reprojected Global Forest Change Product was segmented based on year of disturbance and spatial adjacency using Trimble eCognition Developer 9.1 software.Specifically, forest loss pixels were partitioned into objects using a 4-neighbor adjacency rule: (1) if they occurred in the same year, and (2) if they shared at least one edge with another forest cover loss pixel of the same year.Single-pixel objects were analyzed in the same manner as larger objects.A 4-neighbor adjacency rule was chosen to minimize the risk of forest cover loss objects caused by different disturbances being segmented into the same object.The resulting 15,179,430 forest cover loss objects became the functional units for all subsequent analyses.To match the temporal extent of the WELD dataset, only forest cover loss from 2003 to 2011 was considered.

Temporal Consistency of Forest Loss Year
The WELD annual datasets are generated from all the available Landsat scenes of the year by applying a maximum normalized difference vegetation index (NDVI) compositing algorithm [40].For pixels that experience a forest canopy loss (or more generally a disturbance), either a pre-disturbance or a post-disturbance value might be selected in the composite of the disturbance year depending on the nature of the spectral changes and on the amount of cloud-free data available before and after the disturbance.In some cases, this causes a discrepancy between the forest cover loss year estimated by the Global Forest Change Product, and the WELD time series.For each forest cover loss object defined in Section 3.2, the mean normalized burned ratio (NBR) was computed using the WELD composites from two years before the disturbance year, to the year after the disturbance year.If the disturbance year did not coincide with the year with the largest NBR drop, the year with the largest NBR drop was used in the subsequent steps of the analysis; in 5.6% of the objects the largest NBR drop occurred the year before the forest cover loss year, 23.2% of the objects the largest NBR drop occurred the year after the forest cover loss year.This was a purely nominal change intended to deal with the artifacts of the compositing algorithm: after the disturbance classification, all the objects were reassigned to the original disturbance year.

Object-Level Spectral Indices and Temporal Trajectories
Spectral indices were calculated for each forest loss object each year during the study period (2003-2011) using co-located WELD reflectance data.Spectral indices were first calculated at the pixel level and then averaged to the object.Any pixels labeled as cloudy in the WELD product were removed prior to any object-level spectral index calculation.The average cloud cover in the annual datasets for CONUS during the study period was 0.16%.In the rare instances where an object was comprised entirely of cloudy pixels, a mean value of the year before and after was assigned to the object.Based on previous land cover change studies (e.g., [33,[49][50][51]), the following spectral indices were calculated: normalized difference vegetation index (NDVI), red-green index (RGI), band5 (B5), the Band 5/Band 4 ratio (B54R), normalized difference moisture index (NDMI), normalized burn ratio (NBR), tasseled cap brightness (BRI), tasseled cap greenness (GRE) and tasseled cap wetness (WET) (Table 1).ρ red = red reflectance, ρ green = green reflectance, ρ NIR = near infrared reflectance, ρ SWIR = short − wave infrared reflectance, ρ blue = blue reflectance.
To capture the temporal spectral properties of each forest disturbance type, spectral index data for all 15,179,430 forest loss objects was processed implementing a similar temporal segmentation strategy as presented in Kennedy et al. [39].Spectral trajectories were decomposed to a sequence of straight-line segments that capture slow changes (e.g., forest recovery after disturbance) and abrupt changes (e.g., forest loss due to stand replacing fire).Specifically, a linear regression fit of spectral index versus year was computed using all points and the first and last years as vertices.The point with the largest deviation from the fit became the next vertex and regression fits were computed between this vertex and the first and last data points in the time series (resulting in two segments).This computation was repeated until four segments were fit to the spectral trajectories.This number was assumed to capture the typical forest disturbance lifecycle: pre-disturbance stable forest, forest disturbance, forest recovery, and post-disturbance stable forest.Figure 1 shows typical spectral profiles for stem removal (h03v04 tile, Oregon; WELD tile system displayed in Figure 3), fire (h02v06 tile, California), and stress (h11v09 tile, Colorado).Detection differences between the spectral indices are clearly apparent as well as the temporal characteristics of each disturbance.For example, forest stem removal and fire have abrupt changes in spectral index values (Figure 1, top and middle rows), whereas stress damage is characterized by gradual change over time (Figure 1, bottom row).
To capture the temporal spectral properties of each forest disturbance type, spectral index data for all 15,179,430 forest loss objects was processed implementing a similar temporal segmentation strategy as presented in Kennedy et al. [39].Spectral trajectories were decomposed to a sequence of straight-line segments that capture slow changes (e.g., forest recovery after disturbance) and abrupt changes (e.g., forest loss due to stand replacing fire).Specifically, a linear regression fit of spectral index versus year was computed using all points and the first and last years as vertices.The point with the largest deviation from the fit became the next vertex and regression fits were computed between this vertex and the first and last data points in the time series (resulting in two segments).This computation was repeated until four segments were fit to the spectral trajectories.This number was assumed to capture the typical forest disturbance lifecycle: pre-disturbance stable forest, forest disturbance, forest recovery, and post-disturbance stable forest.Figure 1 shows typical spectral profiles for stem removal (h03v04 tile, Oregon; WELD tile system displayed in Figure 3), fire (h02v06 tile, California), and stress (h11v09 tile, Colorado).Detection differences between the spectral indices are clearly apparent as well as the temporal characteristics of each disturbance.For example, forest stem removal and fire have abrupt changes in spectral index values (Figure 1, top and middle rows), whereas stress damage is characterized by gradual change over time (Figure 1, bottom row).For each spectral index, five features were extracted from the temporal profile, as shown in Figure 2. Average pre-disturbance spectral index values for the two years before the forest cover loss, average post-disturbance values for the disturbance year and the following year, and the difference between post-and pre-disturbance values were extracted.Note that, for the forest loss objects with forest loss year 2003, the pre-disturbance value was calculated from a single year of WELD data; conversely, for the forest loss objects of year 2011, only 2012 WELD data was used as post-disturbance value.The maximum and minimum slopes of the four segments fitted to each spectral index temporal trajectory were added to the predictor variables.The pre-disturbance value, post-disturbance value, difference value, maximum slope and minimum slope were denoted as (with NDVI spectral index as an example) bNDVI, aNDVI, delta_NDVI, maxNDVI, and minNDVI, respectively.The five features extracted for each of the nine spectral indices resulted in a total of 45 predictor variables to be used in the classification (Table 2).
Remote Sens. 2018, 1, x FOR PEER REVIEW 7 of 24 For each spectral index, five features were extracted from the temporal profile, as shown in Figure 2. Average pre-disturbance spectral index values for the two years before the forest cover loss, average post-disturbance values for the disturbance year and the following year, and the difference between post-and pre-disturbance values were extracted.Note that, for the forest loss objects with forest loss year 2003, the pre-disturbance value was calculated from a single year of WELD data; conversely, for the forest loss objects of year 2011, only 2012 WELD data was used as post-disturbance value.The maximum and minimum slopes of the four segments fitted to each spectral index temporal trajectory were added to the predictor variables.The pre-disturbance value, post-disturbance value, difference value, maximum slope and minimum slope were denoted as (with NDVI spectral index as an example) bNDVI, aNDVI, delta_NDVI, maxNDVI, and minNDVI, respectively.The five features extracted for each of the nine spectral indices resulted in a total of 45 predictor variables to be used in the classification (Table 2).

Training Dataset
To create a training dataset for the Random Forest classifier, a stratified sampling approach was employed (Figure 3).Twenty-one 150 × 150 km WELD tiles were selected based on the following criteria: (i) at least one tile should be selected in each EPA Level II Ecoregion [56], and (ii) selected tiles should avoid overlapping multiple ecoregions while also encompassing as many forest disturbance objects as possible.Twelve WELD tiles were selected for stem removal training samples, based on documented timber extraction activity [57].The five tiles with the highest reported tree mortality due to insects in the 2001-2010 period, as reported by the ADS dataset [31], were used to collect stress disturbance training objects.Similarly, the eight tiles with the highest fire occurrence based on the MTBS dataset were used to collect fire training objects (four of those coincided with tiles also used to collect stem removal training objects).
disturbance objects as possible.Twelve WELD tiles were selected for stem removal training samples, based on documented timber extraction activity [57].The five tiles with the highest reported tree mortality due to insects in the 2001-2010 period, as reported by the ADS dataset [31], were used to collect stress disturbance training objects.Similarly, the eight tiles with the highest fire occurrence based on the MTBS dataset were used to collect fire training objects (four of those coincided with tiles also used to collect stem removal training objects).
In all cases, all objects greater than 5 pixels in size were sorted randomly, and the interpreter proceed to interpret them in order.The first 100 objects confidently interpreted as stem removal, insect damage and fire -respectively in the tiles designated to collect training of each classconstituted the training set.The exception were two tiles (h01v06, h07v03) where only 50 stem removal training objects could be identified, and four tiles (h01v06, h07v03, h07v13, h09v04) where only 50 fire objects could be identified.As a result, the training dataset included 1,100 stem removal, 600 fire, and 500 stress forest cover loss objects.The visual interpretation of the training objects was conducted using a methodology adapted from Cohen et al. [58].A trained interpreter relied primarily on high spatial resolution Google Earth imagery, using the MTBS polygons and insect disturbance data as ancillary source of information, with limited input from WELD multi-spectral data and temporal profiles, to manually classify each target forest loss object.It is important to note that Cohen et al. [58] used multiple interpreters to help control measurement error, whereas this study did not.However, like Cohen et al. [58], independent disturbance datasets (e.g., MTBS, insect disturbance data) were employed to aid the interpreter's decisions.In all cases, all objects greater than 5 pixels in size were sorted randomly, and the interpreter proceed to interpret them in order.The first 100 objects confidently interpreted as stem removal, insect damage and fire-respectively in the tiles designated to collect training of each class-constituted the training set.The exception were two tiles (h01v06, h07v03) where only 50 stem removal training objects could be identified, and four tiles (h01v06, h07v03, h07v13, h09v04) where only 50 fire objects could be identified.As a result, the training dataset included 1100 stem removal, 600 fire, and 500 stress forest cover loss objects.The visual interpretation of the training objects was conducted using a methodology adapted from Cohen et al. [58].A trained interpreter relied primarily on high spatial resolution Google Earth imagery, using the MTBS polygons and insect disturbance data as ancillary source of information, with limited input from WELD multi-spectral data and temporal profiles, to manually classify each target forest loss object.It is important to note that Cohen et al. [58] used multiple interpreters to help control measurement error, whereas this study did not.However, like Cohen et al. [58], independent disturbance datasets (e.g., MTBS, insect disturbance data) were employed to aid the interpreter's decisions.

Random Forest Classification
The Random Forest (RF) classifier was used to attribute disturbance type to the forest loss objects and was chosen due to its high performance in remote sensing image classification [59] and large area mapping.The RF classifier uses an ensemble approach, where bootstrap samples of training data are used to create a set, or ensemble, of classification trees [60].Forest loss objects are assigned a disturbance type by majority vote of the ensemble of trees.The R software package "randomForest" was used to implement the RF classifier [61].In the RF classifier, two parameters need to be set: the number of decision trees (Ntree) and the number of variables to be selected for the best split when growing the trees (Mtry).Ntree is commonly set to 500 because the errors stabilize before this number of trees is achieved [59].The Mtry parameter is usually set to the square root of the number of classification variables [62].Following these recommendations, Ntree was set 500, and Mtry was set to 6.In addition to the classification results, predicted class probabilities for each forest loss object are computed, which gives a measure of attribution confidence.Input variable importance was calculated as the mean decrease in classification accuracy.Some of the predictor variables were strongly correlated with each other (~4% of variable pairs had Pearson's correlation coefficient >0.75) and could bias variable importance if not taken into account [63].Consequently, we used a conditional permutation methodology [63] implemented in the "party" R package [64] to calculate mean decrease in classification accuracy.Variable importance can be useful for model reduction (e.g., using only the "important" variables to build simpler, more readily interpretable models), especially in the case where models have hundreds to thousands of predictor variables [65].Since the number of predictor variables in this study is relatively small (n = 45), all variables were input into the Random Forest classification.

Post-Classification of Fire Disturbance
Forest loss objects classified as fire were checked with the MODIS Active Fire Product for temporal and spatial consistency.If objects classified as fire were spatially and temporally co-located with MODIS active fire detections from the forest loss year or previous year, then the object retained its fire designation.Otherwise, the object was reclassified as the disturbance type with the second largest majority vote among the ensemble of classification trees.An object was considered spatially co-located with an active fire detection if at least 20% of its area overlapped at least one active fire detection.

Validation and Accuracy Assessment
We used a stratified cluster sampling strategy to select validation WELD tiles.Following Cohen et al. [13], sampling strata were produced by merging EPA (Environmental Protection Agency) Level II ecoregion data [56] into generalized strata based on common forest types.The resulting regional strata included Mountain West, Lowland West, Central, Northeast, and Southeast regions (Figure 4).In total, 81 WELD tiles were selected for a validation dataset.Selection of tiles within the five strata utilized a weighted approach, where strata with more forested area and a higher number of forest loss objects received more validation tiles.This strategy produced 19, 9, 6, 23, and 24 WELD validation tiles for the Mountain West, Lowland West, Central, Northeast, and Southeast strata, respectively (Figure 4).Within each of the 81 WELD validation tiles, 50 forest loss objects (>5 pixels in size) were randomly selected and manually classified into one of the three disturbance types using the same methodology as the training dataset; replacement samples were drawn if no disturbance could be confidently detected in the high resolution imagery.The final validation set consisted of 3029 stem removal (299,107 pixels), 595 fire (164,932 pixels) and 426 stress (22,320 pixels) forest loss objects, and 106 objects (2432 pixels) that could not be confidently labeled as disturbed.Disturbance type classification accuracy was quantitatively assessed using confusion matrices calculated between the validation dataset and the Random Forest classification result.Three commonly used classification measures (overall accuracy, producer's accuracy, and user's accuracy) are reported.

Image Segmentation
A total of 15,179,430 forest loss objects were segmented for the 415 CONUS WELD tiles with forest loss over the study period.The average forest loss object was ~3.87 hectares in size, however, there was considerable variability within and between forest disturbance types.Representative forest loss object size, shape, and pattern are shown for each major disturbance type in Figure 5.While a large proportion of forest loss objects (71.98%) were less than 5 pixels (Figure 6), they accounted for a small proportion of the total forest loss area (~9.5%).The number of image objects with object size ranging from 1 to 4 pixels accounted for 3.44%, 2.84%, 1.65%, 1.60% of the total forest loss area in the CONUS, respectively.
There was a total of 4,038,889 image objects that were greater than 5 pixels in size.At the regional level (see definition in Section 3.7 and Figure 4), the Southeast accounted for more than half of the total image objects greater than 5 pixels in size (50.22%) and had the largest average segment size (4.41 hectares) among the five regions.The Mountain West accounted for the third largest proportion of objects (up to 20.16%) with an average object size of 4.23 hectares.The Northeast ranked second in the number of objects (up to 21.86%), and objects averaged 2.7 hectares in size.Disturbance type classification accuracy was quantitatively assessed using confusion matrices calculated between the validation dataset and the Random Forest classification result.Three commonly used classification measures (overall accuracy, producer's accuracy, and user's accuracy) are reported.

Image Segmentation
A total of 15,179,430 forest loss objects were segmented for the 415 CONUS WELD tiles with forest loss over the study period.The average forest loss object was ~3.87 hectares in size, however, there was considerable variability within and between forest disturbance types.Representative forest loss object size, shape, and pattern are shown for each major disturbance type in Figure 5.While a large proportion of forest loss objects (71.98%) were less than 5 pixels (Figure 6), they accounted for a small proportion of the total forest loss area (~9.5%).The number of image objects with object size ranging from 1 to 4 pixels accounted for 3.44%, 2.84%, 1.65%, 1.60% of the total forest loss area in the CONUS, respectively.
There was a total of 4,038,889 image objects that were greater than 5 pixels in size.At the regional level (see definition in Section 3.7 and Figure 4), the Southeast accounted for more than half of the total image objects greater than 5 pixels in size (50.22%) and had the largest average segment size (4.41 hectares) among the five regions.The Mountain West accounted for the third largest proportion of objects (up to 20.16%) with an average object size of 4.23 hectares.The Northeast ranked second in the number of objects (up to 21.86%), and objects averaged 2.7 hectares in size.

Effect of the Post-Classification of Fire Disturbance
Across all CONUS WELD tiles, about 5.1% fire areas initially classified by the random forest classifier were reclassified after a spatial and temporal check with the MODIS active fire data (Figure 7).For these reclassified forest loss areas, 73.6% were reclassified as stem removal and 26.4% were reclassified to stress.

Effect of the Post-Classification of Fire Disturbance
Across all CONUS WELD tiles, about 5.1% fire areas initially classified by the random forest classifier were reclassified after a spatial and temporal check with the MODIS active fire data (Figure 7).For these reclassified forest loss areas, 73.6% were reclassified as stem removal and 26.4% were reclassified to stress.

Effect of the Post-Classification of Fire Disturbance
Across all CONUS WELD tiles, about 5.1% fire areas initially classified by the random forest classifier were reclassified after a spatial and temporal check with the MODIS active fire data (Figure 7).For these reclassified forest loss areas, 73.6% were reclassified as stem removal and 26.4% were reclassified to stress.

Predictor Variable Ranking
Figure 8 presents the mean decrease in accuracy (relative variable importance) for all predictor variables input into the Random Forest classifier.Post-disturbance brightness (aBRI), minimum spectral trajectory slope for greenness (minGRE), post-disturbance NBR (aNBR), pre-disturbance redgreen index (bRGI), and pre-disturbance brightness (bBRI) were identified as the top five most important predictor variables for this classification.

Predictor Variable Ranking
Figure 8 presents the mean decrease in accuracy (relative variable importance) for all predictor variables input into the Random Forest classifier.Post-disturbance brightness (aBRI), minimum spectral trajectory slope for greenness (minGRE), post-disturbance NBR (aNBR), pre-disturbance red-green index (bRGI), and pre-disturbance brightness (bBRI) were identified as the top five most important predictor variables for this classification.

Classification of Forest Disturbances
Forest disturbance type classification for CONUS is shown in Figure 9.Of the total 172,686 km 2 of forest loss classified in this study, 83.75% was attributed to forest stem removal, 10.92% to fire, and 5.33% to stress.A cursory visual inspection shows good agreement between the classification results and independent datasets, such as MTBS and the insect mortality-area dataset.
The annual forest loss attributed by disturbance type is shown in Figure 10.In terms of highmagnitude disturbance (>50% forest cover loss) detected in the Global Forest Change Product, stem removal was the dominant forest disturbance agent at the CONUS scale over this time period.The average stem removal area per year was 16,069 km 2 , while forest loss due to fire and stress were much lower at 2096 km 2 and 1022 km 2 , respectively.The stem removal area shows a clear reduction in 2009-2011 period over the previous year, reflecting the reduced demand of timber due to the 2008 economic recession [66].Forest loss due to fire and stress was more variable from year to year, likely owing to strong climatic influence on these disturbance agents.

Classification of Forest Disturbances
Forest disturbance type classification for CONUS is shown in Figure 9.Of the total 172,686 km 2 of forest loss classified in this study, 83.75% was attributed to forest stem removal, 10.92% to fire, and 5.33% to stress.A cursory visual inspection shows good agreement between the classification results and independent datasets, such as MTBS and the insect mortality-area dataset.The annual forest loss attributed by disturbance type is shown in Figure 10.In terms of high-magnitude disturbance (>50% forest cover loss) detected in the Global Forest Change Product, stem removal was the dominant forest disturbance agent at the CONUS scale over this time period.The average stem removal area per year was 16,069 km 2 , while forest loss due to fire and stress were much lower at 2096 km 2 and 1022 km 2 , respectively.The stem removal area shows a clear reduction in 2009-2011 period over the previous year, reflecting the reduced demand of timber due to the 2008 economic recession [66].Forest loss due to fire and stress was more variable from year to year, likely owing to strong climatic influence on these disturbance agents.

Validation and Accuracy Assessment
Table 3 presents the Random Forest classification accuracy at the CONUS scale based on manual validation assessment of 4156 forest loss objects (488,791 pixels).Overall accuracy was 88.1%.While overall accuracy was high, classification performance varied by disturbance type.Stem removal had very high user's and producer's accuracies, demonstrating its relatively distinguishable feature space.Stress disturbance had the lowest user's accuracy (high commission error), mainly due to forest loss caused by fire being incorrectly attributed to stress disturbance.

Overall Performance
The lack of baseline forest disturbance data inhibits our ability to understand how climate change and human activity are altering forests and disturbance regimes [7,15,66].While many efforts have advanced our ability to detect and map forest disturbance [9,[21][22][23][24][25]39,67], the focus of the paper was to further our understanding of forest disturbance regimes by characterizing forest disturbance types, starting from the existing forest cover loss University of Maryland Global Forest Change Product [23].An object-based approach to characterize forest disturbance types across CONUS was prototyped and tested using WELD data.Through this methodology, we classified each forest cover loss object by disturbance type, resulting in a wall-to-wall map of forest disturbance types for CONUS from 2003 to 2011.
The methodology presented in this study advances our ability to characterize forest disturbance regimes at large spatial scales.Quantifying disturbance type using object-based methodology captures the spatial and temporal nature of forest disturbances better than pixel-based methodologies that suffer from noisy outputs due to isolated pixels [68].The use of Landsat WELD data enabled temporal trend analyses, which can provide valuable information regarding magnitude, duration, and direction of forest loss and recovery due to a particular disturbance agent.The high overall accuracy of the classification (88.1%), low training data needs (<1% of the total forest cover loss objects), and semi-automation suggest this methodology could be adapted to other regions where large-scale disturbances are common.It is important to note that, given the higher number of available Landsat scenes within CONUS compared to other regions [43], this methodology may not perform as well outside of CONUS.However, this methodology could potentially be adapted to other areas outside CONUS using satellite sensors with similar spectral and spatial characteristics, such as the European Space Agency's Sentinel sensors [69,70].Combined usage of Landsat and Sentinel data could improve disturbance classification by providing greater opportunities of suitable observations (i.e., higher frequency of cloud-free observations) [71].

Predictor Variable Importance
The top predictors in the disturbance type classification included tasseled cap brightness and greenness, NBR, RGI and NDVI.The performance of these indices is in part due to the large range of forest canopy cover change that can result from these disturbance types.For example, stem removal removes a large proportion of the forest canopy, leading to large magnitude, rapid spectral changes.Other disturbances, such as fire and insects, can remove variable amounts of forest cover depending on the intensity of the disturbance, and result in lower magnitude spectral change [72].Some of the top indices that utilize SWIR wavelengths (BRI, GRE, NBR) have had previous success quantifying different magnitudes of timber harvest (low intensity partial cut to clearcut) [73,74].The performance of indices that use SWIR in this study also matches other studies conducted in boreal forest, where SWIR indices outperformed those that did not use SWIR for delineating between fire and clearcuts [28].It is likely that the more heterogeneous post-fire landscapes, which consisted of mixes of live vegetation, dead standing vegetation, and shadow, contributed to this result.Other top indices that are sensitive to photosynthetically active vegetation (e.g., NDVI) have been shown to also be sensitive to fire intensity-induced stress in remaining live foliage [75], which could result in a lower magnitude (or less rapid) spectral change (compared to clearcuts) depending on the proportion of surviving trees within a forest disturbance object.Overall, it is important to understand that while some predictor variables are ranked higher than others, the classification accuracy is the result of the collective of input variables.

Limitations and Future Work
A limitation of the proposed approach is its grouping of forest disturbance types into only three classes (i.e., stress, fire, stem removal).In recent decades, forest loss in the eastern United States due to other disturbances (e.g., urban development, mining, reservoir construction, storm damage) has been estimated to account for ~43% of the total forest loss, with urban development and mining accounting for 34% and 7%, respectively [10].Forest harvest accounted for >56% of forest loss over the same time period [10].It is likely that these other disturbances were classified in the 'stem removal' class.Many other forest disturbance agents, such as hurricanes [76], wind [77] and ice storms [78] occur across CONUS.However, these disturbances account for a much smaller proportion of forest loss compared to the three disturbance classes considered in this study [8][9][10][11][12][13][14], and are mostly restricted to localized areas in the United States (e.g., hurricanes are most prevalent in areas along the Gulf of Mexico).Including these other disturbance agents would likely result in imbalanced training data and reduced accuracy in the overall classification [79].As a preliminary test to see whether the Global Forest Change Product detects storm damage, GoogleEarth imagery was used to identify forest damage along "major hurricane" (Category 3-6 hurricanes) paths mapped by the National Oceanic and Atmospheric Administration (https://coast.noaa.gov/hurricanes/).Of the major hurricanes that hit the southeastern U.S. during the time period of this study, only two had discernable forest damage; hurricanes Rita (2005) and Katrina (2005).Twenty areas where severe damage occurred (>75% of trees knocked over) were identified in southeast Texas (WELD tile h18v17) and southern Louisiana and Mississippi (WELD tile h21v16).Only ~15% of storm-damaged areas identified in GoogleEarth were captured by the Global Forest Change Product, highlighting potential detection limitations of the product.Future work is needed to assess what disturbances and disturbance intensities are not detected by this product and how classification performance is affected by including more disturbance types.
As illustrated by the preliminary storm damage assessment, this paper's classification approach is limited by the Global Forest Change Product detection capabilities and accuracy.The Global Forest Change Product was designed to capture high-intensity disturbances (>50% forest canopy loss) and thus can miss substantial forest loss resulting from lower-intensity disturbances (<50% forest canopy loss) [13].It is clear that variability in disturbance magnitude (i.e., proportion of forest canopy removed) can result in similar spectral changes and misclassification of disturbance type.For example, stress disturbance had the lowest user's accuracy (high commission error), mainly due to forest loss caused by fire being incorrectly attributed to stress disturbance.It is common for fire disturbance in forests to result in highly variable within-patch severity (i.e., proportion of forest overstory mortality) [80,81], which could account for this error.However, because the MODIS Active Fire product was used in the post-classification phase, it is likely that some small fires were not detected in the coarser resolution MODIS product (1 km spatial resolution) and, therefore, were classified erroneously as stress disturbance rather than fire.Additionally, forest loss resulting from low intensity, multi-year partial cut harvests (few trees harvested per unit area) could be misclassified as stress-induced mortality due to similar changes in SWIR and NIR reflectance [73].It is also quite possible that low intensity partial cut areas were not captured in the Global Forest Change Product that was used in this study, resulting in an underestimation of forest loss.Partial cuts can result in low overall spectral change and variability, leading to misclassification as undisturbed forest cover [74].Likewise, areas of undisturbed forest misclassified as forest loss by the Global Forest Change Product (false positives) could have been erroneously classified into one of the three classes.
While the approach detailed in this study provides a relatively accurate broad-scale approach using the Global Forest Change Product, future work is needed to further refine this methodology and understand its limitations.Our methodology is limited to the temporal range that the Global Forest Change Product is produced (i.e., 2000 to the present) but could easily be adapted to additional years as the product is updated and expanded.The relatively short time series used in this study did not allow for the permanence of forest loss to be assessed for many of the disturbance objects.For example, many forest loss objects that returned to vegetation had spectral index values that were indistinguishable from objects that did not return to vegetation for the first 1-3 years after disturbance.As the Global Forest Change Product expands, longer time series datasets (that include >3 years of post-disturbance observations) could potentially be used to assess the permanence of forest cover

Figure 1 .
Figure 1.Temporal segmentation of BRI, GRE, NBR, and B5 spectral indices for each of the three forest disturbance types.Typical spectral profiles for forest stem removal (top row), fire (middle row), and stress (bottom row) are shown.Vertical lines denote the year when disturbance was detected by the Global Forest Change Product.

Figure 1 .
Figure 1.Temporal segmentation of BRI, GRE, NBR, and B5 spectral indices for each of the three forest disturbance types.Typical spectral profiles for forest stem removal (top row), fire (middle row), and stress (bottom row) are shown.Vertical lines denote the year when disturbance was detected by the Global Forest Change Product.

Figure 2 .
Figure 2. Illustration of the extracted features for every spectral index using the forest loss objects and WELD time series data.Each temporal profile is converted into segments and five features, used in the subsequent classification, are extracted: (1) maximum (max) and minimum (min) slope among all the segments of the trajectory; (3) mean value of the two years before the forest loss year (b), (4) mean value of the two years after the forest loss year (including the disturbance year) (a); and (5) difference between (4) and (3) (delta).

Figure 2 .
Figure 2. Illustration of the extracted features for every spectral index using the forest loss objects and WELD time series data.Each temporal profile is converted into segments and five features, used in the subsequent classification, are extracted: (1) maximum (max) and (2) minimum (min) slope among all the segments of the trajectory; (3) mean value of the two years before the forest loss year (b), (4) mean value of the two years after the forest loss year (including the disturbance year) (a); and (5) difference between (4) and (3) (delta).

Figure 3 .
Figure 3. WELD tiles selected for training samples: green, red, blue, and pink filled tiles represent selected training tiles for stem removal, fire, stress, and both stem removal and fire, respectively.To provide geographic context, ecoregions grouped by common forest type are shown in background as well as the total forest cover extent mapped by Hansen et al. [23] (dark green pixels).

Figure 3 .
Figure 3. WELD tiles selected for training samples: green, red, blue, and pink filled tiles represent selected training tiles for stem removal, fire, stress, and both stem removal and fire, respectively.To provide geographic context, ecoregions grouped by common forest type are shown in background as well as the total forest cover extent mapped by Hansen et al. [23] (dark green pixels).

Figure 4 .
Figure 4. Regional stratification of WELD tiles composing the validation dataset.WELD tiles selected for validation are shown with thicker border lines.Ecoregions grouped by common forest type are shown in background as well as the total forest cover extent mapped by Hansen et al. [23] (dark green pixels).

Figure 4 .
Figure 4. Regional stratification of WELD tiles composing the validation dataset.WELD tiles selected for validation are shown with thicker border lines.Ecoregions grouped by common forest type are shown in background as well as the total forest cover extent mapped by Hansen et al. [23] (dark green pixels).

Figure 5 .
Figure 5. Representative forest loss object size, shape, and pattern for fire and stem removal in northern California (left), stem removal in northwest Oregon (middle), and stress in central Colorado (right).Forest loss objects are displayed with random colors.

Figure 6 .
Figure 6.Forest disturbance object size distribution by region: (a) Mountain West, (b) Lowland West, (c) Central, (d) Northeast, and (e) Southeast.Note log-scale x and y axes.

Figure 5 . 24 Figure 5 .
Figure 5. Representative forest loss object size, shape, and pattern for fire and stem removal in northern California (left), stem removal in northwest Oregon (middle), and stress in central Colorado (right).Forest loss objects are displayed with random colors.

Figure 6 .
Figure 6.Forest disturbance object size distribution by region: (a) Mountain West, (b) Lowland West, (c) Central, (d) Northeast, and (e) Southeast.Note log-scale x and y axes.

Figure 6 .
Figure 6.Forest disturbance object size distribution by region: (a) Mountain West, (b) Lowland West, (c) Central, (d) Northeast, and (e) Southeast.Note log-scale x and y axes.

Figure 7 .
Figure 7. Initial forest disturbance type classification (Left) and the final result checked with the Moderate Resolution Imaging Spectroradiometer (MODIS) active fire product (right) for northwest Oregon (top row) and northwest California (bottom row).Orange, magenta, and blue objects represent stem removal, fire and stress disturbance classes, respectively.Forest cover extent is shown in light grey, non-forest as black, and water as white in background.Each panel is approximately 25 × 25 km.

Figure 7 .
Figure 7. Initial forest disturbance type classification (Left) and the final result checked with the Moderate Resolution Imaging Spectroradiometer (MODIS) active fire product (right) for northwest Oregon (top row) and northwest California (bottom row).Orange, magenta, and blue objects represent stem removal, fire and stress disturbance classes, respectively.Forest cover extent is shown in light grey, non-forest as black, and water as white in background.Each panel is approximately 25 × 25 km.

Figure 8 .
Figure 8. Predictor variable importance based on the mean decrease in accuracy for the Random Forest classification.

Figure 8 .
Figure 8. Predictor variable importance based on the mean decrease in accuracy for the Random Forest classification.

Figure 9 .
Figure 9. Forest loss (2003-2011) classified by major disturbance agent across CONUS.Insets highlight spatial extent and pattern of disturbance agents in the west and southeast United States.Total forest cover extent is shown in light grey as background.

Figure 9 .
Figure 9. Forest loss (2003-2011) classified by major disturbance agent across CONUS.Insets highlight spatial extent and pattern of disturbance agents in the west and southeast United States.Total forest cover extent is shown in light grey as background.

Figure 9 .
Figure 9. Forest loss (2003-2011) classified by major disturbance agent across CONUS.Insets highlight spatial extent and pattern of disturbance agents in the west and southeast United States.Total forest cover extent is shown in light grey as background.

Figure 10 .
Figure 10.Total forest loss (km 2 ) in CONUS attributed to stress, fire, and stem removal from 2003-2011.Data for objects of all sizes is shown.

Figure 10 .
Figure 10.Total forest loss (km 2 ) in CONUS attributed to stress, fire, and stem removal from 2003-2011.Data for objects of all sizes is shown.

Table 3 .
Classification confusion matrix results at the CONUS scale, reporting the number of objects.Correctly classified forest disturbance pixels are highlighted in grey.No disturbance (false positives) could be attributed to 2.5% of the objects.