Shifts in Forest Structure in Northwest Montana from 1972 to 2015 Using the Landsat Archive from Multispectral Scanner to Operational Land Imager

There is a pressing need to map changes in forest structure from the earliest time period possible given forest management policies and accelerated disturbances from climate change. The availability of Landsat data from over four decades helps researchers study an ecologically meaningful length of time. Forest structure is most often mapped utilizing lidar data, however these data are prohibitively expensive and cover a narrow temporal window relative to the Landsat archive. Here we describe a technique to use the entire length of the Landsat archive from Multispectral Scanner to Operational Land Imager (M2O) to produce three novel outcomes: (1) we used the M2O dataset and standard change vector analysis methods to classify annual forest structure in northwestern Montana from 1972 to 2015, (2) we improved the accuracy of each yearly forest structure classification by applying temporal continuity rules to the whole time series, with final accuracies ranging from 97% to 68% respectively for two and six-category classifications, and (3) we demonstrated the importance of pre-1984 Landsat data for long-term change studies. As the Landsat program continues to acquire Earth imagery into the foreseeable future, time series analyses that aid in classifying forest structure accurately will be key to the success of any land management changes in the future.


Introduction
Land managers, foresters, wildlife biologists, and other environmental scientists all need inexpensive and easily accessible information characterizing current and past conditions in order to monitor changes in the age/size of forests under their purview [1][2][3][4][5].Time series analysis tools such as harmonic regression [6], LandTrendr [7] and TimeSync [8], LandsatLinkr [9], season-integrated products [10], and whole time series analyses [11] have great utility for mapping historic and current forest structure.Managers can observe how forests have increased or decreased across structure classes within a landscape context.This can be especially important in multiple-use areas where different uses are often opposed to one another, for instance, on US Forest Service land where tree harvest and including from natural events such as fire and insect infestations and from human influences such as tree harvests.The majority of the study area is covered by six Landsat WRS-1 scenes (MSS) and four WRS-2 scenes (TM, ETM+, OLI) (Figure 1).

Forest Structure Reference Data
We visually classified forest structure at 674 field data sites to be used as reference data for training the initial classification.We used a segmented polygon layer characterizing forest stand that was developed using eCognition [26] to guide us to these locations so that we identified, at minimum, 100 training points per class.Three analysists independently interpreted each of the reference polygons using National Agriculture Imagery Program (NAIP) imagery collected in Montana during 2013 (when the field data sites were collected) along with imagery available in GoogleEarth.The interpreted results from each analyst were compared and the final class for each polygon was chosen as the class identified by two or more of the analysts.Pixel values were extracted from each of these polygons to serve as training data.We defined six structure classes as follows: Alpine (A), Open (O), Stand Initiation (S), Thin (T), Advanced Regeneration (R), and Mature (M).The categories were selected because they captured the gradient in distinct forest structural stages within our study area along with an understanding of wildlife-forest structure relationships in a parallel study [27].

Forest Structure Reference Data
We visually classified forest structure at 674 field data sites to be used as reference data for training the initial classification.We used a segmented polygon layer characterizing forest stand that was developed using eCognition [26] to guide us to these locations so that we identified, at minimum, 100 training points per class.Three analysists independently interpreted each of the reference polygons using National Agriculture Imagery Program (NAIP) imagery collected in Montana during 2013 (when the field data sites were collected) along with imagery available in GoogleEarth.The interpreted results from each analyst were compared and the final class for each polygon was chosen as the class identified by two or more of the analysts.Pixel values were extracted from each of these polygons to serve as training data.We defined six structure classes as follows: Alpine (A), Open (O), Stand Initiation (S), Thin (T), Advanced Regeneration (R), and Mature (M).The categories were selected because they captured the gradient in distinct forest structural stages within our study area along with an understanding of wildlife-forest structure relationships in a parallel study [27].Detailed definitions of these classes can be seen in Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes.
Using imagery available in GoogleEarth that covered the entire study area, we identified 50 random points per reference class for each of the years 2013, 2005, and 1995 to be used for independent validations of the predicted classification maps.In order to perform an accuracy assessment on an early classification, we purchased 1:16,000-to 1:40,000-scale historical aerial photographs from the US Department of Agriculture Aerial Photography Field Office that covered only the northwestern portion of the study area (no historic photos were available for the southern and northeastern portions) and identified a minimum of 50 random points per reference class for the 1975 classification.

Ancillary Data
Several ancillary data sets were acquired for analysis and validation.We acquired a 30-m digital elevation model (DEM) of the study area from the USGS National Elevation Dataset (ned.usgs.gov).We downloaded hydrology data from the USGS National Hydrology Dataset (NHD, nhd.usgs.gov)and used them to mask water bodies from the study area.Finally, we acquired four-band, one-meter-pixel NAIP images from 2013 from the Montana State Library Natural Resource Information System (nris.mt.gov).

Landsat Data Processing through LandsatLinkr
We processed the imagery using R statistical software and the LandsatLinkr (LLR v.0.2.0 (beta)) package [9].The image data that LLR produces are annual cloud-free composites for Tasseled Cap spectral indices [28,29] that are comparable across sensors.Within LLR, TM, and ETM+, surface reflectance data are considered the spatial and spectral standard to which MSS and OLI data are modeled using multiple linear regression based on the relationship between sets of coincident and near-date image pairs between MSS and TM, and ETM+ and OLI.Prior to inter-sensor harmonization, LLR resamples Landsat Pre-Collection Level-1 MSS data to 30-m (nearest neighbor), improves georegistration [30], produces cloud masks [31], and calculates surface reflectance with the COST method [32], where the Lhaze parameter or "dark object value" is estimated from an automated implementation of the histogram method [33].TM, ETM+ (Landsat Collection 1 Level-2), and OLI (Landsat 8 OLI/TIRS C1 Level-2) data required less pre-processing since they were delivered from USGS as surface reflectance (LEDAPS algorithm [34] and L8SR algorithm [35]) and included cloud masks (Fmask algorithm [36]).In the final step of annual compositing, for each year in the time series, all images that spatially intersected our study area were merged using the overlapping per-pixel mean while filling in where pixels were identified as cloud or cloud shadow in the masks.Where there were no unmasked pixels, a null value was assigned.In our case, there were several years that contained large regions of null data because of extensive, consistent cloud cover.These data were produced for each year from 1972 to 2015.
Unpublished data generated during LLR development found that Pearson's correlation (r) between predicted (MSS and OLI) and observed (TM and ETM+) annual tasseled cap index values ranged from 0.82 to 0.96 for a test scene (WRS-2 path 38 row 29-an area very similar in cover types to those included in this study).MSS tasseled cap wetness was about 10% lower in correlation than brightness and greenness, likely because MSS does not include a SWIR band, which is weighed heavily in tasseled cap wetness.OLI correlation to ETM+ was less than 0.9 for all three tasseled cap indices, likely because predictive models between ETM+ and OLI are based on near-date images, unlike MSS and TM, which use coincident images.

Ancillary Data Processing
We derived ancillary datasets from the DEM and NAIP imagery, including slope, aspect, and texture information for the study area.The image texture mean and minimum values were derived from the standard deviation of the first principal component of the 1-m NAIP data, then resampled to 30-m pixels for analysis [26].

Classification of Master Image
We executed a variety of forest structure classifications using several different classification methods (from the Caret package in R (version 3.4.1):http://topepo.github.io/caret/available-models.html).We then chose the classification methods with the highest relative overall accuracy to create each of the initial 6-class to 2-class predictions for 2013 (see Section 2.5.1).

Accuracy Assessment of the Initial 2013 Classification
We randomly selected 1000 pixels, began determining the reference class of each sample pixel, and continued this selection protocol sequentially until exactly 50 pixels were selected in each reference class (except for the Open class-very few random points fell within Open areas, so we used only 36 points for this class in the initial 2013 classifications).Additional pixels (of the 1000) were discarded once the class had reached its quota of 50.We did not conduct an analysis of "from-to" change accuracy because of the absence of reference data for such an analysis.Instead, we followed a newer approach to minimize error propagation in analysis of "from-to" change [37] in terms of its summary of best practices for validation.A weighted error matrix [37] was then calculated for five levels of the initial 2013 classification (Figure 2), however, we estimated the proportion of area for each class based on the reference classification rather than the map classification since our validation data were stratified by reference class.User's and producer's accuracies along with the standard error were calculated for each level of classification as well [37].
Forests 2018, 9, x FOR PEER REVIEW 5 of 21 from the standard deviation of the first principal component of the 1-m NAIP data, then resampled to 30-m pixels for analysis [26].

Classification of Master Image
We executed a variety of forest structure classifications using several different classification methods (from the Caret package in R (version 3.4.1):http://topepo.github.io/caret/availablemodels.html).We then chose the classification methods with the highest relative overall accuracy to create each of the initial 6-class to 2-class predictions for 2013 (see Section 2.5.1).

Accuracy Assessment of the Initial 2013 Classification
We randomly selected 1000 pixels, began determining the reference class of each sample pixel, and continued this selection protocol sequentially until exactly 50 pixels were selected in each reference class (except for the Open class-very few random points fell within Open areas, so we used only 36 points for this class in the initial 2013 classifications).Additional pixels (of the 1000) were discarded once the class had reached its quota of 50.We did not conduct an analysis of "fromto" change accuracy because of the absence of reference data for such an analysis.Instead, we followed a newer approach to minimize error propagation in analysis of "from-to" change [37] in terms of its summary of best practices for validation.A weighted error matrix [37] was then calculated for five levels of the initial 2013 classification (Figure 2), however, we estimated the proportion of area for each class based on the reference classification rather than the map classification since our validation data were stratified by reference class.User's and producer's accuracies along with the standard error were calculated for each level of classification as well [37].

Change Vector Analysis
Open and Stand Initiation were initially difficult to classify separately with sufficient accuracy, so we performed the change analysis using the 2013 5-class prediction as the master map (with O and S combined as one class).We utilized a change vector analysis (CVA) process [38] to classify each year of imagery compared to the master 2013 map.Advantages of using CVA for change analysis are that the process reduces compound error by not reclassifying unchanged locations and it provides reliable training data for years when data are unavailable or difficult to obtain.Compound error in locations that have potentially changed is a concern with this process; however, we can avoid high rates of compound error by utilizing CVA combined with temporal continuity rules applied to the resulting time series (see Section 2.7).CVA was applied to the three tasseled cap components, brightness, greenness, and wetness, to produce a change map for each year as follows:

Change Vector Analysis
Open and Stand Initiation were initially difficult to classify separately with sufficient accuracy, so we performed the change analysis using the 2013 5-class prediction as the master map (with O and S combined as one class).We utilized a change vector analysis (CVA) process [38] to classify each year of imagery compared to the master 2013 map.Advantages of using CVA for change analysis are that the process reduces compound error by not reclassifying unchanged locations and it provides reliable training data for years when data are unavailable or difficult to obtain.Compound error in locations that have potentially changed is a concern with this process; however, we can avoid high rates of compound error by utilizing CVA combined with temporal continuity rules applied to the resulting time series (see Section 2.7).CVA was applied to the three tasseled cap components, brightness, greenness, and wetness, to produce a change map for each year as follows: where 1,2 refer to 2013 (the master image) and each other yearly image respectively, and B = brightness, G = greenness, and W = wetness.A threshold of possible change was assigned to each yearly change map that attempted to insure that all possible changes were classified during the process (on average, approximately 44% of the pixels in each change image were identified as possibly changed).Using C5.0 in R (as it was consistently the best method during the agnostic classification processes (Section 2.5)) and the Tasseled Cap components, elevation, slope, and aspect as predictor variables, only those pixels that fell within the threshold of possible change were reclassified and a random selection of 10,000 pixels (2000 per class) that were within the threshold of unlikely to have changed were used as the training data.

Informing Yearly Classifications with Time Series
We utilized all 44 classified maps to help improve the accuracy of each yearly classification by (1) filling data gaps and (2) using knowledge-based rules to correct for temporal inconsistencies.This process also enabled us to separate Open from Stand Initiation for all of the yearly classifications.This was accomplished through applying several consecutive temporal continuity rules-based on local, ecosystem, and successional knowledge, and expected change behaviors-to the entire time series as follows (code/logic used for each rule can be seen in Supplementary Material: Logic S1: Logic Used for Time-Series-Informed Rules; examples of rules 1 through 6 are shown in Figure 3): 1.
Gap fill #1 (where no data existed in a year-due to clouds, scan line gaps, and/or missing scenes for that year): utilized data from surrounding years by filling with the class of the previous year if it was not Unknown, or with the following year if the class of the previous year was Unknown; remained Unknown if the classes of the surrounding years were Unknown; 2.
Continuity fill: looked for consistent Alpine, Mature, and Open/Stand Initiation (i.e., 1972 = A, 2015 = A, and the majority of the time series at this pixel = A, then fill this pixel with A); 3.
Smoothing: looked for single occurrences within the time series (i.Alpine by elevation (any Open pixel ≥ 2500 m was labeled as Alpine); 8.
Mask low elevation to exclude areas of non-forest (exclude all pixels <500 m).
Resulting in yearly maps from 1972 to 2015, each classified into six forest structure classes as defined in Section 2.2.2.
We calculated trajectories of each of the six structure classes over the 44-year time period using the time-series-informed classifications.The trajectories displayed the percentage of the study area covered by each structure class.These trajectories allowed us to observe trends in regeneration, growth, and loss of forest cover over time within the study area.7) Alpine based on elevation, and (8) Non-forest based on elevation (each rule is described in detail in Section 2.7 and code/logic used for each rule can be seen in Supplementary Material: Logic S1: Logic Used for Time-Series-Informed Rules).

Accuracy Assessment for Initial and Time-Series-Informed Classifications
Due to lack of adequate imagery for each year, we performed independent accuracy assessments on four years spread across the four decades (2013, 2005, 1995, and 1975).Validation data for both the initial and time-series-informed classifications were selected for these years in the same manner as described in Section 2.5.1, with the exception that slightly more than 50 points per class were identified for the 1975 classifications.A weighted error matrix (see Section 2.5.1 [37]) was created for the 6-class predictions to demonstrate the overall accuracies at the most detailed level.The error matrices for the 5-class through 2-class predictions were created by combining classes from the 6class error matrices.Eight combinations of classes were tested (Figure 4).We compared the accuracies of the initial classifications to the time-series-informed classifications to evaluate the influence of the continuity rules on the time series.

Accuracy Assessment for Initial and Time-Series-Informed Classifications
Due to lack of adequate imagery for each year, we performed independent accuracy assessments on four years spread across the four decades (2013, 2005, 1995, and 1975).Validation data for both the initial and time-series-informed classifications were selected for these years in the same manner as described in Section 2.5.1, with the exception that slightly more than 50 points per class were identified for the 1975 classifications.A weighted error matrix (see Section 2.5.1 [37]) was created for the 6-class predictions to demonstrate the overall accuracies at the most detailed level.The error matrices for the 5-class through 2-class predictions were created by combining classes from the 6class error matrices.Eight combinations of classes were tested (Figure 4).We compared the accuracies of the initial classifications to the time-series-informed classifications to evaluate the influence of the continuity rules on the time series.

Comparison to US Forest Service Forest Inventory and Analysis Data
We independently evaluated our final forest structure classification for 2013 using field data collected as part of the US Forest Service Forest Inventory and Analysis (FIA) program during 2006-2012.We spatially assigned each subplot (~170 m 2 from FIA) the forest structure class that was mapped for that pixel during 2013.We assessed differences in FIA-derived forest metrics across forest structure classes by implementing a multivariate analysis of variance (MANOVA) visually displaying boxplots (i.e., median and the interquartile range; IQR).Our structural metrics included basal area weighted diameter at breast height (DBH), derived canopy cover and tree height using the Forest Vegetation Simulator [39], as well as stem density for trees larger than 12.7 cm (i.e., 5 inches) in DBH.Basal area weighted DBH is calculated with Equation (2).
We expected to observe structural differences in our forest structure classes consistent with the aforementioned forest structure class definitions (Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes.).

Classification of Master Image
The best 6-class prediction for the 2013 master image had an overall accuracy of 68% and a standard error of 3%, while the 2-class prediction had an overall accuracy of 98% with a standard error of <1% (Table 1).We found that structure classes Open and Stand Initiation were particularly difficult to distinguish as evidenced by the producer's and user's accuracies (Table 2, 6-Class) so we combined these classes to improve accuracy (Table 2, 5-Class).Structure classes of Thin, Advanced Regeneration, and Mature had some error as well (Table 2).Complete error matrices for all six initial 2013 classifications in Table 1 can be seen in Supplementary Material, Tables S2-S7.

Yearly Forest Structure Classifications
We identified possible change pixels for 43 years of imagery using the 2013 image as the master in 43 CVA processes.On average, approximately 44% of the pixels in each change image were identified as possibly changed (these numbers include large data gaps in the original imagery where, for a particular year, one or more scenes were cloud-covered and thus not included in the analysis).The threshold of change values over-predicted possible change, but we chose to insure that all possible change pixels would be reclassified during the CVA process.Thus, we ended up reclassifying many possibly unchanged pixels, however, nearly all unchanged pixels were likely reclassified to their original 2013 class.
We produced weighted error matrices for the initial classifications of 2005, 1995, and 1975 (complete error matrices can be seen in Supplementary Material, Tables S8-S25).The overall accuracy for each of the 5-class classifications (where Open and Stand Initiation were combined as one class) ranged from 55% for 2005, 58% for 1975, and 60% for 1995 (Table 1) (recall that the initial master 2013 classification had an overall accuracy of 73%).Overall accuracies for each of the 2-class classifications ranged from 97% for 1995 and 2005, and 99% for 1975 (Table 1).

Time-Series-Informed Yearly Classifications
After using the time series and temporal continuity rules to inform each of the yearly classifications, we produced weighted error matrices for every level of classification (eight levels shown in Figure 4) for the final 2013, 2005, 1995, and 1975 predictions (complete error matrices can be seen in Supplementary Material, Tables S26-S57).We observed increased overall accuracies and decreased or equivalent standard errors for most of the classifications (Table 3) and those that did not increase in overall accuracy decreased by less than 1%.Most importantly, the overall accuracy of the 6-class prediction for 2013 increased by 4% and the overall accuracies for the 6   The temporal continuity rules applied to the time series enabled Open to be separated from Stand Initiation in the final 2013 classification resulting in better all-around user's and producer's accuracies for the two classes in the 6-class prediction (Table 4).Thin and Advanced Regeneration are still somewhat confused with one another within the classifications (Table 4).However, the user's and producer's accuracies for the Thin, Advanced Regeneration, and Mature classes in the 6-and 5-class either remained the same or increased from the initial classifications (Tables 2 and 4).
Our assessment of field-derived forest structural metrics across our forest structure classes included a sample of 500 FIA subplots distributed across our 2013 forest structure classes: Alpine (n = 14), Open (n = 71), Stand Initiation (n = 43), Thin (n = 79), Advanced Regeneration (n = 62), and Mature (n = 231).Our MANOVA indicated statistical differences among forest structure classes (F 20,1928 = 11.33,p <0.001) and the patterns were intuitive and consistent (Figure 5).For instance, tree size, tree density, and tree heights were lower in our Alpine and Stand Initiation classes relative to all other classes (Figure 5).This was followed by our Open class, which exhibited slightly larger tree sizes, higher tree densities, and taller trees.Our Thin class exhibited larger tree sizes, higher tree densities, taller trees, and more canopy cover than our Open, Stand Initiation, and Alpine classes.Further, our Advanced Regeneration and Mature classes exhibited more trees, higher canopy cover, and taller trees than all other classes, and Mature was generally greater than Advanced Regeneration for all metrics.

Time Series Analysis
The overall accuracies for the initial master 2013 classification ranged from 68% for the 6-class map to 98% for the 2-class map (Table 1).We observed an increase in the overall accuracy of the final 2013 6-class classification of 4% after applying the temporal continuity rules (Table 3).The final overall accuracies for all years of classifications that separated the Alpine class and separated the more closed-canopy forested classes (Advanced Regeneration and Mature) from sparse or unforested classes (Open, Stand Initiation, and Thin) were all above 85% (Table 3, 2-and 3-class classifications).The clear variation of accuracies from year to year is in part because change vector analysis is likely to have greater errors as time from the master year increases, resulting in larger numbers of changed pixels and therefore increased potential for compound error.Accuracies decreased markedly as we attempted classifications that separated the different closed-or near-closed-canopy forests from each other, especially for the older imagery.The 6-class classifications that separated classes without trees (Open vs. Stand Initiation) resulted in the lowest overall accuracies, and this was especially true with the MSS data.Open and Stand Initiation are characterized by similar vegetation amounts, with the primary differences attributable to species, i.e., Stand Initiation has seedling conifer species, where Open is composed of more grasses and/or shrubs.This type of distinction within a classification could benefit from the greater spectral and radiometric qualities of the later satellite systems.Also, the difficulty distinguishing Open from Stand Initiation might be a function of using only Tasseled Cap components for the classifications and change analyses and not incorporating the thermal data acquired by Landsat sensors into these analyses-the thermal data could have helped identify areas affected by stand-replacing-disturbance [40].
The trajectories that were created with these time series provide critical information to land managers and allow for trend analyses changes that occur over long periods of time (Figure 6).We observed trends in each of the six forest structure classes and compared them to historic knowledge and expected overall changes in the study area after disturbances.We expected to see forest regeneration after harvest or fire [41][42][43], for example, and these trajectories show a general decline in the amount of Mature trees over time along with a general increase in Stand Initiation and Advanced Regeneration trees (Figure 6), supporting the expectation from a pattern of forest harvesting and especially wildfires over the last decade within this study area.The Alpine class remains mostly the same over time, with some decrease shown from 1972 to 1984 (from the MSS data) (Figure 6).This slight decrease could be a function of increases in tree line [44][45][46].Open also remains mostly the same with a larger decrease shown from 1972 to 1984 (Figure 6).The portions of these trajectories before 1984 (i.e., from MSS data) are all more irregular/noisier than the portions from 1984 forward (Figure 6).We believe the forthcoming enhancements and calibration of MSS data [47] will help to smooth the pre-1984 data and will subsequently allow the early classification trajectories to correspond more readily to the later trajectories.The higher variability in early years, however, might also be partially the result of changing forest management practices following the adoption of the National Forest Management Act in 1976.
In addition to overall trends, we can analyze portions of each class trajectory to observe large changes as they happen.For instance, from 2011 to 2013, we see an increase in the Mature class while at the same time we see a similarly rapid decrease in the Advanced Regeneration class and a smaller decrease in the Thin class (Figure 6).This could be a function of the successional processes of the forest where over a long period of time Advanced Regeneration develops into Mature, thus increasing the Mature class and decreasing the Advanced Regeneration class.On the other hand, the Stand Initiation class shows a small increase from 2011 to 2012, corresponding to the decrease of the Advanced Regeneration class at the same time.This could be explained by a large number of fires and harvests in 2011 and 2012 resulting in the loss of some forest cover that is replaced by the Stand Initiation class.While recorded drastic changes and the physical relationships among the forest structure classes can be logical explanations for the changes observed in the trajectories, we must also consider that there might be artifacts in the data that produce faulty information.For example, the final 2015 classification contains some errors (described below in Section 4.2) that are reflected in the last segment of these class trajectories.Changes from one class to another can be observed with a Sankey diagram (Figure 7; decadal diagram in Supplementary Material, Figure S1).The Sankey diagram illustrates the state transition from one class to another over our 44-year study time period.The diagram shows nodes that represent classes and links between those nodes whose heights are proportional to the area of each class at a given period of time and the amount transitioned to the next period of time.Similar to the general decline of the Mature class trajectory (Figure 6), much Mature forest present in 1972 has been transformed nearly evenly to the three other forest types, with just slightly more transition to Stand Initiation, followed by Thin, and finally Advanced Regeneration (Figure 7).This could be due to a combination of clearcut harvests followed by succession and thinning either by harvest or disturbance events.Approximately one quarter of the area of the Open class present in 1972 has transformed to forested classes, spread somewhat evenly between Thin, Advanced Regeneration, and Mature.This could be the result of new forest encroachment as existing stand edges, or possible misclassification in 1972 of the Stand Initialization class, which over-predicted the amount of Open at that time.In general, we observed that over the 44-year study time period the region has become more forested and the forests are younger, less dense, and more uniform (based on our class definitions in Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes) (Figure 7).Changes from one class to another can be observed with a Sankey diagram (Figure 7; decadal diagram in Supplementary Material, Figure S1).The Sankey diagram illustrates the state transition from one class to another over our 44-year study time period.The diagram shows nodes that represent classes and links between those nodes whose heights are proportional to the area of each class at a given period of time and the amount transitioned to the next period of time.Similar to the general decline of the Mature class trajectory (Figure 6), much Mature forest present in 1972 has been transformed nearly evenly to the three other forest types, with just slightly more transition to Stand Initiation, followed by Thin, and finally Advanced Regeneration (Figure 7).This could be due to a combination of clearcut harvests followed by succession and thinning either by harvest or disturbance events.Approximately one quarter of the area of the Open class present in 1972 has transformed to forested classes, spread somewhat evenly between Thin, Advanced Regeneration, and Mature.This could be the result of new forest encroachment as existing stand edges, or possible misclassification in 1972 of the Stand Initialization class, which over-predicted the amount of Open at that time.In general, we observed that over the 44-year study time period the region has become more forested and the forests are younger, less dense, and more uniform (based on our class definitions in Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes) (Figure 7).

Classification of Forest Structure
Forest structure maps often provide a picture of a forest for only one instance in time.From a management perspective, being able to identify forest structure over time using image analysis is of great importance.There is also a pressing need to spatially map changes in forest structure given forest management policies and accelerated disturbances from climate change.
This collection of 44 years of classified maps allows land managers to visualize changes to forest structure that have happened over decades (Figure 8).Users can observe successional trends in areas that have experienced gradual changes such as beetle infestations or drastic, stand-replacing disturbances such as harvest, or fire.Fire scars (red circles in Figure 8, 2007) and clear cuts (blue circle around checkerboard pattern in Figure 8, 2012) are the most obvious visual changes.Insect damage that affects trees slowly is more difficult (but not impossible) to observe as a gradual decrease of the Mature class over the course of the time series (as seen in the Mature trajectory, Figures 6 and 7 and from 1972 to 2015 in black circles in Figure 8).Our maps show clear evidence of changes across the landscape over time.
We were unable to perform accuracy assessments on all 44 classified images (initial and final), but we made assumptions that the accuracies of all the final maps in the time series are similar to those of 2013, 2005, 1995, and 1975, especially after applying the temporal continuity rules.Of note is that several of the continuity rules rely on surrounding dates.The first and last dates in the time series, therefore, might be less reliable, because they cannot take advantage of these particular rules.An example can be seen in the 2015 results where Landsat 7 imagery was used (Figure 8).The effects of the scan line corrector failure remain evident, where in other years the continuity rules might correct many of these errors.

Classification of Forest Structure
Forest structure maps often provide a picture of a forest for only one instance in time.From a management perspective, being able to identify forest structure over time using image analysis is of great importance.There is also a pressing need to spatially map changes in forest structure given forest management policies and accelerated disturbances from climate change.
This collection of 44 years of classified maps allows land managers to visualize changes to forest structure that have happened over decades (Figure 8).Users can observe successional trends in areas that have experienced gradual changes such as beetle infestations or drastic, stand-replacing disturbances such as harvest, or fire.Fire scars (red circles in Figure 8, 2007) and clear cuts (blue circle around checkerboard pattern in Figure 8, 2012) are the most obvious visual changes.Insect damage that affects trees slowly is more difficult (but not impossible) to observe as a gradual decrease of the Mature class over the course of the time series (as seen in the Mature trajectory, Figures 6 and 7 and from 1972 to 2015 in black circles in Figure 8).Our maps show clear evidence of changes across the landscape over time.
We were unable to perform accuracy assessments on all 44 classified images (initial and final), but we made assumptions that the accuracies of all the final maps in the time series are similar to those of 2013, 2005, 1995, and 1975, especially after applying the temporal continuity rules.Of note is that several of the continuity rules rely on surrounding dates.The first and last dates in the time series, therefore, might be less reliable, because they cannot take advantage of these particular rules.An example can be seen in the 2015 results where Landsat 7 imagery was used (Figure 8).The effects of the scan line corrector failure remain evident, where in other years the continuity rules might correct many of these errors.We demonstrated that remote sensing data that are freely available were sufficient to identify forest structure classes across diverse forested landscapes.The fact that field-derived FIA forest structural metrics data (Figure 5) followed expected patterns of forest characteristics across structure classes highlights the achievement that structure classes were ecologically meaningful.These strong We demonstrated that remote sensing data that are freely available were sufficient to identify forest structure classes across diverse forested landscapes.The fact that field-derived FIA forest structural metrics data (Figure 5) followed expected patterns of forest characteristics across structure classes highlights the achievement that structure classes were ecologically meaningful.These strong relationships to one another reveal the legitimacy of the structure class definitions and also corroborate the accuracy assessment of the final 2013 forest structure map.

Change Vector Analysis and Temporal Continuity Rules
CVA is advantageous for mapping change across a long time period, because areas of no change are not reclassified repeatedly, not only saving time, but also avoiding compound error in those places that have not experienced change throughout the time series [48].We over-predicted the potential change in each yearly image by choosing thresholds of change that included, on average, 44% of the pixels in the study area.However, we can be confident that the unchanged pixels used for training the changed pixels were, indeed, accurately classified, and those unchanged pixels that were included within the threshold of change were likely reclassified to their original class.
Looking at changes in forest structure over time allows us to improve each yearly classification using logic and a priori knowledge to fix minor errors within the time series.We applied the temporal continuity rules that we developed for our study area to the entirety of the time series in order to improve upon the initial accuracies of the forest structure classes.Through this process, we were able to improve our ability to separate Open from Stand Initiation and increased the overall accuracy of the master 2013 6-class prediction (initial overall accuracy-68%; final overall accuracy-72%) (Tables 1 and 3) The overall accuracies for the 6-class predictions for 2005, 1995, and 1975 were not very high relative to the 2013 prediction (Tables 1 and 3).The relatively lower accuracies for the maps created through the CVA process might be a result of the lower spatial resolution of the images used to identify the validation points.The 2013 validation points were collected from 1-m NAIP imagery within GoogleEarth, while the 2005 and 1995 validation points were collected from 30-m (Landsat) imagery within GoogleEarth (with comparisons to more recent 1-m NAIP images), and the 1975 validation points were collected from historic aerial photographs at 1:16,000-to 1:40,000-scale.
The LandsatLinkr (LLR) program was invaluable for this project.The processes through which LLR calibrated MSS data to TM data and OLI data to ETM+ data aided in harmonizing the four different types of imagery to be comparable across all 44 years, created a Tasseled Cap Wetness component for the MSS data, and allowed us to use all three Tasseled Cap components for the CVA.Additionally, without the automated methods of LLR, the pre-processing time would have increased the time it took to perform all of the steps in preparing each yearly image for classification through the change vector analysis process by orders of magnitude.

Conclusions
There are three novel outcomes of this study.First, we used standard CVA methods to classify forest structure across a diverse forest landscape from 1972 to 2015 using only M2O and terrain data, explicitly distinguishing forest structure classes for single dates.Second, we were able to apply temporal continuity rules to the time series itself to improve the accuracy of each individual yearly forest structure classification.Finally, we demonstrated that the MSS data (pre-1984) are important within these time series and believe that the USGS goal of calibration and radiometric enhancements of MSS data is timely and very valuable.We argue that there are significant advantages to using the MSS data even if it is lower quality than the TM data and beyond.Using MSS data within long-term time series studies can add critical temporal information with which to learn about our landscapes and how they are changing over time.
The process of using information from the time series to inform each year of any classification can be improved upon by refining the rulesets to each specific classification, paying particular attention to the first and last years of the time series.Additionally, performing accuracy assessment after each ruleset could perhaps identify rules that need adjusting.Finally, incorporating spatial rulesets along with the temporal rulesets might smooth out more "salt and pepper" pixels and clean edges where gradual changes are happening, but should be used with caution, as such smoothing can easily result in the loss of real information.
We classified forest structure in this manner, however, other forest and land use classifications could benefit from this process, particularly in cases where land management practices have changed within the time period of the M2O.Forest managers can utilize this process to estimate the amount of increase or decrease of harvestable trees.Wildlife biologists can utilize this process to estimate the amount of viable habitat that has been lost or gained.Land managers can utilize this process to provide evidence of effects of management practices.We should be using the entire length of the M2O dataset and this process is one way to do that.The Landsat program will continue to acquire Earth imagery into the foreseeable future.Time series analyses that can be performed with these data will only grow in importance for land managers as land management practices adapt and change according to current polices.Indeed, Landsat time series might be used to address additional important issues that cannot be addressed with yearly estimates.The accurate classification of images within these time series can be a key to the success of future land management policies.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/9/4/157/s1,Table S1: Every Landsat image used in the production of yearly composite images for our study area; Dictionary S1: Definitions of Forest Structure Classes; Table S2 S26: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 6-class classification for 2013; Table S27: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 2013 (Open and Stand Initiation combined as one class); Table S28: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 2013 (Advanced Regeneration and Mature combined as one class); Table S29: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 4-class classification for 2013; Table S30: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 2013 (Open, Stand Initiation, and Thin combined as one class); Table S31: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 2013 (Thin, Advanced Regeneration, and Mature combined as one class); Table S32: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 2013 (Advanced Regeneration and Mature combined as one class versus everything else); Table S33: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 2013 (Alpine versus everything else); Table S34: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 6-class classification for 2005; Table S35: Full error matrix of estimated area proportions for the FINAL (time-series-informed  S40: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 2005 (Advanced Regeneration and Mature combined as one class versus everything else); Table S41: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 2005 (Alpine versus everything else); Table S42: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 6-class classification for 1995; Table S43: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 1995 (Open and Stand Initiation combined as one class); Table S44: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 1995 (Advanced Regeneration and Mature combined as one class); Table S45: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 4-class classification for 1995; Table S46: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 1995 (Open, Stand Initiation, and Thin combined as one class); Table S47: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 1995 (Thin, Advanced Regeneration, and Mature combined as one class); Table S48: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 1995 (Advanced Regeneration and Mature combined as one class versus everything else); Table S49: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 1995 (Alpine versus everything else); Table S50: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 6-class classification for 1975; Table S51

Figure 1 .
Figure 1.Location map for the study area within Montana, USA.Boundaries of WRS-1 (MSS) scenes and WRS-2 (TM, ETM+, and OLI) scenes are shown over the study area.Scene boundaries are labeled as path/row.

Figure 1 .
Figure 1.Location map for the study area within Montana, USA.Boundaries of WRS-1 (MSS) scenes and WRS-2 (TM, ETM+, and OLI) scenes are shown over the study area.Scene boundaries are labeled as path/row.

Figure 2 .
Figure 2. Different levels of the initial master 2013 classification from 2 classes to 6 classes.Only the initial 2013 classification had 6 classes-the initial classification of every other yearly classification started with 5 classes because Open and Stand Initiation were too similar to separate at this stage of the study.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.See Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes for detailed definitions of each class.Each level was produced with unique combinations of binary classifications (see Supplementary Material, Tables S2-S7 for further details).

Figure 2 .
Figure 2. Different levels of the initial master 2013 classification from 2 classes to 6 classes.Only the initial 2013 classification had 6 classes-the initial classification of every other yearly classification started with 5 classes because Open and Stand Initiation were too similar to separate at this stage of the study.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.See Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes for detailed definitions of each class.Each level was produced with unique combinations of binary classifications (see Supplementary Material, Tables S2-S7 for further details).

Figure 3 .
Figure 3. Examples of temporal continuity rules applied to a variety of 15-year time series.There were eight temporal continuity rules applied to the time series in order to improve accuracy and separate difficult to distinguish classes such as Open versus Stand Initiation: (1) Gap fill #1, (2) Continuity fill, (3) Smoothing, (4) Separate Open from Stand Initiation, (5) Gap fill #2, (6) Separate Open from Stand Initiation Year 1, (7) Alpine based on elevation, and (8) Non-forest based on elevation (each rule is described in detail in Section 2.7 and code/logic used for each rule can be seen in Supplementary Material: Logic S1: Logic Used for Time-Series-Informed Rules).

Figure 4 .
Figure 4. Different levels of the final 1975, 1995, 2005, and 2013 classifications tested, from 2 classes to 6 classes.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.See Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes for detailed definitions of each class.Each subsequent level was created by combining classes from the 6class prediction.Whereas the initial classifications for every year except 2013 had only 5 classes (Open and Stand Initiation were combined as one class), all of the final yearly predictions were classified into 6 classes.The 5-class through 2-class predictions included different combinations of classes (i.e., Thin combined with the open classes versus Thin combined with the forest classes in the 3-class levels) to allow for the most accurate results.

Figure 3 .
Figure 3. Examples of temporal continuity rules applied to a variety of 15-year time series.There were eight temporal continuity rules applied to the time series in order to improve accuracy and separate difficult to distinguish classes such as Open versus Stand Initiation: (1) Gap fill #1, (2) Continuity fill, (3) Smoothing, (4) Separate Open from Stand Initiation, (5) Gap fill #2, (6) Separate Open from Stand Initiation Year 1, (7) Alpine based on elevation, and (8) Non-forest based on elevation (each rule is described in detail in Section 2.7 and code/logic used for each rule can be seen in Supplementary Material: Logic S1: Logic Used for Time-Series-Informed Rules).

Figure 3 .
Figure 3. Examples of temporal continuity rules applied to a variety of 15-year time series.There were eight temporal continuity rules applied to the time series in order to improve accuracy and separate difficult to distinguish classes such as Open versus Stand Initiation: (1) Gap fill #1, (2) Continuity fill, (3) Smoothing, (4) Separate Open from Stand Initiation, (5) Gap fill #2, (6) Separate Open from Stand Initiation Year 1, (7) Alpine based on elevation, and (8) Non-forest based on elevation (each rule is described in detail in Section 2.7 and code/logic used for each rule can be seen in Supplementary Material: Logic S1: Logic Used for Time-Series-Informed Rules).

Figure 4 .
Figure 4. Different levels of the final 1975, 1995, 2005, and 2013 classifications tested, from 2 classes to 6 classes.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.See Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes for detailed definitions of each class.Each subsequent level was created by combining classes from the 6class prediction.Whereas the initial classifications for every year except 2013 had only 5 classes (Open and Stand Initiation were combined as one class), all of the final yearly predictions were classified into 6 classes.The 5-class through 2-class predictions included different combinations of classes (i.e., Thin combined with the open classes versus Thin combined with the forest classes in the 3-class levels) to allow for the most accurate results.

Figure 4 .
Figure 4. Different levels of the final 1975, 1995, 2005, and 2013 classifications tested, from 2 classes to 6 classes.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.See Supplementary Material: Dictionary S1: Definitions of Forest Structure Classes for detailed definitions of each class.Each subsequent level was created by combining classes from the 6-class prediction.Whereas the initial classifications for every year except 2013 had only 5 classes (Open and Stand Initiation were combined as one class), all of the final yearly predictions were classified into 6 classes.The 5-class through 2-class predictions included different combinations of classes (i.e., Thin combined with the open classes versus Thin combined with the forest classes in the 3-class levels) to allow for the most accurate results.

Table 1 .
Overall accuracies and standard errors (SE) of the initial 2013, 2005, 1995, and 1975 classifications.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.Only the 2013 classifications separated Open from Stand Initiation.The 5-and 6-class classification accuracies were notably lower than those of the classifications with fewer classes, likely because of the difficulty in distinguishing similar cover types from one another, i.e., Open versus Stand Initiation or Advanced Regeneration versus Mature.Validation points for 2013, 2005, and 1995 were interpreted from GoogleEarth imagery.Validation points for 1975 were interpreted from digital aerial photographs.Full error matrices can be seen in Supplementary Material, Tables S2-S25.
-class classifications for 2005 and 1995 were higher than the initial 5-class overall accuracies.The overall accuracy for the 1975 6-class classification was 43%, however, all of the 1975 4-class through 2-class overall accuracies increased from the initial classifications.Most of the errors in the 1975 6-class classification were (1) stand initiation being classified as open and (2) thin and advanced regeneration being classified as mature (see Supplementary Material, S50-S57).

Table 3 .
Overall accuracies and standard errors for final classifications of 2013, 2005, 1995, and 1975.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.We observed increased overall accuracies and decreased or equivalent standard errors for most of the time-series-informed classifications compared to the initial classifications (

Forests 2018, 9 ,
x FOR PEER REVIEW 13 of 21 final 2015 classification contains some errors (described below in Section 4.2) that are reflected in the last segment of these class trajectories.

Figure 6 .
Figure 6.Time series trajectories of six forest structure classes from 1972 to 2015.Percentage of study area of each class for each year is shown on the Y-axis.For instance, Mature has lost area overall, starting with 42% of the study area in 1972 and ending with 35% of the study area in 2015.The transition from MSS to TM data is shown at the 1984 mark.

Figure 6 .
Figure 6.Time series trajectories of six forest structure classes from 1972 to 2015.Percentage of study area of each class for each year is shown on the Y-axis.For instance, Mature has lost area overall, starting with 42% of the study area in 1972 and ending with 35% of the study area in 2015.The transition from MSS to TM data is shown at the 1984 mark.

Figure 7 .
Figure 7. Sankey diagram of the state transition from one class to another of the six forest structure classes from 1972 to 2015.Y-axes show the percent of each class that covers the study area in 1972 and 2015, respectively.For instance, Mature has lost area overall, starting with 42% of the study area in 1972 and ending with 35% of the study area in 2015.If we combine the forested classes (Mature, Advanced Regeneration, and Thin, and Stand Initiation), we observe that there has been a slight increase of forest cover from 1972 to 2015, starting with 67% of the study area and ending with 74% of the study area.

Figure 7 .
Figure 7. Sankey diagram of the state transition from one class to another of the six forest structure classes from 1972 to 2015.Y-axes show the percent of each class that covers the study area in 1972 and 2015, respectively.For instance, Mature has lost area overall, starting with 42% of the study area in 1972 and ending with 35% of the study area in 2015.If we combine the forested classes (Mature, Advanced Regeneration, and Thin, and Stand Initiation), we observe that there has been a slight increase of forest cover from 1972 to 2015, starting with 67% of the study area and ending with 74% of the study area.

Figure 8 .
Figure 8. Changes in forest structure classification from 1972 to 2015 within the southern portion of our study area.Fire scars (red circles in 2007) and clear cuts (blue circle around checkerboard pattern in 2012) are the most obvious visual changes.Insect damage is observed as a gradual decrease of the Mature class over the course of the time series (as seen from 1972 to 2015 in black circles).

Figure 8 .
Figure 8. Changes in forest structure classification from 1972 to 2015 within the southern portion of our study area.Fire scars (red circles in 2007) and clear cuts (blue circle around checkerboard pattern in 2012) are the most obvious visual changes.Insect damage is observed as a gradual decrease of the Mature class over the course of the time series (as seen from 1972 to 2015 in black circles).
) 5-class classification for 2005 (Open and Stand Initiation combined as one class); Table S36: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 2005 (Advanced Regeneration and Mature combined as one class); Table S37: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 4-class classification for 2005; Table S38: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 2005 (Open, Stand Initiation, and Thin combined as one class); Table S39: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 2005 (Thin, Advanced Regeneration, and Mature combined as one class); Table : Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 1975 (Open and Stand Initiation combined as one class); Table S52: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 5-class classification for 1975 (Advanced Regeneration and Mature combined as one class); Table S53: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 4-class classification for 1975; Table S54: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 1975 (Open, Stand Initiation, and Thin combined as one class); Table S55: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 3-class classification for 1975 (Thin, Advanced Regeneration, and Mature combined as one class); Table S56: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 1975 (Advanced Regeneration and Mature combined as one class versus everything else); Table S57: Full error matrix of estimated area proportions for the FINAL (time-series-informed) 2-class classification for 1975 (Alpine versus everything else); Figure S1: Decadal Sankey diagram of the state transition from one class to another of the six forest structure classes from 1972 to 1980, 1980 to 1990, 1990 to 2000, 2000 to 2010, and 2010 to 2015; Logic S1: Logic used for time-series-informed-rules.

Table 2 .
Producer's and User's accuracies for the initial 6-class 2013 classification and 5-class (Open and Stand Initiation combined as one class) 2013, 2005, 1995, and 1975 classifications.A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature.The Open class had the lowest producer's accuracy for the 2015 6-class classification, while the Stand Initiation class had the lowest user's accuracy.When these two classes were combined, their producer's and user's accuracies are notably higher.Validation points for 2013, 2005, and 1995 were interpreted from GoogleEarth imagery.Validation points for 1975 were interpreted from digital aerial photographs.Full error matrices of all the final classifications can be seen in Supplementary Material, TablesS2-S7.

Table 1 )
, and those that did not increase in overall accuracy decreased by less than 1%.Validation points for 2013, 2005, and 1995 were interpreted from GoogleEarth imagery.Validation points for 1975 were interpreted from digital aerial photographs.Full error matrices can be seen in Supplementary Material, Tables S26-S57.

Table S8 :
: Full error matrix of estimated area proportions for the INITIAL 6-class classification for 2013; Table S3: Full error matrix of estimated area proportions for the INITIAL 5-class classification for 2013 (Open and Stand Initiation combined as one class); Table S4: Full error matrix of estimated area proportions for the INITIAL 5-class classification for 2013 (Advanced Regeneration and Mature combined as one class);Table S5: Full error matrix of estimated area proportions for the INITIAL 4-class classification for 2013; Table S6: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 2013; Table S7: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 2013; Full error matrix of estimated area proportions for the INITIAL 5-class classification for 2005 (Open and Stand Initiation combined as one class); Table S9: Full error matrix of estimated area proportions for the INITIAL 4-class classification for 2005; Table S10: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 2005 (Open, Stand Initiation, and Thin combined as one class); Table S11: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 2005 (Thin, Advanced Regeneration, and Mature combined as one class); Table S12: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 2005 (Advanced Regeneration and Mature combined as one class versus everything else); Table S13: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 2005 (Alpine versus everything else); Table S14: Full error matrix of estimated area proportions for the INITIAL 5-class classification for 1995; Table S15: Full error matrix of estimated area proportions for the INITIAL 4-class classification for 1995; Table S16: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 1995 (Open, Stand Initiation, and Thin combined as one class); Table S17: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 1995 (Thin, Advanced Regeneration, and Mature combined as one class); Table S18: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 1995 (Advanced Regeneration and Mature combined as one class versus everything else); Table S19: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 1995 (Alpine versus everything else); Table S20: Full error matrix of estimated area proportions for the INITIAL 5-class classification for 1975; Table S21: Full error matrix of estimated area proportions for the INITIAL 4-class classification for 1975; Table S22: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 1975 (Open, Stand Initiation, and Thin combined as one class); Table S23: Full error matrix of estimated area proportions for the INITIAL 3-class classification for 1975 (Thin, Advanced Regeneration, and Mature combined as one class); Table S24: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 1975 (Advanced Regeneration and Mature combined as one class versus everything else); Table S25: Full error matrix of estimated area proportions for the INITIAL 2-class classification for 1975 (Alpine versus everything else); Table