Next Article in Journal
Consistent Classification of Landsat Time Series with an Improved Automatic Adaptive Signature Generalization Algorithm
Next Article in Special Issue
Mapping Forest Health Using Spectral and Textural Information Extracted from SPOT-5 Satellite Images
Previous Article in Journal
Methodology for Detection and Interpretation of Ground Motion Areas with the A-DInSAR Time Series Analysis
Previous Article in Special Issue
Alpine Forest Drought Monitoring in South Tyrol: PCA Based Synergy between scPDSI Data and MODIS Derived NDVI and NDII7 Time Series
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Landsat Imagery Spectral Trajectories—Important Variables for Spatially Predicting the Risks of Bark Beetle Disturbance

1
Department of Ecosystem Biology, Faculty of Science, University of South Bohemia, Branišovská 31, České Budějovice 37005, Czech Republic
2
Institute of Botany, The Czech Academy of Sciences, Průhonice 25243, Czech Republic
3
Faculty of Environmental Sciences, Czech University of Life Sciences Prague, Kamýcká 129, Praha 6—Suchdol 16521, Czech Republic
4
Department of Ecology, Institute of Entomology, Biology Centre CAS, Branišovská 31, České Budějovice 37005, Czech Republic
5
Institute for Environmental Studies, Faculty of Science, Charles University in Prague, Benátská 2, 12801 Prague 2, Czech Republic
6
College of Earth, Ocean, and Atmospheric Sciences, Oregon State University, 114 Wilkinson Hall, Corvallis, OR 97331, USA
7
Department of Forest Ecosystems and Society, Oregon State University, 321 Richardson Hall, Corvallis, OR 97331, USA
*
Author to whom correspondence should be addressed.
Submission received: 30 May 2016 / Revised: 12 August 2016 / Accepted: 16 August 2016 / Published: 22 August 2016
(This article belongs to the Special Issue Remote Sensing of Forest Health)

Abstract

:
Tree mortality caused by bark beetle infestation has significant effects on the ecology and value of both natural and commercial forests. Therefore, prediction of bark beetle infestations is critical in forest management. Existing predictive models, however, rarely consider the influence of long-term stressors on forest susceptibility to bark beetle infestation. In this study we introduce pre-disturbance spectral trajectories from Landsat Thematic Mapper (TM) imagery as an indicator of long-term stress into models of bark beetle infestation together with commonly used environmental predictors. Observations for this study come from forests in the central part of the Šumava Mountains, in the border region between the Czech Republic and Germany, Central Europe. The areas of bark beetle-infested forest were delineated from aerial photographs taken in 1991 and in every year from 1994 to 2000. The environmental predictors represent forest stand attributes (e.g., tree density and distance to the infested forest from previous year) and common abiotic factors, such as topography, climate, geology, and soil. Pre-disturbance spectral trajectories were defined by the linear regression slope of Tasseled Cap components (Wetness, Brightness and Greenness) calculated from a time series of 16 Landsat TM images across years from 1984 until one year before the bark beetle infestation. Using logistic regression and multimodel inference, we calculated predictive models separately for each single year from 1994 to 2000 to account for a possible shift in importance of individual predictors during disturbance. Inclusion of two pre-disturbance spectral trajectories (Wetness slope and Brightness slope) significantly improved predictive ability of bark beetle infestation models. Wetness slope had the greatest predictive power, even relative to environmental predictors, and was relatively stable in its power over the years. Brightness slope improved the model only in the middle of the disturbance period (1996). Importantly, these pre-disturbance predictors were not correlated with other predictors, and therefore bring additional explanatory power to the model. Generally, the predictive power of most fitted model decreases as time progresses and models describing the initial phase of bark beetle outbreaks appear more reliable for conducting near-future predictions. The pre-disturbance spectral trajectories are valuable not only for assessing the risk of bark beetle infestation, but also for detection of long-term gradual changes even in non-forest ecosystems.

Graphical Abstract

1. Introduction

Forest disturbances caused by wind, fire and insects are prominent forces controlling forest ecosystem dynamics across the world [1]. Insect forest disturbances represent slower events relative to wind and fire, but have strong effects because of irreversible defoliation and subsequent dieback of mature trees [2]. Bark beetles are among the most common insect disturbance agents in temperate coniferous forests, affecting tens of millions of hectares in recent decades [3,4]. Moreover, the affected forest area is predicted to double in the next two decades because of climate change and cultivation of large areas of coniferous plantations outside their natural distribution [5]. Even though bark beetle disturbances are a part of the natural dynamics of coniferous forest [6,7,8] they also represent a significant threat for timber production and certain forest-provided ecosystem services [9]. Therefore, the need for both detection and prediction of bark beetle infestation is increasing.
Bark beetle outbreaks are commonly detected from satellite imagery using spectral indices, which are primarily based on differences and/or ratios between multiple spectral bands [10]. While the first studies of bark beetle infestation used changes in red and near-infrared (NIR) bands (i.e., normalized difference vegetation index (NDVI)), indices based on NIR and shortwave infrared (SWIR) have become the standard, mostly because of their greater sensitivity to forest canopy defoliation caused by bark beetle infestation [11]. The widely used Tasseled Cap (TC) linear transformation has three components: Brightness, Greenness, and Wetness [12,13]. Among these three components, Wetness is often used to identify bark beetle forest disturbances. It is a combination of all Landsat TM or Enhanced Thematic Mapper Plus (ETM+) spectral bands except the thermal band, with the greatest weight given to the SWIR channel (e.g., [14]).
While the use of remote sensing for detection of bark beetle infestation is well established, its application for prediction or risk assessment is rather rare [15,16,17]. Empirical bark beetle risk-assessment models commonly use environmental predictors including topography, climate, or soil types [18,19,20]. While these environmental predictors are important for bark beetle population dynamics and tree vitality, the forest stand attributes like structure, composition and spatial configuration have at least the same importance; however, they differ at various phases of infestation [16,21,22].
Information on tree vigor, which is the crucial factor shaping the risk of infestation, is usually missing in the standard forest inventory data. Earth observation imagery, specifically from the Landsat sensor family, can fill this lack of information by providing data on tree vitality and forest health status [23]. When analyzed as a time series, Landsat can provide information regarding pre-disturbance change events (e.g., windstorm) or gradual changes in forest vitality (e.g., soil depletion or air pollution), which increase the stand’s vulnerability to bark beetle infestation [24,25,26,27]. Besides the tree to stand-level factors, the risk of bark beetle infestation is increased by large-scale drivers, which include climatic variables and/or landscape connectivity [28]. Most of these latter changes are well reflected in spectral indices and can be documented using multi-temporal remote sensing [29]. Free access to the Landsat data archive by the US Geological Survey (USGS) in 2008 accelerated the use of these data for multi-temporal change detection [30]. Subsequently, temporal segmentation algorithms and spectral trajectory processing tools like LandTrendr [31] and TimeSync [32] were developed to extract the main features (e.g., greatest disturbance segment, segment duration, and magnitude) from spectral time series trajectories. These approaches have enabled a wide use of Landsat data for broad scale ecosystem studies [33]. The use of spectral trajectories based on Landsat imagery helped to describe forest disturbance history since 1972 [34], encouraging forest recovery [35], and detecting bark beetle-affected forest stands more precisely [36]. The following studies distinguished different defoliator and bark beetle effects [37], duration and severity of insect disturbances [38], and large-scale drivers of primary forest disturbance events [39]. A similar trajectory-based approach could be used to detect less pronounced pre-disturbance forest state changes which predispose the forest to disturbances [40].
We hypothesize that gradual changes in forest vitality significantly increase forest vulnerability to bark beetle infestation. Therefore our primary objective is to determine if pre-disturbance spectral trajectories are useful predictors of bark beetle infestation. In particular, we examine how explanatory and predictive power of models that include pre-disturbance spectral trajectories change relative to models in which these variables are excluded. For this purpose, we combine pre-disturbance predictors with environmental predictors in a logistic regression model combined with multimodel inference using several sets of models to assess the importance of pre-disturbance spectral trajectories.

2. Data and Methods

2.1. Study Area

The study area (1376 ha, centroid N 48°57.78435′, E 13°28.33847′) is located in the central part of the Šumava Mountains, in the border region between the Czech Republic and Germany, Central Europe (Figure 1). It is part of the Šumava National Park (Czech Republic) and is adjacent to the Bavarian Forest National Park (Germany). The dominant natural vegetation unit is Norway spruce forest (Picea abies L. Karst.) [41,42]. Where the area is not covered by forest, mostly peats and mountain meadows occur. The relief of the study area is created by a peneplain with numerous local hills. The elevation ranges from 700 m to 1453 m with an average elevation of 1093 m.
Over the last 25 years, several outbreak waves of the European spruce bark beetle (Ips typographus L.) have caused extensive spruce forest dieback in this area. The outbreak originated in wind-fallen trees (173 ha) in 1983 and 1984, following two windstorms in the western part of the Bavarian Forest National Park. These fallen and heavily damaged trees provided suitable conditions for the start of the outbreak [43,44]. At the beginning of the 1990s the first infested trees occurred in the Czech region. Because of the subsequent bark beetle dispersal, the Šumava NP administration applied two management approaches to the infested forests in 1995: (1) a small portion of the stands in the core zone was left unmanaged (more than 2.000 ha); and (2) the stands in the surrounding buffer zone were cut to prevent the spread of bark beetles. We analyzed the area in the unmanaged zone where we expect a dominant role of natural predictors of tree infestation. The bark beetle infestation accelerated between 1995 and 1997 when the largest area in the central part of Šumava Mountains was affected [45]. By the end of the year 2000, most of the spruce stands lost all canopy trees and only small patches of living trees persisted, mostly near peat bogs or at forest edges. After the main dieback in the study area the forest regenerated successfully towards a natural mountain spruce forest [46,47].

2.2. Aerial and Satellite Data

Groups of trees infested by bark beetle (grey-attack stage) were delineated from ortho-rectified aerial photographs, which were acquired in 1991 and every year for the period 1994–2000. The smallest polygon was defined as a group of five adjacent dead trees (0.03 ha). Simultaneously, living spruce forest, clear-cuts and non-forest areas (meadows, peat-bogs) were also delineated. To avoid their influence on values of mixed pixels, we excluded all pixels up to 100 m to clear-cuts and non-forest areas. The resulting land cover polygon layers were transformed to UTM Zone 33 North coordinate system and overlaid on a 30 m × 30 m grid to make them compatible with Landsat imagery. Pixels were defined as declined when the dieback was higher than 50% of its area. We used 16 Landsat TM images from 1984 to 1999 (Table 1). For some years (1988, 1994), multiple images were used to cover no data areas because of clouds and their shadows. To avoid phenological differences among the images, only data acquired during the high summer season were used (nominally July and August with an exception in 1999). The imagery was obtained from the USGS data archive (glovis.usgs.gov). All images were radiometrically normalized to a common image (acquired on 4 August 1990) using the MADCAL approach of Canty et al. [48], following Schroeder et al. [49].

2.3. Spectral Trajectories

Based on the study of Hais et al. [40], which shows the potential of spectral trajectories to describe forest vigor decline before bark beetle infestation, we used this approach to test the significance of these predictors for modelling bark beetle infestation risk. The abovementioned approach uses linear regression to fit multi-temporal spectral data. In our study, we consider all three TC components (Wetness, Brightness and Greenness), though mainly Wetness and Brightness have been shown as the appropriate multi-spectral vegetation indices to describe forest structure [50] and forest health in relation to insect disturbance [14,51,52,53]. To obtain the pre-disturbance spectral trajectories, we calculated the linear regression slope of each TC component from 1984 to one year before the year for which the model was calibrated. The slope was calculated for each 30 m × 30 m pixel except those already infested in previous years. Slopes for non-affected stands were calculated from the same years as the affected ones in the one-year model. An example of Wetness slope is shown in Figure 2. These slopes, which represent historical stand condition (e.g., stress, stability, growth), were included as a model predictor together with commonly used environmental predictors in models of bark beetle infestation risk.

2.4. Environmental Predictors of Bark Beetle Infestation

Despite a large amount of environmental predictors tested for the risk of bark beetle infestation [19,21], we selected representatives of each group (see Table 2) to cover all types of possible influences. Forest stand susceptibility to bark beetle infestation is determined by forest condition itself or abiotic environmental factors, which in turn influence forest vigor. One of the most significant forest attributes for bark beetle dispersal is distance to the nearest infested forest [21]. Despite just 100 m buffer zones of forest in each year, it is still possible to test the distance to the nearest declined forest inside these buffer zones. In addition, we tested the direction to the nearest declined forest in the previous year. The direction azimuth value was transformed to eight 45° wide classes (N, NE, E, ES, S, SW, W and NW). Distance to the nearest clear-cut area and distance to the nearest natural forest edge were taken as a static variable from aerial photograph acquired in 1991. We hypothesize that the clear-cut edges might ease the entry of bark beetles into the stand while the natural forest edges represent a physical barrier for the beetles. Moreover, the trees in forest edge adjacent to clear-cut may suffer from water stress due to higher amount of solar radiation [20,54]. The degree of forest “naturalness” was calculated as a difference between current forest inventory and potential natural tree species composition [55], derived from modified Czech typological maps (natural forest types [56]), which followed a similar concept as potential natural vegetation [57]. Stand age represents an important attribute and we used age of dominant (>50%) tree layer because bark beetles preferentially attack trees older than 60 years [44]. If any tree layer does not reach 50%, we calculated age as a mean value from the first and second dominant tree layer age. Similarly, the dominant tree layer was used for determining the tree density. In fact, there was mostly one dominant tree layer and one small layer of young trees. Soil is characterized by edaphic category, using the Czech typological system. It combines soil texture, nutrient richness and soil moisture [58]. The age, tree density and edaphic category are from the Forest Management Institute (FMI) database (www.uhul.cz). The geology parameters are from Czech Geological Survey maps (www.geology.cz). We calculated the distance to two merged geological classes (7th moor, peat and 8th sediment), because we anticipated that the distance to these geology classes reflect wet or aquatic conditions or water supply. Altitude and other topographic variables were derived from a digital elevation model (DEM) with high spatial resolution (5 m) and vertical error of 1 m (DMR 4G from ČÚZK, www.cuzk.cz). We have used the heat load index (HLI) which describes the variation in potential solar radiation and increasing temperature on SW slopes [59], and the area solar radiation (ASR12), which is a cumulative value of potential solar radiation for entire year (watt hours per square meter—Wh/m2) with maximum value on south slopes. Topographical wetness index (TWI), reflecting potential wetness was calculated using SAGA GIS 2.1.0. We expect that water availability may positively influence the tree resistance against the bark beetle infestation.

2.5. Statistical Models

Bark beetle dispersal results in gradual infestation of adjacent stands from year to year, which carries the possibility of changes in predictor significance during the outbreak period [21]. To assess predictor significance changes, we built an individual statistical model of bark beetle infestation for every year within the studied period. The model for a given year included declined and living forest limited by the 100 m buffer zone around declined forest in the previous year because more than 90% of affected trees for every year were located in this zone. Because of gaps in data on bark beetle dispersal for 1992–1993, we used a 250 m buffer zone between the years 1991 and 1994.
To assess if pre-disturbance spectral trajectories improve the ability of a model to predict bark beetle infestation, we first composed a maximum model that contains all environmental variables plus all spectral trajectories (slopes). We used the square root transformation of the ASR12 variable to get values comparable in order to those of the other variables. To assess potential collinearity between any of the candidate predictors, we looked at the correlations between all combinations of variable pairs, as well as variance inflation factors (VIFs), calculated from the corvif function, which quantifies the severity of multicollinearity [61]. The highest correlation we found was between Brightness and Greenness slopes, ranging between −0.89 and −0.97 for various years. Also, Greenness slope achieved the highest VIF value every year, ranging between 80 and 130, but reaching as high as 430 in one year. A strategy to avoid multicollinearity is to remove the predictor with the highest VIF value and recalculate the VIFs, repeating this process until all VIFs are smaller than a threshold value commonly set to 3–10 [61]. Both approaches suggest removing Greenness slope from our subsequent analyses. Doing this, maximum observed correlations between any two remaining variables reach 0.6–0.7, and all VIFs are mostly below 3, with one or two VIFs reaching 3.1–3.5 in some years. Thus, except for Greenness slope, we keep all explanatory variables for the subsequent statistical modeling.
Since we fix three explanatory variables as mandatory parts of any model (DEM, Dist-source, ASR12), 12 other explanatory variables are left as candidate predictors. From the previous analyses we found that these three variables have relatively stable contribution to the models across years. While the DEM has relatively low importance, the latter two variables played dominant roles for models. Because these variables do not show any trend between the years, we decided to fix them. Using the all-subset approach, we investigate the effect of the pre-disturbance variables by running a total of M = 212 = 4096 candidate models.
Since our response variable is binary (30 m × 30 m forest on pixel is declined or not), we use binary logistic regression for fitting a regression curve. Our data for years 1995 and 2000 suffer from the problem of complete separation, meaning that a predictor or a subset of predictors is associated with only one outcome value when the predictor is greater than a threshold. This occurs commonly in small data sets when an event is rare, and this is also the case here. To cope with this, we use the Bayesian generalized linear model approach (bayesglm function) implemented in the library arm in the statistical software R version 3.1.1 [62,63]. We used the default procedure settings, i.e., an independent Cauchy prior distribution for the model coefficients.

2.6. Multimodel Inference

Because of high uncertainty in model selection (several candidate models or even several dozens of candidate models may describe data similarly well), we used the multimodel inference method as a basis for investigating explanatory and predictive performance of the considered model variables. In this approach, inference is based not on one, but several models in the candidate set [64,65,66]. In particular, the Akaike information criterion (AIC) is calculated for each model in the candidate set (we actually calculate AICc, since our sample sizes are relatively small) and these models are ranked according to their Akaike weights, defined for a model i as:
w i = exp ( 1 2 Δ AIC ( model  i ) ) k = 1 M exp ( 1 2 Δ AIC ( model  k ) )
In this formula, ΔAIC is the difference between AIC of the model with the minimum AIC value and AIC of each respective model. The Akaike weight wi expresses the likelihood that the model i is the best approximating model [66]. We note that AIC (and any of its variants) provides only a relative measure of how good a model is, given the data and the candidate model set.
The Akaike weights allow estimating relative importance of every predictor variable in the full model by summing the Akaike weights for each model in which an individual variable occurs [64,66]. We calculated this for all predictor variables, including with and without the pre-disturbance Wetness and Brightness slopes, allowing assessment of importance of the remote sensing predictors relative to the environmental variables. Similar to the Akaike weight, this relative importance (or predictor weight) can be interpreted as equivalent to a likelihood that the respective variable is a part of the best model [64,66].
We conducted the above-described procedure for four nested model groups. The first group consists of all M candidate models (full model set). The second group consists of all candidate models in which the predictor variable Brightness slope is absent and the only candidate pre-disturbance predictor is the Wetness slope (Wetness slope model set). Similarly, the third group consists of all candidate models in which the Wetness slope is absent and the only candidate pre-disturbance predictor is Brightness slope (Brightness slope model set). Our final set of models did not contain any pre-disturbance variable (plain model set). We used these four model sets to assess predictive importance of the Wetness slope and Brightness slope variables. All calculations involving multimodel inference were performed using the R package MuMIn [67].

2.7. Predictions

From each of the four model sets (the full, Wetness slope, Brightness slope and plain model sets), candidate models with the highest Akaike weights (with summed Akaike weight of at least 0.95) were selected. These were then used to predict probability of bark beetle attack in the years different from those for which the models were fitted. This was done to assess whether models fitted for a specific year in the outbreak cycle could be used to predict dynamics of the bark beetle attack in its other part of the cycle or how robust the fitted models are across the outbreak. These “best” model sets contained tens of candidate models, as opposed to thousands of models in the original model sets. We used models fitted for the year F (referred to as F-models) to calculate probability of bark beetle attack in the year P (referred to as P-predictions) as follows. Any F-model was used individually to calculate P-predictions, using data from the year P as its “validation” data. These individual P-predictions were then averaged, weighting specific P-predictions by the respective F-model Akaike weights [64,66].
Binary logistic regression models allow modeling binary responses (presence/absence of infestation) by linear regression via a logistic link function. As such, they can predict the probability of infestation, which can then be classified into two discrete classes choosing a probability threshold [15]. Receiver Operating Characteristics (ROC) curves are a way of evaluating such binary models, since they are threshold-independent [68,69]. The area under the ROC curve (AUC) is often taken as a typical performance measure for a binary response variable since it provides a single measure of classification accuracy [68,69]. We quantified a predictive power of our model-averaged predictions using the ROC curves and AUC values. In particular, to assess an impact of considering the pre-disturbance spectral trajectories for making predictions, we statistically tested the null hypothesis that relevant pairs of the ROC curves (e.g., when both pre-disturbance variables are used vs. only one such variable is used) are statistically insignificant. We used the R package pROC [70] for calculating the ROC curves, AUC values, and statistical differences between any two ROC curves.

3. Results

3.1. Pre-Disturbance Forest Condition

Both pre-disturbance predictors Wetness slope and Brightness slope show a similar spatial pattern (Figure 3). The spatial variability of both of these predictors is quite high and the slope values vary from negative to positive, with a clear spatial autocorrelation. Except for the highest Wetness slope and Brightness slope in affected stands at the beginning of the outbreak relative to latter years of infestation [40], the most obvious observation here is that the slopes of both indices vary most near non-forested areas (Figure 3 and Figure 4). While values of the Wetness slope decrease near clear-cuts and locally increase near natural non-forested areas, the Brightness slope values demonstrate quite opposite behavior (Figure 3 and Figure 4). Its values increase near clear-cuts and decrease near natural non-forested areas. Still, the correlation between Wetness and Brightness slopes is not that high (r = −0.54). The linear regression of Wetness from 1984 to 1994 has mean values of coefficient of determination for pixels with affected forest higher (R2 = 0.35; SD 0.32) than pixels with unaffected forest until 1994 (R2 = 0.17; SD 0.33). Similar results show the linear regression of Brightness from 1984 to 1994 for affected forest, which has higher coefficient of determination (R2 = 27; SD 0.37) than unaffected (R2 = 0.14; SD 0.38). More detailed information about the linear fitting of pre-disturbance trajectories in our study area is given in [40].

3.2. Predictors of Bark Beetle Infestation

Statistical modelling based on the all-subsets approach and multimodel inference suggests that when both pre-disturbance predictors are taken into account, Wetness slope is an important variable to be considered for understanding the bark beetle infestation patterns (Table 3). Conversely, the likelihood that Brightness slope is a part of the best model is often as low as 0.3 (Table 3). In addition, whereas the high importance weights of Wetness slope are relatively stable over the years (with the exception of 2000), the Brightness slope weights show an increase in the years 1994–1996 and then a steady decline to lower values (≥0.28, see Table 3). The weights of other environmental predictors show relatively high variability over the years. However, similar to the pre-disturbance predictors, we recognize two basic trend-groups (Table 3). One group of predictors has a relatively high weight in the initial years of bark beetle infestation, while gradually losing weight towards the last modelled year. This group includes the edaphic category, which has high weight values until 1998 (with the exception of 1995); heat load index (HLI) and distance to wet or aquatic conditions (Geol-wet) show similar trends. The second group of predictors (TWI, AGE, and Dg-natur) reach the highest importance weights in the middle of bark beetle outbreak (1996–1997) when the largest forest areas are declined in our study area [45]. A special case of the opposite trend is represented by distance to clear-cut (D-clear-cut), which has the highest weights at the beginning and the end of the outbreak. The remaining predictors show somewhat fluctuating weights over the years.
To examine how explanatory power of models that include our two considered remote sensing variable changes relative to models in which one of these variables is excluded, we repeated the above analysis considering all environmental predictors but only with either Wetness or Brightness slope. With just Wetness slope (i.e., without Brightness slope), its importance weight and those of all environmental variables remained virtually unchanged (compare Table 3 and Table A1). On the contrary, with just Brightness slope (i.e., without Wetness slope), its importance weight increased in several years (1994, 1996), but decreased in other ones (1995, 1997; compare Table 3 and Table A2). These changes in the Brightness slope importance weights were accompanied by changes in the HLI importance weights. In particular, the HLI importance weight generally increased, with the highest increase in the years 1996 to 1998 being in the range of 10% to 20%. Similar importance weight changes were observed when the results based on the Wetness model set were compared with those based on the plain model set (high importance of Wetness slope; Table A3), and when the results based on the Brightness model set were compared with those based on the plain model set (low to medium importance of Brightness slope; Table A3).

3.3. Model Predictions

We assessed the predictive power of our fitted models by calculating the ROC curves and the corresponding AUC values for all years except the training year of the model. Regardless of the modelling scenario used, the trend is quite similar. Table 4 shows the results of models including both Wetness and Brightness slope. Generally, the predictive power of most fitted models decreases as time progresses. Further, models fitted at the beginning of bark beetle outbreak have higher predictive power than models fitted after outbreak culmination. The models describing an initial phase of bark beetle outbreaks thus appear more reliable for conducting near-future predictions. The higher predictive power of a few models for the year 2000 breaks this rule, but at the end of outbreak, reliability of models somewhat decreases because of low sample size.
Inclusion of the pre-disturbance indices generally increases model predictive power (Table 5). Although the improvement may seem marginal, it is statistically significant in many cases (Table A4). Interestingly, an increase is generally higher with backward prediction (Table 5, values below the diagonal) than forward one (Table 5, values above the diagonal). Differences in Table 5 would not change much if we replace the plain model set with the Brightness model set. On the other hand, if we consider differences between AUC values corresponding to the full model set and the Wetness model set, they are all virtually zero and statistically insignificant. This also suggests that with respect to the model predictive power, the Wetness slope plays a more important role than the Brightness slope.

4. Discussion

4.1. Spectral Indices Reflecting the Forest Health Change

Our results clearly demonstrate that long-term changes of TC Wetness and Brightness are related to the probability of bark beetle infestation during outbreak. This corroborates our hypothesis that bark beetles preferably affect the weakened forest stands, and that these changes in forest health or tree vigor are detectable in spectral signal change. In general the nature of spectral signal change is connected with the tree defoliation [37], chlorophyll fluorescence change and/or water content decrease in needles due to water stress [71]. This spectral change is well detected in NIR and SWIR spectral region [11]. In particular the SWIR-based indices (NDMI, Wetness, Brightness, NBR) are often used to assess change in forest health due to acid rain [72], defoliation following insect attack (e.g., [51]), or forest mortality (e.g., [73]). For this study, we selected Wetness and Brightness indices. These two indices show almost opposite behavior in reaction on forest decline [45]. In contrast, some authors have reported success using NIR-based indices. For example, Senf et al. [37] found stronger association of spruce budworm disturbance with Greenness than two other TC indices. However, the same study describes the highest sensitivity of Wetness and Brightness to mountain pine beetle disturbance. We removed Greenness early during our statistical analysis due to correlation and collinearity. Moreover, previous results from the same study area [45] support this removal, because they showed the inconsistent response of Greenness to the forest decline, when the magnitude of change was lower and the values changed with delay after forest decline. A possible reason could be the variability in understory composition; a dying canopy may uncover either soil and litter, or a green shrub/herbaceous layer. The sensitivity in spectral response to these contrasting cover types may be greater in Greenness than Wetness or Brightness and lead to inconsistent spectral-temporal trends associated with infestation-related tree mortality.
Further understanding of ecological processes underlying changes in Wetness and Brightness can be inferred from index trends we observed near forest edges. The decrease of Wetness and increase of Brightness and/or high variability of both predictors near clear-cuts may be related to water stress due to higher amount of solar radiation, particularly in southern and western forest edges [20,54]. Similar results were reported by Dantas de Paula et al. [74] in a Landsat-based study, where the tree cover variation increases towards the forest edge. The authors suggested the underlying causes for this variability to be degradation debt, hyperdynamism and environmental factors (topography, soil conditions etc.). In contrast to clear-cuts, the natural forest edges do not show obvious change of variability in the pre-disturbance indices. Nevertheless, some of these edges have higher Wetness slope values than do areas inside the forest. This effect is pronounced at the margin between forest and meadows (e.g., Lusen valley), where the high canopy closure protects the soil against drying, or denser canopies may serve as a physical defense to bark beetle infestation.

4.2. Pre-Disturbance Forest Conditions

If the trend in both indices indicates changes in forest health status or stand vigor, we should be able to describe their main causes. One of the strongest factors responsible for forest decline in Central Europe and North America in the second half of the last century was anthropogenic air pollution causing acid deposition. Even though the Šumava Mountains do not belong to the most heavily polluted areas in Czech Republic (e.g., the “black triangle” in northern Czech Republic with large historical forest decline), there is evidence from ground observations of forest decline due to air pollution in the last century [75,76]. Other retrospective methods suggest the consequences of acid deposition on forest health can be still recognized. Šantrůčková et al. [76] found a proportional relation between Δ 13C in tree rings and soil acidification (pH decrease). Polák et al. [77] describe tree defoliation followed by crown regeneration using secondary shots that were detectable even after 30 years. It is therefore likely that multi-temporal imagery captures this forest health decline due to acid deposition, which is in concordance with results reported by Rock et al. [72]. Another possible long-term stress factor in the Šumava Mountains is stand damage due to repeated windstorms or snow [7]. The repeated windstorms may damage the tree roots and cause tree vigor decrease. Drought is another possible long-term factor [28], and its role probably increases during bark beetle infestation.
All the aforementioned factors may reflect the spatial variability resulting from climate and topography differences from local to regional scales. Even though some studies mention no significant relation between forest health and bark beetle attack (e.g., [78]), Bytnerowicz et al. [79] suggested that long-term elevated levels of atmospheric N and S depositions and elevated O3 predispose trees to insect attacks and other stresses.

4.3. Spectral Trajectories

The use of linear regression for Wetness and Brightness fit impose simple trends over time. This assumption is consistent with results of Hais et al. [40], who reported lower Wetness slope values in the first years of the bark beetle outbreak in comparison to latter years, which were more stable. This suggests that the initial bark beetle attack occurred in the most vulnerable stands that show long-term Wetness/Brightness change before the disturbance. Nevertheless, the resulting spectral trajectory may be influenced by other factors, which can act antagonistically. While long-term stress factors decrease Wetness values (or increase Brightness) over time, canopy growth may support an opposite trend, in particular for fast-growing tree species. For example, partial forest defoliation may result in small spectral change, which is quickly followed by spectral recovery after disturbance [37,38]. Similarly, forest thinning is reflected in Wetness change [16] when the values first decrease due to lower canopy closure (and higher role of ground cover in the spectral signal), but after several years, the values increase due to faster tree growth [16]. Users should therefore be careful when using these indices in managed forests. Additionally, recent forest history often includes high magnitude changes that are not always linear in their spectral-temporal response. Because of this complexity, the non-linear fit of spectral trajectory is often needed to describe forest history for a given pixel [35]. Nevertheless, much easier interpretable linear fits are valuable for description of many long-term processes with unaltered trend like gradual systematic change in natural vegetation communities [80]. Similarly we used the linear fit to describe pre-disturbance forest condition assuming it mostly follows simple trend. However, to avoid complexity in the course of spectral-temporal response we use just forest pixels without any type of disturbance during assessed pre-disturbance time.

4.4. Risk of Bark Beetle Infestation and Model Predictions

Among all predictors used in our study, Wetness slope provided the greatest improvement of bark beetle infestation risk modelling, while Brightness slope was a significant predictor only for a few years. Moreover, Wetness slope and Brightness slope are not correlated with the environmental predictors, and thus are able to explain additional observed variability. The Landsat spatial resolution (30 m) provides a good dataset to describe the spatial variability on a stand level. Relatively low weights of importance for some environmental variables may be caused by the small dataset, especially when the spatial resolution of the data is relatively coarse and/or the predictor is categorical. This may be the case for very low weights of broad edaphic categories in models for years 1995 and 2000 (lowest number of cases), which are high in other years. This limitation has lower importance for continuous variables like pre-disturbance indices. This supports the potential of the Landsat data for local to global-level modelling in general. In agreement with Lausch et al. [21], we conclude that the importance of individual predictors changes during the bark beetle infestation. The predictors with the highest model weights (edaphic category, Geol-wet, HLI) in the first years of the bark beetle outbreak document the selective infestation of the most susceptible forest stands. The selective infestation is critical until bark beetle population size reaches the threshold where tree defense is overwhelmed [78]. After successful selective infestation, bark beetle propagation produces population sizes so great that selective behavior is an unnecessary strategy to seize trees, which shifts increased importance to other predictors. This was demonstrated in 1997, when the largest forest area was affected, and tree-related predictors had the highest weights (age, density, DG-natur). The trees were probably more susceptible individually, reflecting stand structure and not the environment predictors as in the beginning of the bark beetle outbreak. The decrease in almost all predictor weights further into the future, shows that the high bark beetle density led to random-like infestation, where the beetles were able to overcome any tree’s defense. The predictor D-clear-cut has high importance in the first and last years of bark beetle outbreak, which might be connected with low bark beetle numbers in both periods and need for selective behavior.
Predictive power of our fitted models (including all predictors) calculated using ROC curves and corresponding AUC values was quite similar. The high AUC values when data were used for fitting and also for prediction (0.76–0.91) documented reasonably good fit to the data, when assessed using AUCs. The decrease of predictive power in fitted models over the years documents not only the general decrease of fits to the distant years, but it also supports the explanation of predictor importance changes during a bark beetle outbreak [21]. Moreover, the slower decrease of predictive power (higher predictive power for distant years) for models fitted at the beginning of the bark beetle outbreak nicely corresponds with the hypothesis that at the beginning of the bark beetle outbreak beetles are more selective, while later, when beetle densities rise and/or the most susceptible stands are already destroyed, beetles need to be more opportunistic.

5. Conclusions

Using pre-disturbance slopes of Tasseled Cap components Wetness and Brightness, spectral trajectories improved our understanding and prediction of bark beetle infestation. Linear regression slopes of these two variables calculated for the time period before the disturbance significantly increase the predictive power of bark beetle infestation models from year to year and even for longer predictions. Moreover, these predictors were not correlated with environmental predictors and therefore are able to explain additional observed variability. This also suggests that long-term gradual decrease of forest vitality significantly increases forest vulnerability to bark beetle infestation. When comparing the importance weights of both pre-disturbance predictors, the Wetness slope weights are high and relatively stable over the years (with the exception of 2000), while the weights of Brightness slopes show an increase only in the middle of the bark beetle outbreak (1996) and steadily decline to low values. We distinguished two components of relationship between these two indices and tree susceptibility to bark beetle infestation. First, they reflect long-term stress factors (acid deposition, drought, snow, windstorms). The second component is their change near forest clear-cuts edges. This probably relates to higher solar radiation at open edges or other disturbances caused during the clear-cut.
We have also identified two groups of environmental factors important for prediction of bark beetle infestation risk. The first group has high importance weights in the initial years of the outbreak (edaphic category, distance to aquatic conditions, heat load index), which might be connected with high pest selectivity for the most susceptible stands. The second group of tree-related predictors has the highest weights in the middle of the outbreak, when the largest forest area was affected, and the pest numbers were highest (age of dominant tree layer, tree density of dominant tree layer, Distance to the natural non-forested area). The first can be used to identify potentially vulnerable areas, and the second for predicting bark beetle spread after the initial phase of infestation.
But forest ecosystems also face long-term, gradual changes due to hidden stress factors that contribute to their vulnerability and susceptibility to bark beetle infestation. And this is where non-correlated pre-disturbance spectral trajectories, based on Landsat data with its vast and accessible archive, can really help.

Acknowledgments

The work was supported by COST LD 15158 and was part of long-term research development project No. RVO 67985939. We thank Zdeňka Konopová for help with vectorization of land-cover units from aerial photographs, Karel Matějka for data support, Stanislav Grill for discussions focused on insect forest disturbance and four anonymous reviewers for helpful comments.

Author Contributions

Martin Hais conceived the general idea of the paper, processed satellite data, made GIS analyses and wrote the manuscript. Jan Wild authored ecological interpretation of modelling results. He wrote the manuscript and analyzed the GIS data. Luděk Berec performed the main statistical analyses and wrote the manuscript. Josef Brůna processed the GIS data and wrote the manuscript. Robert Kennedy formed the idea of pre-disturbance spectral trajectories and supervised the manuscript writing. Justin Braaten made the radiometric normalization of the Landsat imagery, wrote and corrected the manuscript. Zdeněk Brož performed the aerial images rectification, delineated the trees affected by bark beetles and calculated some of the environmental variables.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Appendix A

A1. Predictor Importance Weights

Table A1. Predictor importance weights based on multimodel inference when only the Wetness slope is used with the environmental predictors. The darkest/lightest green show the highest/lowest predictor importance weights.
Table A1. Predictor importance weights based on multimodel inference when only the Wetness slope is used with the environmental predictors. The darkest/lightest green show the highest/lowest predictor importance weights.
YearWet-slDir-SourceD-Clear-CutD-Natur-EdgDg-NaturAGEDensityEdaphicGeol-WetHLITWILegend
19940.991.001.000.960.520.650.361.000.981.000.62≤0.2
19950.680.090.330.520.560.360.480.010.390.840.650.21–0.40
19961.000.070.371.001.000.670.301.001.000.701.000.41–0.60
19971.000.950.480.291.001.001.001.000.540.741.000.61–0.80
19981.001.000.261.000.300.290.370.990.290.550.610.81–1.00
19990.990.450.680.340.860.310.940.470.480.290.30
20000.260.010.950.320.380.260.320.010.370.670.44
Table A2. Predictor importance weights based on multimodel inference when only the Brightness slope is used with the environmental predictors. The darkest/lightest green show the highest/lowers predictor importance weights.
Table A2. Predictor importance weights based on multimodel inference when only the Brightness slope is used with the environmental predictors. The darkest/lightest green show the highest/lowers predictor importance weights.
YearBright-slDir-SourceD-Clear-CutD-Natur-EdgDg-NaturAGEDensityEdaphicGeol-WetHLITWILegend
19940.901.001.000.960.520.710.351.000.981.000.62≤0.2
19950.270.070.290.550.590.360.490.010.380.860.650.21–0.40
19961.000.140.391.001.000.670.321.001.000.821.000.41–0.60
19970.330.990.360.291.001.001.001.000.690.941.000.61–0.80
19980.361.000.261.000.270.280.330.990.350.710.590.81–1.00
19990.290.220.740.340.840.310.930.670.460.320.28
20000.280.010.950.320.380.270.320.010.370.670.44
Table A3. Predictor importance weights based on multimodel inference when no pre-disturbance index is used as a predictor. The darkest/lightest green show the highest/lowest predictor importance weights.
Table A3. Predictor importance weights based on multimodel inference when no pre-disturbance index is used as a predictor. The darkest/lightest green show the highest/lowest predictor importance weights.
YearDir-SourceD-Clear-CutD-Natur-EdgDg-NaturAGEDensityEdaphicGeol-WetHLITWILegend
19941.001.000.980.560.770.331.000.971.000.61≤0.2
19950.070.290.550.590.360.490.010.380.860.650.21–0.40
19960.160.681.001.000.600.331.001.000.881.000.41–0.60
19970.990.350.291.001.001.001.000.690.941.000.61–0.80
19981.000.261.000.270.280.320.990.350.720.590.81–1.00
19990.220.740.340.840.310.920.680.460.330.28
20000.010.950.320.380.270.320.010.370.670.44

A2. Predictive Power of the Remote Sensing Variables

Table A4. Results of comparison tests of two ROC curves, one corresponding to the full model set (with both Wetness and Brightness slopes) and the other to the plain model set (with no remote sensing variable). p-values are shown. The ROC curve differs significantly for most comparisons.
Table A4. Results of comparison tests of two ROC curves, one corresponding to the full model set (with both Wetness and Brightness slopes) and the other to the plain model set (with no remote sensing variable). p-values are shown. The ROC curve differs significantly for most comparisons.
1994199519961997199819992000
1994-0.090.000.000.000.000.00
19950.00-0.000.000.000.810.93
19960.000.31-0.010.000.860.03
19970.000.010.00-0.000.000.04
19980.000.010.000.00-0.040.46
19990.000.090.000.000.01-0.41
20000.000.590.020.420.490.36-

References

  1. Frelich, L.E. Forest Dynamics and Disturbance Regimes; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar]
  2. Waring, R.H.; Running, S.W. Forest Ecosystems: Analysis at Multiple Scales, 3rd ed.; Elsevier: San Diego, CA, USA, 2007. [Google Scholar]
  3. Seidl, R.; Schelhaas, M.-J.; Lexer, M.J. Unraveling the drivers of intensifying forest disturbance regimes in Europe. Glob. Chang. Biol. 2011, 17, 2842–2852. [Google Scholar] [CrossRef]
  4. Schelhaas, M.J.; Nabuurs, G.J.; Schuck, A. Natural disturbances in the European forests in the 19th and 20th centuries. Glob. Chang. Biol. 2003, 9, 1620–1633. [Google Scholar] [CrossRef]
  5. Seidl, R.; Schelhaas, M.-J.; Rammer, W.; Verkerk, P.J. Increasing forest disturbances in Europe and their impact on carbon storage. Nat. Clim. Chang. 2014, 4, 806–810. [Google Scholar] [CrossRef] [PubMed]
  6. Wild, J.; Kopecký, M.; Svoboda, M.; Zenáhlíková, J.; Edwards-Jonášová, M.; Herben, T. Spatial patterns with memory: Tree regeneration after stand-replacing disturbance in Picea abies mountain forests. J. Veg. Sci. 2014, 25, 1327–1340. [Google Scholar] [CrossRef]
  7. Čada, V.; Morrissey, R.C.; Michalová, Z.; Bače, R.; Janda, P.; Svoboda, M. Frequent severe natural disturbances and non-equilibrium landscape dynamics shaped the mountain spruce forest in central Europe. For. Ecol. Manag. 2016, 363, 169–178. [Google Scholar] [CrossRef]
  8. Svoboda, M.; Janda, P.; Nagel, T.A.; Fraver, S.; Rejzek, J.; Bače, R. Disturbance history of an old-growth sub-alpine Picea abies stand in the Bohemian Forest, Czech Republic. J. Veg. Sci. 2012, 23, 86–97. [Google Scholar] [CrossRef]
  9. Seidl, R.; Rammer, W.; Jäger, D.; Lexer, M.J. Impact of bark beetle (Ips typographus L.) disturbance on timber production and carbon sequestration in different management strategies under climate change. For. Ecol. Manag. 2008, 256, 209–220. [Google Scholar] [CrossRef]
  10. Lambert, N.J.; Ardo, J.; Rock, B.N.; Vogelmann, J.E. Spectral characterization and regression-based classification of forest damage in Norway spruce stands in the Czech Republic using Landsat Thematic Mapper data. Int. J. Remote Sens. 1995, 16, 1261–1287. [Google Scholar] [CrossRef]
  11. Rullan-Silva, C.D.; Olthoff, A.E.; Delgado de la Mata, J.A.; Pajares-Alonso, J.A. Remote monitoring of forest insect defoliation. A review. For. Syst. 2013, 22, 377–391. [Google Scholar] [CrossRef]
  12. Kauth, R.J.; Thomas, G.S. The tasselled cap—A graphic description of the spectral-temporal development of agricultural crops as seen by Landsat. In Proceedings of the Symposium Machine Processing of Remotely Sensed Data, West Lafayette, IN, USA, 29 June–1 July 1976; pp. 41–51.
  13. Crist, E.P.; Cicone, R.C. Application of the Tasseled Cap concept to simulated thematic mapper data. Photogramm. Eng. Remote Sens. 1984, 50, 343–352. [Google Scholar]
  14. Skakun, R.S.; Wulder, M.A.; Franklin, S.E. Sensitivity of the thematic mapper enhanced wetness difference index to detect mountain pine beetle red-attack damage. Remote Sens. Environ. 2003, 86, 433–443. [Google Scholar] [CrossRef]
  15. Wulder, M.A.; White, J.C.; Bentz, B.; Alvarez, M.F.; Coops, N.C. Estimating the probability of mountain pine beetle red-attack damage. Remote Sens. Environ. 2006, 101, 150–166. [Google Scholar] [CrossRef]
  16. Coops, N.C.; Waring, R.H.; Wulder, M.A.; White, J.C. Prediction and assessment of bark beetle-induced mortality of lodgepole pine using estimates of stand vigor derived from remotely sensed data. Remote Sens. Environ. 2009, 113, 1058–1066. [Google Scholar] [CrossRef]
  17. Lausch, A.; Heurich, M.; Gordalla, D.; Dobner, H.-J.; Gwillym-Margianto, S.; Salbach, C. Forecasting potential bark beetle outbreaks based on spruce forest vitality using hyperspectral remote-sensing techniques at different scales. For. Ecol. Manag. 2013, 308, 76–89. [Google Scholar] [CrossRef]
  18. Overbeck, M.; Schmidt, M. Modelling infestation risk of Norway spruce by Ips typographus (L.) in the Lower Saxon Harz Mountains (Germany). For. Ecol. Manag. 2012, 266, 115–125. [Google Scholar] [CrossRef]
  19. Netherer, S.; Nopp-Mayr, U. Predisposition assessment systems (PAS) as supportive tools in forest management—Rating of site and stand-related hazards of bark beetle infestation in the High Tatra Mountains as an example for system application and verification. For. Ecol. Manag. 2005, 207, 99–107. [Google Scholar] [CrossRef]
  20. Dutilleul, P.; Nef, L.; Frigon, D. Assessment of site characteristics as predictors of the vulnerability of Norway spruce (Picea abies Karst.) stands to attack by Ips typographus L. (Col., Scolytidae). J. Appl. Entomol. 2000, 124, 1–5. [Google Scholar] [CrossRef]
  21. Lausch, A.; Fahse, L.; Heurich, M. Factors affecting the spatio-temporal dispersion of Ips typographus (L.) in Bavarian Forest National Park: A long-term quantitative landscape-level analysis. For. Ecol. Manag. 2011, 261, 233–245. [Google Scholar] [CrossRef]
  22. Lausch, A.; Heurich, M.; Fahse, L. Spatio-temporal infestation patterns of Ips typographus (L.) in the Bavarian Forest National Park, Germany. Ecol. Indic. 2013, 31, 73–81. [Google Scholar] [CrossRef]
  23. Wulder, M.A.; Franklin, S.E. Understanding Forest Disturbance and Spatial Pattern: Remote Sensing and GIS Approaches; CRC Press, Taylor & Francis Group: Boca Raton, FL, USA, 2007. [Google Scholar]
  24. Bjørnstad, O.N.; Økland, B. Synchrony and geographical variation of the spruce bark beetle (Ips typographus) during a non-epidemic period. Popul. Ecol. 2003, 45, 213–219. [Google Scholar] [CrossRef]
  25. Brůna, J.; Wild, J.; Svoboda, M.; Heurich, M.; Müllerová, J. Impacts and underlying factors of landscape-scale, historical disturbance of mountain forest identified using archival documents. For. Ecol. Manag. 2013, 305, 294–306. [Google Scholar] [CrossRef]
  26. Aakala, T.; Kuuluvainen, T.; Wallenius, T.; Kauhanen, H. Tree mortality episodes in the intact Picea abies-dominated taiga in the Arkhangelsk region of northern European Russia. J. Veg. Sci. 2011, 22, 322–333. [Google Scholar] [CrossRef]
  27. Schroeder, L.M.; Weslien, J.; Lindelöw, Å.; Lindhe, A. Attacks by bark- and wood-boring Coleoptera on mechanically created high stumps of Norway spruce in the two years following cutting. For. Ecol. Manag. 1999, 123, 21–30. [Google Scholar] [CrossRef]
  28. Seidl, R.; Müller, J.; Hothorn, T.; Bässler, C.; Heurich, M.; Kautz, M. Small beetle, large-scale drivers: How regional and landscape factors affect outbreaks of the European spruce bark beetle. J. Appl. Ecol. 2015, 53, 530–540. [Google Scholar] [CrossRef] [PubMed]
  29. Wang, J.; Sammis, T.W.; Gutschick, V.P.; Gebremichael, M.; Dennis, S.O.; Harrison, R.E. Review of satellite remote sensing use in forest health studies. Open Geogr. J. 2010, 3, 28–42. [Google Scholar] [CrossRef]
  30. Wulder, M.A.; Masek, J.G.; Cohen, W.B.; Loveland, T.R.; Woodcock, C.E. Opening the archive: How free data has enabled the science and monitoring promise of Landsat. Remote Sens. Environ. 2012, 122, 2–10. [Google Scholar] [CrossRef]
  31. Kennedy, R.E.; Yang, Z.; Cohen, W.B. Detecting trends in forest disturbance and recovery using yearly Landsat time series: LandTrendr—Temporal segmentation algorithms. Remote Sens. Environ. 2010, 114, 2897–2910. [Google Scholar] [CrossRef]
  32. Cohen, W.B.; Yang, Z.; Kennedy, R. Detecting trends in forest disturbance and recovery using yearly Landsat time series: TimeSync—Tools for calibration and validation. Remote Sens. Environ. 2010, 114, 2911–2924. [Google Scholar] [CrossRef]
  33. Kennedy, R.E.; Andréfouët, S.; Cohen, W.B.; Gómez, C.; Griffiths, P.; Hais, M.; Healey, S.P.; Helmer, E.H.; Hostert, P.; Lyons, M.B.; et al. Bringing an ecological view of change to Landsat-based remote sensing. Front. Ecol. Environ. 2014, 12, 339–346. [Google Scholar] [CrossRef] [Green Version]
  34. Pflugmacher, D.; Cohen, W.B.; Kennedy, R.E. Using Landsat-derived disturbance history (1972–2010) to predict current forest structure. Remote Sens. Environ. 2012, 122, 146–165. [Google Scholar] [CrossRef]
  35. Kennedy, R.E.; Yang, Z.; Cohen, W.B.; Pfaff, E.; Braaten, J.; Nelson, P. Spatial and temporal patterns of forest disturbance and regrowth within the area of the Northwest Forest Plan. Remote Sens. Environ. 2012, 122, 117–133. [Google Scholar] [CrossRef]
  36. Meddens, A.J.H.; Hicke, J.A.; Vierling, L.A.; Hudak, A.T. Evaluating methods to detect bark beetle-caused tree mortality using single-date and multi-date Landsat imagery. Remote Sens. Environ. 2013, 132, 49–58. [Google Scholar] [CrossRef]
  37. Senf, C.; Pflugmacher, D.; Wulder, M.A.; Hostert, P. Characterizing spectral–temporal patterns of defoliator and bark beetle disturbances using Landsat time series. Remote Sens. Environ. 2015, 170, 166–177. [Google Scholar] [CrossRef]
  38. Meigs, G.W.; Kennedy, R.E.; Cohen, W.B. A Landsat time series approach to characterize bark beetle and defoliator impacts on tree mortality and surface fuels in conifer forests. Remote Sens. Environ. 2011, 115, 3707–3718. [Google Scholar] [CrossRef]
  39. Cohen, W.B.; Yang, Z.; Stehman, S.V.; Schroeder, T.A.; Bell, D.M.; Masek, J.G.; Huang, C.; Meigs, G.W. Forest disturbance across the conterminous United States from 1985 to 2012: The emerging dominance of forest decline. For. Ecol. Manag. 2016, 360, 242–252. [Google Scholar] [CrossRef]
  40. Hais, M.; Kennedy, R.; Braaten, J.; Pflugmacher, D. Landsat spectral trajectories may aid in detecting precursors to spruce bark beetle outbreak. Unpublished Manuscript. 2016. [Google Scholar]
  41. Neuhäuslová, Z. (Ed.) The Map of Potential Natural Vegetation of the Šumava National Park; Explanatory Text; Suppl. 1/2001; Silva Gabreta: Vimperk, Czech Republic, 2001.
  42. Fischer, H.S.; Winter, S.; Lohberger, E.; Jehl, H.; Fischer, A. Improving transboundary maps of potential natural vegetation using statistical modeling based on environmental predictors. Folia Geobot. 2013, 48, 115–135. [Google Scholar] [CrossRef]
  43. Heurich, M.; Reinelt, A.; Fahse, L. Die buchdruckermassenvermehrung im nationalpark bayerischer wald. In Waldentwicklung im Bergwald nach Wind, und Borkenkäferbefall; Bayer, Staatsforstverwaltung Wiss. Reihe: Grafenau, Germany, 2001; Volume 14, pp. 9–48. [Google Scholar]
  44. Skuhravý, V. Lýkožrout Smrkový Ips typographus (L.) a Jeho Kalamity; Agrospoj: Praha, Czech Republic, 2002. [Google Scholar]
  45. Hais, M.; Jonášová, M.; Langhammer, J.; Kučera, T. Comparison of two types of forest disturbance using multitemporal Landsat TM/ETM+ imagery and field vegetation data. Remote Sens. Environ. 2009, 113, 835–845. [Google Scholar] [CrossRef]
  46. Nováková, M.H.; Edwards-Jonášová, M. Restoration of Central-European mountain Norway spruce forest 15 years after natural and anthropogenic disturbance. For. Ecol. Manag. 2015, 344, 120–130. [Google Scholar] [CrossRef]
  47. Zeppenfeld, T.; Svoboda, M.; DeRose, R.J.; Heurich, M.; Müller, J.; Čížková, P.; Starý, M.; Bače, R.; Donato, D.C. Response of mountain Picea abies forests to stand-replacing bark beetle outbreaks: Neighbourhood effects lead to self-replacement. J. Appl. Ecol. 2015, 52, 1402–1411. [Google Scholar] [CrossRef]
  48. Canty, M.J.; Nielsen, A.A.; Schmidt, M. Automatic radiometric normalization of multitemporal satellite imagery. Remote Sens. Environ. 2004, 91, 441–451. [Google Scholar] [CrossRef] [Green Version]
  49. Schroeder, T.A.; Cohen, W.B.; Song, C.; Canty, M.J.; Yang, Z. Radiometric correction of multi-temporal Landsat data for characterization of early successional forest patterns in western Oregon. Remote Sens. Environ. 2006, 103, 16–26. [Google Scholar] [CrossRef]
  50. Cohen, W.B.; Spies, T.A.; Fiorella, M. Estimating the age and structure of forests in a multi-ownership landscape of western Oregon, U.S.A. Int. J. Remote Sens. 1995, 16, 721–746. [Google Scholar] [CrossRef]
  51. Bonneau, L.R.; Shields, K.S.; Civco, D.L. Using Satellite Images to Classify and Analyze the Health of Hemlock Forests Infested by the Hemlock Woolly Adelgid. Biol. Invasions 1999, 1, 255–267. [Google Scholar] [CrossRef]
  52. Wulder, M.A.; Skakun, R.S.; Kurz, W.A.; White, J.C. Estimating time since forest harvest using segmented Landsat ETM+ imagery. Remote Sens. Environ. 2004, 93, 179–187. [Google Scholar] [CrossRef]
  53. Jin, S.; Sader, S.A. Comparison of time series tasseled cap wetness and the normalized difference moisture index in detecting forest disturbances. Remote Sens. Environ. 2005, 94, 364–372. [Google Scholar] [CrossRef]
  54. Jakuš, R. A method for the protection of spruce stands againstIps typographus by the use of barriers of pheromone traps in north-eastern Slovakia. Anz. Schädlingskd Pflanzenschutz Umweltschutz 1998, 71, 152–158. [Google Scholar] [CrossRef]
  55. Matějka, K. Forest Altitudinal Zones with Dominant Norway Spruce in the Czech Republic. Available online: http://www.infodatasys.cz/public/Lesnik21_2014km.pdf (accessed on 25 May 2016).
  56. Randuška, D. Forest typology in Czechoslovakia. In Application of Vegetation Science to Forestry. Handbook of Vegetation Science; Jahn, G., Ed.; Dr. W. Junk Publication: The Hague, The Netherlands, 1982; Volume 12, pp. 147–178. [Google Scholar]
  57. Tüxen, R. Die Heutige Potentielle Natürliche Vegetation Als Gegenstand der Vegetationskartierung; Zentralstelle für Vegetationskartierung: Remagen, Germany, 1956. [Google Scholar]
  58. Viewegh, J.; Kusbach, A.; Mikeska, M. Czech forest ecosystem classification. J. For. Sci. 2003, 49, 85–93. [Google Scholar]
  59. McCune, B.; Keon, D. Equations for potential annual direct incident radiation and heat load. J. Veg. Sci. 2002, 13, 603–606. [Google Scholar] [CrossRef]
  60. Matějka, K. Hodnocení Přirozenosti Lesů NP Šumava (Evaluation of Naturalness of Forests in the Šumava National Park); Report for the Ministry of Environment; Ministry of Environment: Prague, Czech Republic, 2014.
  61. Zuur, A.F.; Ieno, E.N.; Elphick, C.S. A protocol for data exploration to avoid common statistical problems. Methods Ecol. Evol. 2010, 1, 3–14. [Google Scholar] [CrossRef]
  62. Gelman, A.; Su, Y.-S. Data Analysis Using Regression and Multilevel/Hierarchical Models; Cambridge University Press: Cambridge, UK, 2015. [Google Scholar]
  63. R Development Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2014. [Google Scholar]
  64. Burnham, K.P.; Anderson, D.R.; Huyvaert, K.P. AIC model selection and multimodel inference in behavioral ecology: Some background, observations, and comparisons. Behav. Ecol. Sociobiol. 2011, 65, 23–35. [Google Scholar] [CrossRef]
  65. Grueber, C.E.; Nakagawa, S.; Laws, R.J.; Jamieson, I.G. Multimodel inference in ecology and evolution: Challenges and solutions. J. Evol. Biol. 2011, 24, 699–711. [Google Scholar] [CrossRef] [PubMed]
  66. Symonds, M.R.E.; Moussalli, A. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behav. Ecol. Sociobiol. 2011, 65, 13–21. [Google Scholar] [CrossRef]
  67. Bartoń, K. MuMIn: Multi-Model Inference; R Package Version 1.15; Bartoń: Krakow, Poland, 2016. [Google Scholar]
  68. Fawcett, T. An introduction to ROC analysis. Pattern Recognit. Lett. 2006, 27, 861–874. [Google Scholar] [CrossRef]
  69. Fielding, A.H.; Bell, J.F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 1997, 24, 38–49. [Google Scholar] [CrossRef]
  70. Robin, X.; Turck, N.; Hainard, A.; Tiberti, N.; Lisacek, F.; Sanchez, J.-C.; Müller, M. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. 2011, 12. [Google Scholar] [CrossRef] [PubMed]
  71. Yoshida, Y.; Joiner, J.; Tucker, C.; Berry, J.; Lee, J.-E.; Walker, G.; Reichle, R.; Koster, R.; Lyapustin, A.; Wang, Y. The 2010 Russian drought impact on satellite measurements of solar-induced chlorophyll fluorescence: Insights from modeling and comparisons with parameters derived from satellite reflectances. Remote Sens. Environ. 2015, 166, 163–177. [Google Scholar] [CrossRef]
  72. Rock, B.N.; Vogelmann, J.E.; Williams, D.L.; Vogelmann, A.F.; Hoshizaki, T. Remote Detection of forest damage. Bioscience 1986, 36, 439–445. [Google Scholar] [CrossRef]
  73. Collins, J.B.; Woodcock, C.E. An assessment of several linear change detection techniques for mapping forest mortality using multitemporal Landsat TM data. Remote Sens. Environ. 1996, 56, 66–77. [Google Scholar] [CrossRef]
  74. Dantas de Paula, M.; Groeneveld, J.; Huth, A. The extent of edge effects in fragmented landscapes: Insights from satellite measurements of tree cover. Ecol. Indic. 2016, 69, 196–204. [Google Scholar] [CrossRef]
  75. Kopácek, J.; Veselý, J.; Stuchlík, E. Sulphur and nitrogen fluxes and budgets in the Bohemian Forest and Tatra Mountains during the Industrial Revolution (1850–2000). Hydrol. Earth Syst. Sci. 2001, 5, 391–406. [Google Scholar] [CrossRef]
  76. Šantrůčková, H.; Šantrůček, J.; Šetlík, J.; Svoboda, M.; Kopáček, J. Carbon isotopes in tree rings of Norway spruce exposed to atmospheric pollution. Environ. Sci. Technol. 2007, 41, 5778–5782. [Google Scholar] [CrossRef] [PubMed]
  77. Polák, T.; Cudlín, P.; Moravec, I.; Albrechtová, J. Macroscopic indicators for the retrospective assessment of Norway spruce crown response to stress in the Krkonoše Mountains. Trees 2007, 21, 23–35. [Google Scholar] [CrossRef]
  78. Wermelinger, B. Ecology and management of the spruce bark beetle Ips typographus—A review of recent research. For. Ecol. Manag. 2004, 202, 67–82. [Google Scholar] [CrossRef]
  79. Bytnerowicz, A.; Badea, O.; Barbu, I.; Fleischer, P.; Frączek, W.; Gancz, V.; Godzik, B.; Grodzińska, K.; Grodzki, W.; Karnosky, D.; et al. New international long-term ecological research on air pollution effects on the Carpathian Mountain forests, Central Europe. Environ. Int. 2003, 29, 367–376. [Google Scholar] [CrossRef]
  80. Vogelmann, J.E.; Xian, G.; Homer, C.; Tolk, B. Monitoring gradual ecosystem change using Landsat time series analyses: Case studies in selected forest and rangeland ecosystems. Remote Sens. Environ. 2012, 122, 92–105. [Google Scholar] [CrossRef]
Figure 1. Study area located in the Šumava Mountains, Czech Republic–Germany border region.
Figure 1. Study area located in the Šumava Mountains, Czech Republic–Germany border region.
Remotesensing 08 00687 g001
Figure 2. Spectral trajectory of Wetness TC component in a given pixel calculated as the slope of linear regression of Wetness for consecutive years from 1984 to one year before the year for which the model was calibrated. Slope was calculated for both infested and non-infested pixels; the infested one is shown as an example.
Figure 2. Spectral trajectory of Wetness TC component in a given pixel calculated as the slope of linear regression of Wetness for consecutive years from 1984 to one year before the year for which the model was calibrated. Slope was calculated for both infested and non-infested pixels; the infested one is shown as an example.
Remotesensing 08 00687 g002
Figure 3. Spatial distribution of Wetness slope and Brightness slope values. Similar to Figure 2, the pixel values were calculated as the slope of linear regression of Wetness/Brightness for consecutive years from 1984 to one year before the year for which the model was calibrated. The maps are composed from pixels combined among all evaluated years (1994–2000) according to the year of bark beetle infestation. The lowest Wetness slope values and highest Brightness slope values are mostly adjacent to clear-cuts, while the slope values in the vicinity of natural non-forested areas show opposite behavior, especially in the Lusen valley.
Figure 3. Spatial distribution of Wetness slope and Brightness slope values. Similar to Figure 2, the pixel values were calculated as the slope of linear regression of Wetness/Brightness for consecutive years from 1984 to one year before the year for which the model was calibrated. The maps are composed from pixels combined among all evaluated years (1994–2000) according to the year of bark beetle infestation. The lowest Wetness slope values and highest Brightness slope values are mostly adjacent to clear-cuts, while the slope values in the vicinity of natural non-forested areas show opposite behavior, especially in the Lusen valley.
Remotesensing 08 00687 g003
Figure 4. Wetness slope and Brightness slope values as they vary with distance to the nearest clear-cuts (A,B) and natural non-forested areas (C,D). The trend (blue line) of slope among distances is visualized using local polynomial regression (loess) ±1 standard error mean (grey area).
Figure 4. Wetness slope and Brightness slope values as they vary with distance to the nearest clear-cuts (A,B) and natural non-forested areas (C,D). The trend (blue line) of slope among distances is visualized using local polynomial regression (loess) ±1 standard error mean (grey area).
Remotesensing 08 00687 g004
Table 1. List of Landsat TM5 images (path 192, row 25).
Table 1. List of Landsat TM5 images (path 192, row 25).
Acquisition Date
3 August 1984
9 August 1986
11 July 1987
6 August 1988
5 July 1988
8 July 1989
4 August 1990
7 August 1991
9 August 1992
31 August 1994
14 July 1994
30 July 1994
1 July 1995
7 August 1997
10 August 1998
14 September 1999
Table 2. List of predictors for bark beetle infestation including pre-disturbance indices describing forest susceptibility and predictors of bark beetle dispersal.
Table 2. List of predictors for bark beetle infestation including pre-disturbance indices describing forest susceptibility and predictors of bark beetle dispersal.
Variable Name and DescriptionAbbreviationStateUnitsData Source
Forest structure
Distance to the nearest infested forest in the previous yearDist-sourcedynamic(m)map of land cover
Direction to the nearest infested forest in the previous yearDir-sourcedynamic8 classesmap of land cover
Distance to the nearest clear-cutD-clear-cutstatic(m)map of land cover
Distance to the natural non forested areaD-natur-edgstatic(m)map of land cover
Degree of naturalness of forest standDg-naturstatic6 classes[55,60]
Age of dominant tree layerAgestaticyearsFMI database
Tree density of dominant tree layerDensitystatictrees no./haFMI database
Pedology and geology
Edaphic categoryEdaphicstatic11 classesFMI database
Distance to wet or aquatic conditionsGeol-wetstatic(m)Geological map
Topography and climate
AltitudeDEMstatic(m)DMR 4G DEM
Heat load indexHLIstaticunitlessDMR 4G DEM
Area solar radiation for entire yearASR12staticWh/m2DMR 4G DEM
Topography wetness indexTWIstaticunitlessDMR 4G DEM
Pre-disturbance forest conditions
Slope of Brightness indexBright-sldynamicunitlessLandsat imagery
Slope of Wetness indexWet-sldynamicunitlessLandsat imagery
Slope of Greenness indexGreen-sldynamicunitlessLandsat imagery
Table 3. Predictor importance weights based on multimodel inference when both pre-disturbance indices are used as predictors. The darkest/lightest green show the highest/lowers predictor importance weights.
Table 3. Predictor importance weights based on multimodel inference when both pre-disturbance indices are used as predictors. The darkest/lightest green show the highest/lowers predictor importance weights.
YearWet-slBright-slDir-SourceD-Clear-CutD-natur-EdgeDg-NaturAGEDensityEdaphicGeol-WetHLITWILegend
19940.930.38110.960.530.650.3610.9810.62≤0.2
19950.810.570.090.330.530.560.360.480.010.420.840.620.21–0.40
199610.860.070.32110.610.3110.710.41–0.60
199710.540.950.450.2911110.530.7410.61–0.80
199810.2710.2610.30.30.370.990.290.540.610.81–1.00
19990.990.280.450.680.40.860.310.930.470.480.290.3
20000.260.280.010.950.320.380.260.320.010.370.660.43
Table 4. Ability of fitted full models for each individual year (rows) to predict bark beetle forest infestation in the years 1994 to 2000 (columns), except the same year for which the model was fitted, evaluated using AUC. Higher AUC values mean higher predictive power of the respective models; value 0.5 represents random prediction.
Table 4. Ability of fitted full models for each individual year (rows) to predict bark beetle forest infestation in the years 1994 to 2000 (columns), except the same year for which the model was fitted, evaluated using AUC. Higher AUC values mean higher predictive power of the respective models; value 0.5 represents random prediction.
Years of PredictionLegend
1994199519961997199819992000
Years of Model Fit1994-0.760.750.660.490.500.62≤0.5
19950.79-0.770.710.600.520.470.51–0.60
19960.810.83-0.730.540.510.660.61–0.70
19970.730.750.79-0.650.570.530.71–0.80
19980.650.680.630.69-0.570.760.81–0.90
19990.690.620.630.670.65-0.82
20000.700.660.630.670.650.65-
Table 5. Difference between AUC values corresponding to the full model set (with both Wetness and Brightness slopes) and the plain model set (with no remote sensing variable).
Table 5. Difference between AUC values corresponding to the full model set (with both Wetness and Brightness slopes) and the plain model set (with no remote sensing variable).
Year of Prediction
1994199519961997199819992000
Year of Model Fit1994-0.010.010.020.01−0.01−0.04
19950.03-0.020.020.03−0.01 0.01
19960.040.01-0.010.01 0.01−0.02
19970.060.040.03-0.02 0.03 0.03
19980.060.070.050.01- 0.12 0.02
19990.100.060.060.020.02-−0.02
20000.000.000.000.000.00 0.00-

Share and Cite

MDPI and ACS Style

Hais, M.; Wild, J.; Berec, L.; Brůna, J.; Kennedy, R.; Braaten, J.; Brož, Z. Landsat Imagery Spectral Trajectories—Important Variables for Spatially Predicting the Risks of Bark Beetle Disturbance. Remote Sens. 2016, 8, 687. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080687

AMA Style

Hais M, Wild J, Berec L, Brůna J, Kennedy R, Braaten J, Brož Z. Landsat Imagery Spectral Trajectories—Important Variables for Spatially Predicting the Risks of Bark Beetle Disturbance. Remote Sensing. 2016; 8(8):687. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080687

Chicago/Turabian Style

Hais, Martin, Jan Wild, Luděk Berec, Josef Brůna, Robert Kennedy, Justin Braaten, and Zdeněk Brož. 2016. "Landsat Imagery Spectral Trajectories—Important Variables for Spatially Predicting the Risks of Bark Beetle Disturbance" Remote Sensing 8, no. 8: 687. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080687

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop