Differentiation of Alternate Harvesting Practices Using Annual Time Series of Landsat Data

Abstract

Sustainable forest management practices allow for a range of harvest prescriptions, including clearcut, clearcut with residual, and partial or selective cutting, which are largely distinguished by the amount of canopy cover removed. The different prescriptions are aimed to emulate natural disturbance, encourage regeneration (seed trees), or offer other ecosystem services, such as the maintenance of local biodiversity or habitat features. Using remotely sensed data, stand-replacing disturbance associated with clearcutting is commonly accurately detected. Novel time series-based change detection products offer an opportunity to determine the capacity to detect and label a wider range of harvest practices. In this research, we demonstrate the capacity of time series imagery, spectral metrics, and related attributed change products, to distinguish between different harvesting practices over a study area in central British Columbia, Canada. Producer’s accuracy of harvest attribution was 79%, with 93% of harvest blocks >5 ha accurately identified. In relation to the amount of canopy cover removed, clearcut harvesting was the most accurately classified (84%), followed by clearcut with residual (79%), and partial cut (64%). Applying detailed spectral metrics derived from Landsat data revealed clearcut and partial cuts to be spectrally distinct. The annual nature of the Landsat time series also offers spatial harvest information within typical, often decadal, forest inventory update cycles. The statistically significant (p < 0.05) relationship between harvest practices and Landsat spectral information indicates a capacity to add increased attribution richness to remote sensing depictions of forest harvest.

1. Introduction

Approximately 65% of Canada is occupied by forest-dominated ecosystems [1], of which ~350 million hectares are represented by trees and other wooded land [2], with the remainder dominated by lakes and wetlands. Knowledge on the extent, condition, and management of forest resources is critical to meet both national [3,4] and international [5] expectations of sustainable forest management and to support reporting requirements. In addition, the magnitude of harvest events, and the type of harvest (clearcut or selective cut), has implications to carbon budget modeling. As such, information on the extent and type of forest harvesting is captured and maintained through spatial forest inventories to assist with management and reporting activities. Due to the nature of natural resource stewardship in Canada, provincial and territorial resource management agencies are responsible for managing forest resources. An element of forest resource management is the allocation of timber harvest permits, including the type of harvesting that may be undertaken [6,7]. Given the processes in place, many jurisdictions do not maintain spatial harvest records in a comprehensive and standardized fashion, with statistics often reported aspatially and/or on a jurisdictional basis (see [7]). Traditionally, forest inventories in Canada have been based upon the collection, delineation, and attribution of air photos on a periodic basis, often of a decadal interval. Ground data are used to calibrate and ensure the consistency of the information in the forest inventory. Jurisdictional forest inventories are strategic in nature, aimed at ensuring an understanding of the nature and extent of forest resources across a given jurisdiction [8]. At the operational level, new technologies for forest inventory are undergoing investigation [9] with lidar data showing the greatest promise and operational uptake [10]. Increasingly large areas are being flown to acquire lidar, yet provincial and territorial standards have yet to be altered to entrench lidar as the sole or favored data source. In contrast and at a different categorical and spatial scale, space-based remote sensing technology offers an opportunity to capture wall-to-wall information on forest harvesting. The opening of the Landsat archive has enabled the generation of time series informed, annual depictions of forest change by type [11]. Canada conducts a National Forest Inventory (NFI) to produce nationally consistent reporting on forest resources. Undertaken in partnership with provincial and territorial resource management agencies, the NFI is based upon a systematic grid of 2 × 2 km photo plots that are, in turn, supported by a subset of ground plots. The photo plots are analogous in information content to polygonal forest inventory data. Based upon the systematic grid and photo plot representation, a 1% sample of Canada’s forest resources allows for statistically robust reporting. Given the 1% sample and a 10-year update cycle, there are numerous opportunities for the integration of wall-to-wall satellite derived information on forest dynamics [12]. The frequent and consistent acquisition of high quality, spatially precise, calibrated, analysis-ready remotely-sensed imagery, from sensors such as Landsat and Sentinel 2, facilitates improved management and stewardship, which subsequently supports integrated monitoring and reporting information needs [13].

Free and open access to the United States Geological Survey (USGS) Landsat archive [14] has led to an increase in the use of Landsat data for both large area and temporal analyses (e.g., [15,16,17,18]). Landsat imagery has imaging characteristics that allow for detection of human interactions with the landscape [19] based upon spatial, spectral, temporal, and radiometric qualities that enables implementation of robust and transparent analyses [19,20]. Further, as noted by Belward and Skøien [21], the temporal depth of the Landsat archive is unique among Earth observation satellites at the moderate or medium spatial resolution (10–100 m), typically required to capture anthropogenic activities. Utilizing the spatially complete and temporally deep database of Landsat data available for Canada [12] time series analyses have been conducted utilizing the 30 m spatial resolution data available from Landsat Thematic Mapper ™ and Enhanced Thematic Mapper Plus (ETM+) for the period 1984 to 2012 [18]. These data both define baseline knowledge of historic patterns of forest disturbances, as well as portray current forest disturbance and recovery trends that are needed to make informed management decisions. The 1984 baseline offered by Landsat is earlier than most existing spatial data records, hence otherwise unavailable historic information on forest harvesting is now available as well as a present day accounting that can augment jurisdictional records.

Furthermore, an opportunity exists to better understand the forest change attributable to harvest as reported by Hermosilla et al. [18], and to determine the capacity to increase the detail of the current attribution present in the change hierarchy used by Hermosilla et al. [11]. This potential is increasingly important as partial cuts become more commonplace and the ability to report on, or track the outcomes of, different harvest practices is desired. In 2014, only 12% of timber harvest on public land in British Columbia (by area) was clearcut, the remainder being some variation of partial cut technique or clearcut with reserves [22]. This increasing trend is expected to continue as large clearcut blocks are subject to public scrutiny [23]. The ability to distinguish these harvest practices using remotely sensed data presents a potential avenue to governments focused on implementation of sustainable forest management practices.

The science of detecting harvest events from remotely sensed data has evolved over time. Regression modeling using Landsat data by Healey et al. [24] highlights that the shortwave infrared bands and associated metrics are the best indicators for detecting forest removals. Before the opening of the Landsat archive, Kennedy et al. [25] focused on using Landsat metrics to detect the year of change for different harvest methods. While the year of change was accurate (overall accuracy = 90%), the classification of the harvest types was considerably lower (overall accuracy = 74%). Use of the Landsat time series is useful for differentiating forest disturbances, as demonstrated by Schroeder et al. [26], who classified clearcuts and fire events with an overall accuracy of 93%.

In this research, we compared ground truth information of specific harvest practices to disturbances detected using a time series of Landsat imagery [18]. In so doing, we determined the detection capacity for harvests implemented with differing levels of intensity (that is, canopy cover removal as related to a given harvest practice). Using spectral and temporal metrics, we determined the capacity of Landsat time series to relate how harvest practices manifest on a spectral level. We implemented statistical tests on the most descriptive spectral metrics to gain an improved understanding of the significance of the differences exhibited across the differing harvest practices present. We also investigated the accuracy of the time series disturbance classification against the type of harvest practice and assessed the influence of harvest polygon size on classification accuracy.

2. Study Area

The study area is a 50 × 35 km zone centered on the city of Prince George in north-central British Columbia, Canada (Figure 1). The region is situated in the interior plateau of the Western Cordillera, between the Coast and Columbian mountain ranges. Elevation in the area ranges between 550 and 800 meters with a gently rolling terrain. Classified within the Sub-boreal Spruce biogeoclimatic zone, the dominant climax tree species in the area are Picea glauca × engelmannii, Picea glauca, and Abies lasiocarpa. Pinus contorta is also common in mature stands and is the main pioneer species for the region along with Populus tremuloides [27].

Prince George and the surrounding Cariboo developmental region are strongly invested in logging and forestry, accounting directly for 16.6% of all jobs in the region and 2.5% of GDP provincially [28]. Various harvest practices have been applied in the area. Harvest history in the region is reflective of trends occurring across the rest of North America, with clearcuts being the most common and partial or residual prescriptions being more preferential in recent years [29]. Forest management in the region has also been influenced by the recent mountain pine beetle outbreak, resulting in high mortality of the economically important Pinus contorta subsp. contorta [30].

3. Data

3.1. Landsat Time-Series Data

To identify harvest polygons, we utilized the Landsat-derived disturbance information dataset produced by Hermosilla et al. [18] following the Composite2Change (C2C) approach (Figure 2). This approach is based on the production of best-available pixel (BAP) image composites from Landsat imagery by selecting the most suitable observations for each pixel location within a specific date range (August 1, ±30 days) based on the scoring functions defined by White et al. [13]. These functions rank every pixel based on the presence and distance to clouds and their shadows, the atmospheric quality, and the acquisition sensor. These image composites are further refined by removing noisy observations and infilling of data gaps using spectral trend analysis of pixel time series [31], which results in seamless annual surface-reflectance composites for all of Canada from 1984 to 2012, as well as the detection and characterization of forest change events greater than 0.5 ha (minimum map unit). Following the object-based image analysis approach presented in Hermosilla et al. [11], detected changes are attributed to a change type (i.e., fire, harvesting, road, or non-stand-replacing), based on their spectral, temporal, and geometrical characteristics using a Random Forests classifier. The overall accuracy of the change attribution process for all classes is 92%, with harvesting having a producer’s and user’s accuracy of 88%.

3.2. Validation Data

Information on the location, extent and type of harvest activity occurring within the region was compiled by the Pacific Forestry Center (PFC) of the Canadian Forest Service (CFS), and consists of a sample of harvest blocks selected from the Prince George forest region. Data were compiled from a combination of ground site visits (on 8–10 September 2003) and aerial photo interpretation. More details on the reference data are provided in Leckie et al. [32]. Ground visits were given priority at sites that were deemed to be the most likely to be misclassified from aerial photography. For this study, we utilized all reference harvest polygons that intersected with the C2C change objects, with more than 50% overlap, and were greater than the C2C minimum map unit (0.5 ha). The size of the analyzed harvest polygons derived from C2C varied from 0.5 to 96 ha, with an average size of 8.8 ha (σ = 13). A total of 160 reference harvest polygons were used for validation, distributed by harvest type as follows: clearcut (80 polygons), clearcut with residual (47), and partial cut (33). Clearcut implies a canopy closure less than 1% post-harvest. Clearcut with residual indicates a canopy cover between 1% and 10% post-harvest. Finally, partial cut implies a selective harvesting technique resulting in a stand with >10% canopy cover remaining (Figure 3). The average polygon size for clearcuts was 7.8 ha (σ = 10), while clearcuts with residual polygons were 9.5 ha on average (σ = 15), and partial cut polygons were 10.3 ha (σ = 15.6). The spatial distribution of the harvest polygons is shown in Figure 4.

4. Methods

To compare the changes detected from the Landsat time series and the forest harvest data, our approach was as follows. We assessed the capacity of the Landsat time-series-derived metrics to describe and discriminate between the three different harvest practices identified in the reference data. This was undertaken by matching the reference harvest polygons with the change pixels detected from the C2C Landsat time series. Then, we evaluated the C2C attribution of harvest, to assess (i) how accurately the C2C data identified harvest; and then (ii) how that accuracy varied by harvest type.

4.1. Spectral Characterization of Harvesting Practices

We analyzed the capacity of the Landsat time-series spectral metrics to discriminate between the three different harvesting practices. To do so, six spectral metrics, obtained from the time series, were calculated and extracted for each of the analyzed harvest polygons [31]. The six spectral metrics were change persistence (duration), average change magnitude variation, average spectral value post-change, average spectral value pre-change, range of pixel series values, and standard deviation of pixel series values (see full description in Hermosilla et al., [31]). Persistence represents the length in years following the harvest event that spectral values remain at post-change levels.

Magnitude variation refers to the spectral difference of Normalized Burn Ratio (NBR) [33] before and after the harvest event. The NBR is a spectral index used to map burn severity, and is calculated using Landsat TM/ETM+ bands 4 and 7 as follows:

$$NBR=\frac{B4-B7}{B4+B7}$$

(1)

Pre-change represents the average spectral NBR value of the years leading up to the year of greatest change and post-change represents the average NBR value directly following the harvest event. The standard deviation of pixel series values represents the spectral variation of NBR within each harvest practice. The time series range measured using tasseled cap greenness (TCG) [34] represents the variability of pixel values for each harvest practice. An ANOVA together with Tukey’s post hoc test were performed on these spectral metrics in order to determine significant differences (α = 0.05, p < 0.05) in their values between the harvesting practices. Additionally, we examined the complete trends of the spectral values centering the trends on the change year [35]. In this way, trends can be compared regardless of the year of disturbance, with trends referenced relative to the number of years before (y − 1, y − 2, etcetera) and after (y + 1, y + 2, etcetera).

4.2. Change Attribution

We evaluated the accuracy of the C2C change type attribution using the reference harvest polygons, and summarized these results in a modified confusion matrix [36]. As the reference data only characterizes harvest, we were not able to evaluate the accuracy of the other C2C change types. In cases of multiple C2C disturbance type classifications within the harvest polygon, the most frequently occurring class in terms of overlapping area was chosen as representative. Any change attribution value other than harvesting was considered as misclassified. Producer’s accuracies were estimated by harvest practice. Lastly, the influence of the size of the harvesting polygons on the classification accuracy was analyzed.

5. Results

5.1. Discriminating Harvest Type

Figure 5 shows the distribution of the Landsat-derived spectral metrics for the three harvest practices (i.e., clearcut, clearcut with residual, and partial cut).

Table 1. Statistically significant different groups determined with Tukey’s post hoc test (p < 0.05).

Metric	Clearcut	Clearcut with Residual	Partial Cut
Change persistence	a	a	b
Average change magnitude variation NBR	b	ab	a
Average spectral value post-change NBR	b	a	a
Average spectral value pre-change NBR	a	a	a
Range of pixel series values TCG	b	ab	a
Standard deviation of pixel series values NBR	a	a	b

Homogenous groups per metric are identified using the same letter. NBR, Normalized Burn Ratio; TCG, tasseled cap greenness.

shows the significant differences between the spectral values for the different harvest types as determined with the ANOVA. The average pre-change NBR showed no statistical difference, and confirmed uniformity of the forest stands prior to harvest. Thus, significant differences in spectral metrics post-harvest were assumed to be the result of the harvest method applied and not influenced by different pre-harvest forest conditions.

The results of the ANOVA indicated that partial cuts were significantly different from clearcuts, and also different from clearcuts with residual. The average spectral magnitude variation of the change mirrors the expected intensity gradient in harvest practices. As shown in Figure 6, clearcuts show the largest spectral change magnitude variation, followed by clearcuts with residual, with partial cuts having the lowest spectral change values. Statistically, clearcut and partial cut categories are significantly different from each other, and the clearcut with residual statistically similar to either clearcut or partial cut, depending on the metric. The average spectral value post-change NBR complements the magnitude variation metric by indicating that more intense harvesting methods lead to significantly lower spectral values directly after the harvest. These intense changes with low spectral values are more likely to have a steeper recovery trend over the following years compared to less intense changes as shown by the persistence.

The change persistence (duration) shows significant differences between partial cut and the other harvest practices. Persistence for partial cuts ranged from 1 to 5 years, resulting in the recovery trend that is less steep and thus comparatively slower. All three harvest practices have a median persistence of 1 year, yet the partial cuts had a noticeably longer range of durations that can extend up to 5 years. The 5-year range reflects successive cuts removing mature trees to attain an even-aged new stand, which is the defining feature of partial cuts. The relative spectral heterogeneity of recovery time for partial cuts, combined with its criteria of being over 10% canopy cover (which may be difficult to accurately estimate in the field) explains the comparatively weaker classification percentage.

Another explanatory metric which detects difference between harvesting practices is the TCG. The results clearly show large magnitude reductions in TCG values as a consequence of clearcuts. This is in contrast to the smaller magnitude reductions caused by partial cuts. As a result, partial cuts are less likely to be classified as harvesting than clearcuts as the magnitude of the change is comparatively small and likely to be confused with other types of disturbance or missed completely. Once again, the smaller magnitude change can partially be attributed to the variable intensity (anything above 10% cover remaining) of the partial cut.

The mean of NBR values for each of the harvesting practices is shown in Figure 6a–c. Figure 6d shows the average NBR trajectory for each harvest type.

5.2. Change Attribution

Overall accuracy for all reference harvest polygons correctly classified as harvesting was 78.8% (

Table 2. Accuracy of C2C (Composite2Change) harvest attribution.

	Reference Data
	Clearcut	Residual Cut	Partial Cut	Total
Harvest	67	38	21	126
Non-stand replacing	7	5	9	21
Road/infrastructure	6	4	3	13
Fire	0	0	0	0
Total	80	47	33	160
Producer’s accuracy (%)	83.8	79.2	63.6	78.8
Omission error (%)	16.2	20.8	36.4

).

Table 3. Producer’s accuracy, by object size.

Size	Polygons Misclassified	Accuracy (%)
>5 ha	5 of 70	93.0
2–5 ha	14 of 41	66.0
<2 ha	18 of 53	66.0s

displays the attribution accuracy of the harvest polygons based on size. The results showed that 93% of harvest polygons over 5 ha were accurately classified. As the size of the polygons decreases, the classification errors increase, leading to an accuracy of 66% for polygons less than 5 ha. Of the different harvesting practices identified in the reference data, clearcut harvesting is the most accurately classified (83.8%), followed by clearcut with residual with 79.2% of correctly classified. Partial cut forestry practices showed the lowest attribution accuracy (63.6%). Clearcut with residual and partial cut practices had a higher percentage of misclassification as non-stand replacing changes. No harvest polygons were labeled by the C2C algorithm as fire.

6. Discussion

In this research, we assessed the capability of Landsat time series to identify and characterize different harvest practices reported in established harvest records, including clearcut, clearcut with residual, and partial cut. A range of different metrics were analyzed to examine which spectral properties differentiate the harvesting practices. The accuracy of the change attribution varied by harvest practice. Clearcuts were the most accurately classified harvesting practice. Clearcutting removes the majority of vegetation from site and exposes soil. Indeed, we found that the average spectral post-change value was the most statistically important metric for discriminating clearcuts, in keeping with the results of Lambert et al. [37] who identified clearcuts assessing the post-disturbance scenario. A large and fast spectral change in a long-term trajectory, depicted by a steep post change spectral evolution, is commonly representative of a major stand replacing disturbance [38]. Some clearcut plots were erroneously identified as roads, which also are characterized by a steep drop in the spectral response through time. In areas with harvest activity, however, roads are sometimes temporary and can have various surface materials (e.g., dirt, gravel) rendering them difficult to attribute due to lack of uniformity and a spectral similarity to clearcuts [31]. Specific spatial rules to distinguish between roads and other features based on geometric properties, such as linearity, are possible [39]. Moreover, roads tend to persist and will have different post-disturbance spectral characteristics that would allow them to be distinguished from harvesting, given additional observation years.

Clearcut with residual harvest practices were also correctly identified by the C2C algorithm in the majority of cases. The forest stands that are harvested in this manner typically undergo a near complete removal of trees, with between 1% and 10% of canopy cover retained. As a result, this harvesting practice is expected to have a similar spectral response to clearcuts. The remaining trees impact the averaged spectral values after the change event, thus reducing the magnitude, standard deviation, and other metrics. This is best exemplified through the range in greenness values of the time series, which show clearcut with residual practices to be not significantly different from partial cuts. As clearcut residual polygons retain large mature trees to aid in regeneration, the average greenness of the polygon has less deviation from the pre-change spectral response [40]. Depending on the configuration of the residual cover, it may confound the automated detection and attribution of the change event using the Landsat time series.

Partial cuts were the most challenging practice to identify, with lower accuracy relative to clearcuts and residual cuts. This detection and differentiation difficulty is explained by the spectral response of partial cut polygons, which show comparatively low change magnitude and standard deviation values through the time series. The majority of the confusion is from partial cuts being labeled as non-stand replacing changes, which are typically gradual changes in vegetation condition that do not lead to a change in land cover class (i.e., disease, insects, water stress, and decline; [41]). Partial cuts are defined as having >10 percent canopy cover post-harvest. Since the aim of partial cut practices is to mimic low severity natural disturbances [42,43] or maintain biodiversity [44] these changes have a high probability of being associated with non-stand replacing disturbances, such as drought stress or insects. The lower range of temporal TCG values for partial cuts made them distinguishable from clearcuts with this metric. This contrasts with Healey et al. [24], which found that TCG values were weakly correlated with forest removal using basal area as the forest-removal measurement, whereas we used canopy cover.

The analysis of the influence of polygon size on the change attribution accuracy indicates that smaller harvest polygons were more likely to be misclassified. The smallest polygons were 0.5 ha which is equivalent to 6 Landsat-pixels. The Landsat 30 m pixel would likely be influenced by factors outside the harvest plots, as exact delineation of polygon edges is difficult at this spatial resolution. These surrounding spectral values become more influential as polygon sizes decrease, increasing error from a proportional perspective. As expected, the large polygons were classified more reliably, due to less opportunity for edge effects and, more importantly, a greater number of pixels (and related spectral metric values) to sample per harvest event.

7. Conclusions

Despite some minor misclassification, the results of this study showed that partial cuts were significantly different from clearcuts in most spectral metrics characterizing the disturbance or the post change properties. This represents the potential ability to further distinguish between harvesting practices in the disturbance classification hierarchy presented in Hermosilla et al. [31]. Due to the opening of the Landsat archive, the analysis of spectral time series of Landsat data covering large jurisdictions is an increasingly operational approach for verifying harvest activity at a reduced cost compared to field visits or photo digitization approaches for traditional strategic forest inventories. Further differentiation of harvest practices will provide a rich baseline for forest management, analysis, and decision-making moving forward. Our results along with those of others [26] demonstrate the capability of utilizing Landsat time series to further classify forest change events into more specific categories. The methodological opportunities here applied are a direct consequence of the richness of the Landsat archive over Canada, in contrast to many nations and regions globally where the application of these methods would be more limited [12,45]. Spectral metrics derived from Landsat time series provide a quantitative and spatially explicit opportunity to generate refined attribution detail related to forest harvest practices. Traditional photo-based forest inventory techniques, while modernizing with increasingly digital methods, remain limited in their ability to provide a regularly updated and extensive representation of entire management areas. With the combination of the Landsat program’s archival depth, calibrated sensors, analysis ready data, increased acquisition frequency, as well as growing synergies with complementary satellites (e.g., Sentinel 2), new opportunities exist for augmenting and developing forest inventories from regional to national scales.