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

Abstract

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.

1. 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 threatened or endangered species give management conflict. A time series analysis of forest structure can provide empirical evidence of loss or gain of different ages of forest, which sets the backdrop for discussion concerning current and future management policies.

Most successful forest structure mapping has been performed utilizing Light Detection and Ranging (lidar) data [12,13,14,15,16,17]. Lidar data are acquired by using pulsed laser measurements to provide a 3-dimensional point cloud from which precise measures of tree size and stand structure can be derived. However, lidar data are prohibitively expensive and difficult to process and store, especially for large areas. In addition, there is no long-term archive of lidar data with which to perform a time series analysis. This is where the utility of easily accessible Landsat data, and more specifically, the temporal length of the entire Landsat archive dataset, surpasses other datasets.

The Landsat archive is a priceless Earth data repository that is increasingly used for long-term, time-series-based change analyses of our environment, especially since the United States Geological Survey (USGS) began providing the data free of charge (1 October 2008; https://landsat.gsfc.nasa.gov/landsat-7) [18,19]. Most recent classification-based time series studies, however, have not used the entire Landsat archive length; rather they have utilized annual data starting with the Thematic Mapper (TM) sensor through the Enhanced Thematic Mapper Plus (ETM+) and the Operational Land Imager (OLI) sensors that are currently in orbit [20,21,22,23,24]. Although providing useful insights on temporal changes in our environment, these previous studies have not taken advantage of an additional 12 years of Multispectral Scanner (MSS) data from 1972 to 1984. Some studies have taken advantage of the entire Landsat data archive (MSS to OLI or M2O) [12,25], but these analyzed continuous variable data rather than classified thematic data and utilized the M2O time series to inform a single date analysis. Forest-based studies in the United States could greatly benefit from analysis of 1970s data because the National Forest Management Act of 1976 had a profound effect on how federal forests are managed in the United States. We acknowledge that the MSS data have some weaknesses—including, fewer bands and narrow spectral resolution, coarser spatial resolution, lower radiometric resolution, and large gaps in data—but we believe these data contain considerable amounts of information that can be incorporated in time series change analyses, and as such should be exploited when reasonably possible.

A substantial amount of image pre-processing needs to be performed on each scene in order to produce an M2O time series for analysis. LandsatLinkr, a program developed by the Laboratory for Applications of Remote Sensing in Ecology at the USDA Forest Service, PNW Research Station, and Oregon State University, Department of Forest Ecosystems and Society, provides an efficient and automated system “for processing large volumes of Landsat imagery and building long spectrally consistent chronologies across MSS, TM, ETM+, and OLI sensors” [9]. This system not only automates image pre-processing, but also spectrally calibrates MSS images to TM images and OLI images to ETM+ images so that spectral indices can be compared among the four image sensor types.

Our objective was to use the M2O dataset to accurately classify forest structure in a diverse forest landscape in northwestern Montana. Only Landsat imagery and several ancillary data layers were utilized to classify forest structure, reducing the overall costs as compared to a lidar-based study. We observed changes in annual forest structure from 1972 to 2015 and used the temporal information within the time series to improve upon each year of the structure classifications. Additionally, trends in growth and recovery from changes such as harvest, fire, and pests were examined. Finally, we demonstrated the importance of utilizing pre-1984 Landsat data (i.e., MSS) for long-term change studies.

2. Materials and Methods

2.1. Study Area

Our study area covers approximately 3.6 million ha in northwestern Montana (Figure 1). It is a mountainous region ranging from 549 to 3188 m in elevation with a variety of alpine, open, and mixed forest types. These complex forests have gone through a variety of changes since the early 1970s, 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).

2.2. Data Collection

2.2.1. Landsat Imagery

We downloaded 684 Landsat images including MSS (1972–1984) (230), TM (1984–2012) (252), ETM+ (1999–2015) (164), and OLI (2013–2015) (38) summer images (June through September) from the USGS (for the full list of images used in the production of each yearly composite image, see Supplementary Material, Table S1). Images with clouds were included, because the pre-processing procedures were intended to fill those gaps with surrounding data (either temporally or spatially within the annual season).

2.2.2. 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.

2.2.3. 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).

2.3. 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.

2.4. 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].

2.5. 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).

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].

2.6. 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:

$$\Delta \mathrm{Magnitude}=\sqrt{{\left({\mathrm{B}}_1-{\mathrm{B}}_2\right)}^2+{\left({\mathrm{G}}_1-{\mathrm{G}}_2\right)}^2+{\left({\mathrm{W}}_1-{\mathrm{W}}_2\right)}^2},$$

(1)

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.

2.7. 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):

• 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;
• 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);
• Smoothing: looked for single occurrences within the time series (i.e., if 1973 = Thin, 1974 = Open/Stand Initiation, and 1975 = Thin, fill the 1974 pixel with Thin);
• Separate Open from Stand Initiation: based on whether previous years were Mature, Advanced Regeneration, or Thin (we expect that Open would not follow Mature, Advanced Regeneration, or Thin directly, while Stand Initiation can follow them directly); does not include any changes to 1972;
• Gap fill #2: 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;
• Separate Open from Stand Initiation for the 1972 classification: based on 1973 values;
• Alpine by elevation (any Open pixel ≥ 2500 m was labeled as Alpine);
• 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.

2.8. 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 6-class 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.

2.9. 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² 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).

$$\sum \frac{\mathrm{tree}\mathrm{Basal}\mathrm{Area}\ast \mathrm{DBH}}{\mathrm{Total}\mathrm{Basal}\mathrm{Area}},$$

(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.).

3. Results

3.1. 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. 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.

	Overall Accuracy 2013 (%)	SE 2013 (%)	Overall Accuracy 2005 (%)	SE 2005 (%)	Overall Accuracy 1995 (%)	SE 1995 (%)	Overall Accuracy 1975 (%)	SE 1975 (%)
A vs. O vs. S vs. T vs. R vs. M (6-class)	68	3	---	---	---	---	---	---
A vs. OS vs. T vs. R vs. M (5-class)	73	2	55	3	60	3	58	2
A vs. O vs. S vs. T vs. RM (5-class)	78	2	---	---	---	---	---	---
A vs. OS vs. T vs. RM (4-class)	82	2	78	2	75	2	70	2
A vs. OST vs. RM (3-class)	---	---	90	2	82	2	80	2
AOST vs. RM (2-class)	---	---	93	2	87	2	83	2
A vs. OS vs. TRM (3-class)	91	1	84	2	87	2	83	2
A vs. OSTRM (2-class)	98	<1	97	1	97	<1	99	<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. 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, Tables S2–S7.

Classification	Class Name	Producer’s Accuracy 2013 (%)	User’s Accuracy 2013 (%)	Producer’s Accuracy 2005 (%)	User’s Accuracy 2005 (%)	Producer’s Accuracy 1995 (%)	User’s Accuracy 1995 (%)	Producer’s Accuracy 1975 (%)	User’s Accuracy 1975 (%)
6-class	A	98	76	---	---	---	---	---	---
	O	36	67	---	---	---	---	---	---
	S	82	50	---	---	---	---	---	---
	T	49	98	---	---	---	---	---	---
	R	54	80	---	---	---	---	---	---
	M	94	63	---	---	---	---	---	---
5-class	A	96	76	96	78	96	63	96	94
	OS	93	75	82	67	80	88	75	84
	T	43	98	32	54	44	53	20	33
	R	54	80	30	42	28	35	8	12
	M	94	61	60	45	58	40	74	33

, 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.

3.2. 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).

3.3. 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. 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 (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.

	Overall Accuracy 2013 (%)	SE 2013 (%)	Overall Accuracy 2005 (%)	SE 2005 (%)	Overall Accuracy 1995 (%)	SE 1995 (%)	Overall Accuracy 1975 (%)	SE 1975 (%)
A vs. O vs. S vs. T vs. R vs. M (6-class)	72	2	62	3	58	3	43	3
A vs. OS vs. T vs. R vs. M (5-class)	75	2	65	3	61	3	57	3
A vs. O vs. S vs. T vs. RM (5-class)	82	2	80	2	77	3	63	2
A vs. OS vs. T vs. RM (4-class)	85	2	83	2	80	2	76	3
A vs. OST vs. RM (3-class)	93	1	91	2	86	2	87	2
AOST vs. RM (2-class)	96	1	94	2	89	2	88	2
A vs. OS vs. TRM (3-class)	90	1	85	2	88	2	88	2
A vs. OSTRM (2-class)	97	<1	97	1	98	<1	100	<1

) 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-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, Tables S50–S57).

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. Producer’s and User’s accuracies for the final 6-class and 5-class (Open and Stand Initiation combined as one class) 2013, 2005, 1995, and 1975 classifications.

Classification	Class Name	Producer’s Accuracy 2013 (%)	User’s Accuracy 2013 (%)	Producer’s Accuracy 2005 (%)	User’s Accuracy 2005 (%)	Producer’s Accuracy 1995 (%)	User’s Accuracy 1995 (%)	Producer’s Accuracy 1975 (%)	User’s Accuracy 1975 (%)
6-class	A	98	71	96	79	96	81	96	100
	O	94	56	82	52	80	61	87	39
	S	72	82	72	80	68	80	0	0
	T	43	100	34	66	46	65	20	50
	R	56	80	28	50	18	50	6	34
	M	96	66	72	48	82	47	91	40
5-class	A	98	71	96	79	96	82	96	100
	OS	94	75	88	69	85	76	84	77
	T	43	100	34	66	46	65	20	46
	R	56	81	28	50	18	50	6	35
	M	96	66	72	48	82	47	91	40

A = Alpine, O = Open, S = Stand Initiation, T = Thin, R = Advanced Regeneration, M = Mature. Producer’s accuracy for Open and user’s accuracy for Stand Initiation have improved for the final 2013 classification compared to the initial classification (Table 2). We observe that Thin and Advanced Regeneration are difficult to distinguish from one another, while the producer’s and user’s accuracies for the Thin, Advanced Regeneration, and Mature classes have either remained the same or increased from the initial classifications (Table 2). 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, Tables S26–S57.

). 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 (Table 2 and Table 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.

4. Discussion

4.1. 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).

4.2. 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, Figure 6 and Figure 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 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.

4.3. 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%) (Table 1 and Table 3) The overall accuracies for the 6-class predictions for 2005, 1995, and 1975 were not very high relative to the 2013 prediction (Table 1 and Table 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.

5. 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.