Evaluation of the Influence of Disturbances on Forest Vegetation Using the Time Series of Landsat Data: A Comparison Study of the Low Tatras and Sumava National Parks

Abstract

This study focused on the evaluation of forest vegetation changes from 1992 to 2015 in the Low Tatras National Park (NAPANT) in Slovakia and the Sumava National Park in Czechia using a time series (TS) of Landsat images. The study area was damaged by wind and bark beetle calamities, which strongly influenced the health state of the forest vegetation at the end of the 20th and beginning of the 21st century. The analysis of the time series was based on the ten selected vegetation indices in different types of localities selected according to the type of forest disturbances. The Landsat data CDR (Climate Data Record/Level 2) was normalized using the PIF (Pseudo-Invariant Features) method and the results of the Time Series were validated by in-situ data. The results confirmed the high relevance of the vegetation indices based on the SWIR bands (e.g., NDMI) for the purpose of evaluating the individual stages of the disturbance (especially the bark beetle calamity). Usage of the normalized Landsat data Climate Data Record (CDR/Level 2) in the research of long-term forest vegetation changes has a high relevance and perspective due to the free availability of the corrected data.

1. Introduction

The issue of the time series has been a highly discussed subject recently, particularly in the field of observations of local and global change. The time series from Landsat gives us insights from the 1980s that can be used for a variety of analyses. This is the only continuous mission in the field of high-resolution data. This fact brings many problems that need to be addressed through the pre-processing of data: from the non-homogeneity of the sensed environment, through the different types of sensors with different radiometric and spectral resolutions, to the aging of the sensors. Several possible methods of processing the long time series are currently being used. Due to the low temporal resolution (for example, unlike Sentinel-2) of the Landsat images, it is necessary to pre-process the data for a possible quantitative comparison. We can use, for example, spectral unmixing methods [1], fractional methods [2], radiometric normalization and linear image regression [3,4,5,6,7,8,9,10,11,12,13]. Other methods such as cross-calibration methods that combine sets of multiple data types [14] are also used.

Thanks to the correct pre-processing of the raw data, it is possible to monitor several local and global phenomena. From the local phenomena, changes in forest ecosystems have been the most observed. In forest ecosystems, it is possible to observe the decrease of forests in tropical areas, the cultivation of palm trees or the range of biotic and abiotic disasters, among which bark beetle calamities are associated, which is also pertinent to this article. The massive expansion of bark beetle calamities was mainly due to the planting of non-indigenous species of coniferous forests in Central Europe. In the last few years, there have been large droughts that evoke dissemination of the pests to damage the weakened vegetation [15] and also allows the pests to dig out several times in a year. The second major factor is global change and global warming, which helps to exacerbate calamities.

Current studies [16] help us determine the parameters of the environment into which the pests are spreading. The quality of the trees, altitude, soil quality, groundwater content, or aspect belong among the factors that affect dissemination. Besides medium multispectral data, hyperspectral data [17] or very high-resolution multispectral data [18] are also used. These data help us determine the individual phases of damage more accurately, such as the Red Attack phase. High-resolution multispectral data help us locate the areas by various indicators and detect the particular phase or quantity of damage [19,20,21]. For medium resolution, only the wider changes [22,23,24,25] can be localized. On this basis, this study documented the process of bark calamities in various phases in the selected localities of two national parks using a time series of satellite imagery that was pre-processed with relative radiometric normalization for the purpose of comparison of values over time. As an indicator of change, specific types of vegetation indices were selected.

This work was focused on evaluating of changes in the forest areas in the selected localities of the Low Tatras NP (National Park) and Sumava NP using TS methods. The creation of the TS is based on the Landsat images. The changes in the forest vegetation induced by various disturbance processes were evaluated based on the selected vegetation indices, their tentative capabilities being tested and validated by comparing the in-situ data. An important part of this study was the evaluation of the data, methods and results achieved for the landscape and nature management and nature conservation institutions.

The main objectives of this work were:

• To evaluate and compare changes in the forest areas in the selected localities of the Low Tatras National Park and Sumava National Park using TS methods, which use normalized Landsat data
• To evaluate the suitability of individual vegetation indices for the detection of different types of biotic and abiotic disturbances
• To validate and interpret the results using in-situ data
• To discuss and recommend the suitability of the Earth Observation for nature conservation and management of the Low Tatras National Park and Sumava National Park.

Observed Area

The Low Tatras National Park (in Slovakia) is 73 km². The dominant land cover is forest; thus, we can find extensive forest ecosystems. In many parts of the Low Tatras, spruce forests dominate. Their development takes place under the influence of natural and anthropogenic factors. The integrity of the forest ecosystems has been disturbed by several abiotic events (wind, snow, frost, and avalanches), biotic agents (sub-insects) or anthropogenic influences (e.g., air pollution).

On 19 November 2004, from 15:00 to 24:00, a windstorm passed through the area of the Low Tatras with a maximum speed of about 175 km/h. The territory of the Low Tatras National Park and its protection zone was damaged and extensive forest stands in several areas and localities were destroyed. The wind storm in November 2004 hit the territory of NAPANT very strongly and, as a result, led to the reproduction of the European bark beetle (Ips typographus). During the summer of 2007, extremely favorable conditions for grafting of subcortical insects were created in almost all non-native and natural spruces in Slovakia (Figure 1). The observed sites of this study are located in the Dumbier part of the Low Tatras (Western Carpathians), in the Mlynna valley. It is in an about 7 km long valley. After the wind calamity in 2004, the valley area was affected by a biotic bark beetle calamity, mainly culminating in 2009 (mass pine devastation). The sites of interest were categorized according to the type of disturbance. The first site represents the territory with the wind disturbance in 2004, the second represents the prevailing beetle disturbance between 2006 and 2009, the third site has minimal disturbance, the fourth site represents a place with the beetle disturbance (2005) immediately after wind calamity in 2004 and the last site represents a place with minimal disturbance (Figure 2 and

Table 1. Overview of the selected localities in the Low Tatras NP with the GPS position in WGS 84.

ID	Type of Disturbance	Year	Northing	Easting
1	Wind calamity	2004	48.9047556	19.6364333
2	Bark beetle calamity	2006–2009	48.9155019	19.6520547
3	Minimal disturbance	–	48.9010178	19.6532992
4	Bark beetle with wind calamity	2004–2005	48.903718	19.716530
5	Minimal disturbance	–	49.0019478	19.6275692

).

The Sumava is an extensive mountain range on the border of Czechia, Austria and Germany (Bavaria). The Sumava (Bohemian) Forest together with the neighboring Bavarian Forest creates the most extensive forest landscape in Central Europe, called the “Green Roof of Europe”. The area of the park is more than 900 km². The Sumava forests are dominated by spruce.

During the 1990s in the last century, a bark beetle calamity occurred in the Czech part of Sumava NP with the culmination during 1995–2001. The Kyrill gales had a strong impact on the forest vegetation in 2007. In 2008, there was a sharp rise in the bark beetle affected forest areas and it continued to grow until 2010. In the period from 2011 to 2014, there was a slight decline in the intensity of the disturbance with a stabilization period after 2014 (Figure 3). The first site of interest in this study represents the territory with the wind and bark beetle disturbance, the second and fourth represent the bark beetle disturbance, and the third and fifth site are with minor disturbances (Figure 4 and

Table 2. Overview of the selected localities in the Sumava NP with the GPS position in WGS 84.

ID	Type of Disturbance	Year	Northing	Easting
6	Bark beetle calamity with wind and non-natural recovery	since 2007	48.973628	13.561746
7	Bark beetle calamity and natural recovery	since 2008	48.983655	13.561720
8	Minimal disturbance	–	48.989884	13.558225
9	Bark Beetle with natural recovery	since 2009	48.9847661	13.5227772
10	Minimal disturbance	–	49.0451933	13.4761361

).

2. Materials and Methods

2.1. Data

Satellite and in-situ data were used to process the study. Landsat Climate Data Record (CDR) satellite imagery was used as the remote sensing data. The distributed Landsat CDR data were pre-processed with geometric and atmospheric corrections using the CDR database. For more information on the methodologies, see, e.g., [26,27].

A sufficient number of usable Landsat satellite images (data from the Landsat 4–8 missions) were found for our areas of interest. Eleven images were finally selected for each of the Low Tatras NP (

Table 3. Overview of the used Landsat images in the Low Tatras NP.

Name	Date	Sensor
LC81880262013219LGN00	7 August 2013	Landsat 8
LC81880262015193LGN00	12 July 2015	Landsat 8
LE71880261999221SGS01	9 August 1999	Landsat 7
LE71880262001242SGS00	30 August 2001	Landsat 7
LT51880261994183XXX02	2 July 1994	Landsat 5
LT51880262005245KIS00	2 September 2005	Landsat 5
LT51880262006200KIS01	19 July 2006	Landsat 5
LT51880262007203MOR00	22 July 2007	Landsat 5
LT51880262009240KIS00	28 August 2009	Landsat 5
LT51880262011198MOR00	17 July 2011	Landsat 5
LT41880261992202XXX02	20 July 1992	Landsat 4

) and the Sumava NP (

Table 4. Overview of the used Landsat images in the Sumava NP.

Name	Date	Sensor
LC81920262013215-SC20160702161259	3 August 2013	Landsat 8
LC81920262015221-SC20160702161912	9 August 2015	Landsat 8
LE71920262002209-SC20160702153836	28 July 2002	Landsat 7
LT51920261994211-SC20160702152112	30 July 1994	Landsat 5
LT51920261998222-SC20160702151843	10 August 1998	Landsat 5
LT51920262004223-SC20160702151849	10 August 2004	Landsat 5
LT51920262005241-SC20160702151622	29 August 2005	Landsat 5
LT51920262006196-SC20160702151531	15 July 2006	Landsat 5
LT51920262007231-SC20160702151832	19 August 2007	Landsat 5
LT51920262009236-SC20160702151537	24 August 2009	Landsat 5
LT51920262010191-SC20160702151637	10 July 2010	Landsat 5

). A crucial condition for the selection of the data was their minimum cloud cover and no significant haze infection. Another important aspect of the selection was the date of the acquisition. According to many authors (e.g., Vogelmann et al. [28] and Griffiths et al. [29,30]), the summer to the first autumn phase is an appropriate time for the assessment of the forest vegetation. Based on these works, it was decided to select data in the months of July to September.

The in-situ data, used for the validation of the results from the Landsat data analysis, were obtained from the administration of the national parks and from their own field survey. The information obtained from the national parks included the field records of foresters, information from forest management plans and archive records of nature conservation documentation. The forest management plan is legislatively enshrined in the Forest Law on Forest Economic Planning and it provided data on the status of the forests and their past management: the area and category of the forest, the age of the stock and the age, stocking and representation of the individual trees, the average height, stockpile and method of management and proposal of the economic measures.

2.2. Data Processing

To determine the state and changes in the vegetation from the satellite images, the key information was placed on the spectral characteristics of the vegetation species studied. Each species of vegetation has specific characteristics, and we can detect them and determine their properties (e.g., health status) based on this knowledge. Spectral reflection can be represented by a curve that can be divided into three basic parts according to the main structural properties: the area of pigmentation absorption, the cellular structures and the area of water absorption [31]. The time series analyses monitor the changes in the vegetation over long periods of time. The spectral characteristics of studied objects are often used for estimating the health phase of vegetation or it could be based on the spectral characteristics of the derived data, such as the vegetation indices.

To ensure compatibility, the Landsat data source was used. As mentioned above, the atmospheric and radiometric corrections did not need to be performed on the CDR data because the images are distributed with the correction and conversion to surface reflectance (i.e., surface reflectance). The Fmask [32] feature was used to create cloud-free images and to mask the clouds, shadows, snow and other disturbing elements.

It is clear from many studies that a significant problem in creating the time series is to ensure compatibility between the types of data that are often taken by multiple types of sensors [33]. A different sensor type and a different acquisition time may have a large influence on the result of the time series analyses. For this reason, relative radiometric normalization was used to eliminate the influence of the different acquisition times, the different phenological phases of the vegetation and the influence of the different spectral and radiometric characteristics of the different sensors. Thanks to normalization, we can reduce the effects of the different time and place of acquisition, the different positions of the sun, and the different radiometric and spectral differences.

Based on previous empirical testing, the PIF (Pseudo-Invariant Features) Linear Based method was selected for the normalization purposes. This normalization was selected based on the recommendations of many studies [33,34] and the testing of relative radiometric normalizations in the article by Lastovicka et al. (2017) [35]. It is a relatively standard proven method, which already has its origins in the 1990s, and the PIF Linear Based method is probably the most commonly used [34]. There are many variations of PIF and it is used in many programs for time series (IR-MAD/MAD CAL algorithm). The linear-based method is based on linear alignment using linear regression. For the PIF linear-based method, a classical linear regression (y = b₁ * x + b₂ + ε) was used, which represents the approximation of the known values using the smallest square method.

Thus, by linear regression, we looked for the most suitable parameters, b₁ and b₂, to approximate the two observed frames in the correlation diagram. The values b₁ and b₂ represent the slope and intercept, respectively [34].

The method consists of finding the pseudo-invariant objects and the invariant objects in two observed images that occur as clusters [33,36,37,38]. Using a scatter plot, the values of both images and the search for the unmodified objects were displayed manually or automatically using the algorithm [39]. The invariant areas required for this method were selected based on the NDVI index, where the minimum area of the index was the appropriate area. Roads, parking areas, roofs and other long-standing elements were among the selected invariant objects. These selected elements are suitable for normalization calculations. Then, the regression parameters were computed [33]. The equalization was performed for each spectrum band separately—the band-by-band method [39].

For the data normalization of the CDR Landsat data, custom applications were developed in the MATLAB environment. The developed application was inspired by the TimeSync web application (made by Oregon State University—http://timesync.forestry.oregonstate.edu [40]), which has limitations in accessing custom data as well as in bulk data over the Internet. Our application allows one to work with any type of satellite data in an off-line mode. Thanks to the relatively standardized MATLAB encryption, this application is transferable to a wide range of software and operating systems. The application itself consists of several algorithms for which two user interfaces were created. The first user interface allows seeing the differences between the two spectral bands that are needed for normalization purposes. The scatter graph and interdependence were displayed using a linear regression curve.

After the data standardization process, the selected vegetation indices were calculated, and the time series charts were created. To create the curves, the developed applications were used. For the specified location, the application created time series charts. The calculations take the pixel values of the site or its surroundings (3 × 3 pixels) into account. In this study’s case, the closest neighborhood of the site was used (3 × 3 pixels). Nine pixels were selected for one pixel in the neighborhood of the central pixel (the localities were chosen from a larger area of disturbed or non-disturbed places, thus the values were averaged for the approximate value errors or outlying values). For a more detailed description of this application, see the work of Lastovicka et al. [35].

Based on several elaborated studies [34,41,42,43,44], it was decided to use the list of vegetation indices shown in

Table 5. List of the vegetation indices and their expressions.

Vegetation Index	Shortcut	Equation
Foliar Moisture Index	FMI (value * 10)	(NIR)/(RED * SWIR)
Normalized Difference Moisture Index	NDMI	(NIR − SWIR)/(NIR + SWIR)
Normalized Difference Vegetation Index	NDVI	(NIR − RED)/(NIR + RED)
Simple Ratio Index	SR (value/100)	(NIR)/(RED)
Transformed Vegetation Index	TVI	sqrt ((NIR-RED)/(NIR + RED) + 0.5)
Wide-band Normalized Difference Infrared Index	wNDII	(2 * NIR − SWIR)/(2 * NIR + SWIR)

Note: The FMI and SR values were optimized for the value level to be comparable with the other chosen indices in the graphs.

.

The masking and calculation of the vegetation indices was performed in ENVI 5.3. The final database created for the time series analysis included calculations of six vegetation indices from the Landsat CDR normalized data by the PIF Linear Based method.

For the statistical approach, the minimum (MIN), the maximum (MAX), the difference between maximum and minimum (MAX-MIN) and the standard deviation (St. Deviation) were calculated.

For both areas of interest, the aggregate mathematical statistical characteristics were calculated. From the locations where calamities occurred, the average values for the given period were calculated and the graphs for all the interest areas were plotted.

3. Results

This section focuses on the description and interpretation of the value of the vegetation indices in the 10 different localities in the Low Tatras and Sumava National Parks that have been under different development influencing the forest areas (damage from wind calamity, bark beetle and forest vegetation without any significant influence). The results describe and interpret the values of the studied vegetation indices during the observed period (1992–2015) and evaluate their suitability for the assessment of the forest changes. The results should indicate the applicability of the different vegetation indices or their combinations for the detection of the different types of disturbances.

3.1. Locality 1

A significant part of the Low Tatras National Park was damaged by the wind calamity of Elizabeth, which took place on 19 November 2004. The first site represents the area heavily damaged by this event, followed by the natural regeneration of the vegetation in the following years after the calamity.

The results of the development of the individual monitored indices are documented in Figure 5. The evolution of the values of the monitored indices shows that the wind calamity that occurred in 2004 significantly affected the state of the forest vegetation. The examined indices strongly reflected this wind calamity. It can be seen in Figure 5 that the most significant change (decrease) is in the values of the following vegetation indices NDVI, NDMI and wNDII; the decline in the curve is less significant for the FMI, SR and TVI indices. These results can be also seen from the statistical information in the Figure 5b, similar figures were calculated for all sites. For areas with disturbance, there was detected a significant difference between minimum and maximum. In comparison, this statistical indicator was lower in the areas with minimal disturbance. These results were also achieved with standard deviation. Places with minimal disturbance had lower standard deviation. The area gradually began to be covered by the ascending deciduous woods (after harvesting the dead trees), especially Sorbus aucuparia. These young deciduous trees helped with the regeneration of the vegetation cover in the damaged area.

In 2005, there is missing value in Figure 5a. This is in case of clouds in the satellite image (Figure A1 and Figure A2 in Appendix A). In this situation, linear interpolation was used, and the curve was connected to the nearest value point. The same method was used in other localities, if there was a problem with missing values due to cloud cover.

3.2. Locality 2

The second observed site is a site affected by the strongly expanding, devastating bark beetle after 2007. The forest vegetation after the attack was dying and left to spontaneously develop. The TS graph of the vegetation indices in this area can be seen in Figure 6.

In Figure 6, we can see some obvious trends. At the initiation phase of the bark beetle calamity, the vegetation indices wNDII, NDMI and FMI declined, while the values of the NDVI, SR and TVI indices slightly increased in the period 2005–2007. During 2007–2009, the location was characterized by the significant degradation of the spruce forest. All the monitored indices responded with a decrease in the values continuing until 2012. Then, the natural regeneration of the territory with the dominance of the deciduous trees began, which successively occupied the deforested areas in the territory of the calamity. This was mainly reflected by the NDMI, wNDII and NDVI indices.

3.3. Locality 3

The next observed locality in the Low Tatras represents a place where, unlike the two previous localities, there have been no significant changes in the vegetation during the observed period. The time series of the vegetation indices can be seen in Figure 7.

It can be seen in the figure that relatively small value changes occurred in the observed area. Insignificant oscillations are mainly evident in the wNDII and NDMI indices, as well as in the FMI. A slight fluctuation was recorded in the TVI index, less so in the NDVI and SR indices. In terms of interpretation, it is necessary to realize that changes in the index values may reflect specific, local conditions, e.g., weather conditions: precipitation vs. drought [39].

3.4. Locality 4

The fourth observed site is a site affected by the strongly expanding, devastating bark beetle after 2001. The forest vegetation after the attack was dying and left to spontaneously develop. The TS graph of the vegetation indices in this area can be seen in Figure 8. All monitored indices responded with a decrease in the values with minimal values in 2005. Then, the natural regeneration of the territory was reflected by an increase of the values, mainly in NDMI, wNDII and NDVI indices.

3.5. Locality 5

Locality 5 (Figure 9) in the Low Tatras represents a place where, unlike Localities 1, 2 and 4, there have been no significant changes in vegetation during the observed period. As can be seen in the graphs, values have a fixed level. Only in 2007 we can see an insignificant decline of values of the indices.

3.6. Locality 6

The first site from the Sumava NP is a site that has been hit by a wind calamity (2007). After disturbing the tree structures, the forest was much more affected by bark beetle attack. The result can be seen in Figure 10. Here, the steep decline in the curves after 2008 is especially visible, caused by the disturbance of the trees by the wind and the subsequent bark beetle. The destruction process is so much faster than, for example, in Locality 7.

3.7. Locality 7

The next location from the Sumava NP is similar to Locality 2 from the Low Tatras region. This is a place where a biotic disturbance has occurred, i.e., a bark beetle calamity. Figure 11 shows the results of the time series analyses. The NDMI and wNDII indices reflected the whole disturbance event with a decreasing value. The NDVI, SR, and TVI indices had similar trends in the first years of the calamity (decrease). The dead trees were not harvested.

3.8. Locality 8

The next site of interest is similar to Locality 3 in the Low Tatras NP: without any significant disturbance. This location is covered by a spruce forest with an approximate age of 50 years according to the in-situ data. In Figure 12, we can see the development with the oscillation in the values without any significant breaks.

3.9. Locality 9

The next location from the Sumava NP is similar to Locality 2 from the Low Tatras region and Locality 7 from the Sumava National Park. This is the place where a biotic disturbance has occurred, i.e., bark beetle calamity. Figure 13 shows the results of the time series analyses. The NDVI, NDMI and wNDII indices reflected the whole disturbance event with a decreasing value. The dead trees were not harvested after the disturbance.

3.10. Locality 10

The last site of interest is similar to Localities 3 and 5 in the Low Tatras NP (without any significant disturbance) and Locality 8 in the Sumava National Park (with minimal disturbance). This location is covered by a spruce forest with an approximate age of 60 years according to the in-situ data. In Figure 14, we can see the development with the oscillation in the values with only small breaks.

3.11. Statistical Results and Comparison of Both Specific Areas

Figure 15, Figure 16, Figure 17 and Figure 18 show the average values from the Low Tatras National Park and the Sumava National Park. The graphs 16 and 18 show the minimal variance caused by the different site conditions. However, in Figure 15 and Figure 17, we can see the graphs of both sites affected by the calamities with average values. For Figure 17 (the Sumava National Park), it can be stated that, due to the gradual spread of the disturbance in the individual years, the longer span of the time of the calamity breakdown shows the calamity break is slower. In addition to the summary graphs, the individual aggregate values of the individual vegetation indices are described in

Table 6. Statistical information about the average values in the disturbed areas in the Low Tatras NP.

Vegetation Index	MIN	MAX	MAX-MIN	Standard Deviation
FMI (value * 10)	0.0412	0.1676	0.1264	0.0421
NDMI	0.1612	0.4191	0.2579	0.0883
NDVI	0.6712	0.8249	0.1537	0.0551
SR (value/100)	0.0560	0.1078	0.0518	0.0177
TVI	1.0805	1.1510	0.0705	0.0252
wNDII	0.4683	0.6555	0.1872	0.0710

and

Table 7. Statistical information about the average values in the areas with minimal disturbance in the Low Tatras NP.

Vegetation Index	MIN	MAX	MAX-MIN	Standard Deviation
FMI (value * 10)	0.0932	0.1980	0.1048	0.0273
NDMI	0.3283	0.4213	0.0930	0.0262
NDVI	0.8108	0.8812	0.0705	0.0227
SR (value/100)	0.0984	0.1678	0.0694	0.0212
TVI	1.1448	1.1752	0.0304	0.0098
wNDII	0.5948	0.6614	0.0666	0.0200

for the Low Tatras, and

Table 8. Statistical information about the average values in the disturbed areas in the Sumava NP.

Vegetation Index	MIN	MAX	MAX-MIN	Standard Deviation
FMI (value * 10)	0.0334	0.1411	0.1077	0.0374
NDMI	0.1341	0.4411	0.3070	0.1065
NDVI	0.6812	0.7658	0.0846	0.0309
SR (value/100)	0.0532	0.0762	0.0231	0.0079
TVI	1.0868	1.1250	0.0383	0.0141
wNDII	0.4454	0.6745	0.2291	0.0797

and

Table 9. Statistical information about the average values in the areas with minimal disturbance in the Sumava NP.

Vegetation Index	MIN	MAX	MAX-MIN	Standard Deviation
FMI (value * 10)	0.0716	0.1694	0.0978	0.0300
NDMI	0.3433	0.4737	0.1304	0.0395
NDVI	0.7051	0.7922	0.0871	0.0290
SR (value/100)	0.0593	0.0866	0.0273	0.0094
TVI	1.0976	1.1367	0.0391	0.0130
wNDII	0.6061	0.6952	0.0891	0.0271

for the Sumava National Park.

4. Discussion and Conclusions

The forests in the Low Tatras and Sumava National Parks have been severely affected by wind calamities and, subsequently, by bark beetle insects. For the selected locations, which differed in the type of disturbance, the calculation of the vegetation indices and the comparison/evaluation of their trajectories were performed. The main aim was to determine the suitability of the individual vegetation indices for the detection of the different types of biotic and abiotic disturbances using a TS analysis. Based on the results of this study, it is obvious that each type of disturbance, whether biotic or abiotic, has a specific development and different consequences on the state (health) of the vegetation/spectral expression. In the case of a wind disaster, a large part of the forest vegetation is devastated in a quick period. This corresponds to a large impact decrease (more than with the bark beetle calamity, because the biotic calamity is decreasing more slowly) in the indices using RED, NIR and SWIR bands [45]. A forestry intervention or succession of new vegetation causes an increase in the reflection in the green and especially in the near infrared part of the spectrum in a relatively short time after the disturbance. A wind calamity has a strong effect on the vegetation, however the recovery in the following years is quick, both in a spontaneous way and a forest-controlled way. When the three sites affected by wind calamities in Low Tatras and Sumava NPs were compared, similar result could be seen. The vegetation indices based on the use of the SWIR band (e.g., the NDMI and wNDII index) show similar trends as well [46]. A bark beetle calamity has a much more complicated progress with a specific effect on the vegetation. Several generations of beetles gradually occupy the forest with a strong culmination phase. There can be three generations of beetles in one year, which greatly disturb the resistance of the forest stand. Our results confirmed the high sensitivity of the wNDII and NDMI vegetation indices based on the SWIR band in the case of the bark beetle disturbance (see Figure 17). These vegetation indices were able to reflect all the phases of the disturbance using the time series analysis. The reflectivity of the healthy forest vegetation in the SWIR was lower than in the NIR. However, during the disturbance/damage, the SWIR reflectivity tended to increase. It is clear from the results that this aspect played an important role in detecting the disturbed period (compare the NDMI or wNDII vegetation indices that use the SWIR band with the vegetation indices based on NIR such as NDVI). An important result of this study is that the NDMI and wNDII vegetation indices can detect all the stages of the bark calamity (before, during and after the culmination) and are thus suitable indicators in monitoring the health status of the forest vegetation in such specific/environmentally complex processes as a bark beetle calamity with its causes and consequences. On the other hand, it can be stated that indices using SWIR are more sensitive to more local factors [45]. The values of the vegetation indices can be affected by many factors out of our current investigation, e.g., meteorological factors such as precipitation or soil moisture. For a more detailed interpretation of these changes, it would be necessary to include several factors, such as meteorological, which was beyond the scope of this study.

In terms of comparison of the two observed areas (Sumava NP and Low Tatras NP), statistical summaries were made. It can be said that the calamity conditions of the observed sites took place over a longer period in Sumava areas. This can also be seen from the individual figures. However, if we look at aggregate indicators, both locations have a very different evolution of recovery phase. In particular, in the Low Tatras, trees were harvested after the calamities. This results in a different development from the Sumava areas, where the trees were left. The spread of new trees and bushes in Sumava is slower and the values of the selected indices after overcoming the recovery phase do not reach their initial values from before the calamity, but in the Low Tatras in several cases indices approximately reach their initial values from before the calamity. However, this does not indicate anything about an unnatural interference with nature. In Figure 15, Figure 16, Figure 17 and Figure 18 and Table 6, Table 7, Table 8 and Table 9, the NDMI index can be picked up again as the best index to detect the disturbance of the forest ecosystems. Secondly, the wNDII index also had comparable results. The NDVI index also corresponded to the disturbances, but the problem of the NDVI index is that it cannot differentiate the recovery phase in comparison with the wNDII and NDMI indices.

A long-term TS can be used to track the development and current state of diverse forest ecosystems. The Landsat data are a very useful data source in this respect [47]. On the other hand, it is necessary to pay attention to the compatibility of the data coming from various types of Landsat sensors [40]. Although Landsat CDR data contain atmospheric and radiometric corrections, it is still useful to perform radiometric normalization. Thanks to the normalization, we can reduce the effects of the different times (phenology period) and places of acquisition, the different positions of the sun, and radiometric and spectral differences. The PIF Linear Based method appears to be a suitable standardization method. For data normalization of the CDR Landsat data, custom applications were developed in the MATLAB environment. The application allows one to work with any type of satellite data and it is transferable to a wide range of software and operating systems.

This study should serve as an example when the development of the impacts of the disturbances on the state of the forest vegetation was analyzed by traditional methodological procedures using vegetation indices. These results and methods should be useful, inspirational to the institutions for the management of the protected areas and to determine the appropriate forest management practices implemented in them. The time series results obtained can be applied to generate predictions of future states and to observe the spontaneous and economically managed forest regeneration. The use of satellite data provides progressive opportunities to monitor and evaluate the state and changes in the forest vegetation for the purpose of the protection and management of national parks [46,47].