Sentinel-2 data in an evaluation of the impact of the disturbances on forest vegetation

remote sensing Article Sentinel-2 Data in an Evaluation of the Impact of the Disturbances on Forest Vegetation Josef Lastovicka 1, Pavel Svec 2, Daniel Paluba 1, Natalia Kobliuk 1, Jan Svoboda 1, Radovan Hladky 1,3 and Premysl Stych 1,* 1 Department of Applied Geoinformatics and Cartography, EO4Landscape Research Team, Faculty of Science, Charles University, 12843 Prague, Czech Republic; [email protected] (J.L.); [email protected] (D.P.); [email protected] (N.K.); [email protected] (J.S.); [email protected] (R.H.) 2Department of Geoinformatics, Faculty of Mining and Geology, VSB–Technical University of Ostrava, 70800 Ostrava-Poruba, Czech Republic; [email protected] 3Administration of the Low Tatras National Park, State Nature Conservation of the Slovak Republic (NAPANT), 97401 Banska Bystrica, Slovakia *Correspondence: [email protected] Received: 8 April 2020; Accepted: 23 May 2020; Published: 13 June 2020   Abstract: In this article, we investigated the detection of forest vegetation changes during the period of 2017 to 2019 in the Low Tatras National Park (Slovakia) and the Sumava National Park (Czechia) using Sentinel-2 data. The evaluation was based on a time-series analysis using selected vegetation indices. The case studies represented five different areas according to the type of the forest vegetation degradation (one with bark beetle calamity, two areas with forest recovery mode after a bark beetle calamity, and two areas without significant disturbances). The values of the trajectories of the vegetation indices (normalized difference vegetation index (NDVI) and normalized difference moisture index (NDMI)) and the orthogonal indices (tasseled cap greenness (TCG) and tasseled cap wetness (TCW)) were analyzed and validated by in situ data and aerial photographs. The results confirm the abilities of the NDVI, the NDMI and the TCW to distinguish disturbed and undisturbed areas. The NDMI vegetation index was particularly useful for the detection of the disturbed forest and forest recovery after bark beetle outbreaks and provided relevant information regarding the health of the forest (the individual stages of the disturbances and recovery mode). On the contrary, the TCG index demonstrated only limited abilities. The TCG could distinguish healthy forest and the gray-attack disturbance phase; however, it was difficult to use this index for detecting different recovery phases and to distinguish recovery phases from healthy forest. The areas affected by the disturbances had lower values of NDVI and NDMI indices (NDVI quartile range Q 2 –Q 3 : 0.63–0.71; NDMI Q 2 –Q 3 : 0.10–0.19) and the TCW index had negative values (Q 2 –Q 3 : − 0.06– − 0.05)). The analysis was performed with a cloud-based tool—Sentinel Hub. Cloud-based technologies have brought a new dimension in the processing and analysis of satellite data and allowed satellite data to be brought to end-users in the forestry sector. The Copernicus program and its data from Sentinel missions have evoked new opportunities in the application of satellite data. The usage of Sentinel-2 data in the research of long-term forest vegetation changes has a high relevance and perspective due to the free availability, distribution, and well-designed spectral, temporal, and spatial resolution of the Sentinel-2 data for monitoring forest ecosystems. Keywords: time series; Sentinel-2; vegetation index; bark beetle disturbance; Czechia; Slovakia; cloud computing; big data; temporal resolution Remote Sens. 2020,12, 1914; doi:10.3390/rs12121914 www.mdpi.com/journal/remotesensing Remote Sens. 2020,12, 1914 2 of 26 1. Introduction One possibility to monitor and evaluate forest vegetation is to use Earth Observation (EO) methods. EO provides unique information for the purpose of observing dynamic phenomena on the Earth’s surface [ 1 ]. Satellites are equipped with various sensors that provide images with information that the human eye cannot see. This is especially useful for detecting the early stages of forest disturbance when there are no visual signs of damage [ 2 ]. One of the most beneficial outputs based on EO is the worldwide databases of status and changes of forest areas [ 3 ], the databases of forest area changes for Eastern Europe [ 4 ], or the databases of systematically addressed changes in forest areas, both in North and South America (e.g., in Brazil and Colombia) [ 1 , 5 – 9 ]. Other authors evaluated a forest affected by bark beetles with vegetation indices from Landsat images, and the subsequent pixel classification achieved a total accuracy of 80–82% [ 10 ] or focused on the classification of three categories of the forest using the maximum likelihood classifier: healthy trees, damaged trees, and grasslands [11]. The forest bark beetle calamity itself can be divided into the green-attack (without visible damage), followed by a red-attack (the first visible phase of the death of trees), and finally a gray-attack (dead trees). This method successfully separated the red-attack trees from the healthy trees and the grasslands using the RGI (red green index) and reflectivity in the green band with an overall classification accuracy of 86%. The SWIR (short-wavelength infrared) band effect was demonstrated in detecting forest disturbances using time series (TS) methods [ 12 ]. In Central Europe, evaluations of the state and changes in forest vegetation using Landsat have been presented by many studies [ 13 – 16 ]. EO data have often been used in the evaluation of forests in Czechia, which have been struggling with long-term disturbances (wind calamities or the subsequent bark beetle invasions) [17–21]. To determine the state and the changes of the forest vegetation from the satellite images, the key information was found in the spectral properties of the vegetation species studied [ 22 ]. The suitability of the multispectral data and the vegetation indices calculated from these could differ in the monitoring of the forest vegetation and how it was proven in various studies [ 20 , 23 – 28 ]. The unanswered question is of which real spectral differences measured during the observed years are relevant to determine the characteristics and individual phases of the forest disturbance. TS methods are often used for the evaluation of the forest changes. For TS purposes, a whole range of satellite data can be used. For choosing a data type, the availability and suitability of the images, especially in terms of time and the radiometric and spatial resolution, need to be considered [29]. Due to the availability of the free access archive, Landsat’s satellite imagery archive is one of the most commonly used archives for a wide range of disciplines [ 30 ]. Landsat data provides visible, NIR (near-infrared) and SWIR bands that have often been used in the evaluation of forest degradation affected by a disturbance, e.g., the bark beetle calamity. Landsat TM data have proved to be effective for detecting the red-attack phase of the bark beetle disturbance [ 31 ]. From a temporal point of view, it is appropriate to use the Landsat mission data for a longer TS, especially as their time coverage is longer than 40 years. However, looking at the number of acquisitions per year, the 16-day temporal resolution is a limiting factor that is even lower due to cloud coverage, especially in mountainous areas [ 1 , 30 ]. In the case when there are only a few images during the year, then only one chosen reference image per year is often used [ 19 ]. Some studies tried to combine a range of different types of satellite data, e.g., Landsat with MODIS data [ 32 ]. In the case of a long TS period (more than 20 years), various radiometric and spectral resolutions of the acquired data can influence the results of the TS due to different types of sensors [33]. The Copernicus program has brought a new revolution in EO monitoring. ESA is developing new missions of satellites called Sentinels specifically for the operational needs of the Copernicus program. The images are received via two parallel missions 2A and 2B and in the case of the overlapping scenes, the temporal resolution is less than five days [ 34 ]. The Sentinel-2 multispectral optical dataset is now available with the ambition to provide data with better resolutions (spatial, temporal, and spectral) than the traditional data like Landsat. The Sentinel-2 data have been available since 2015; therefore, the archive is ready to be tested with the TS analyses. The spectral vegetation indices calculated Remote Sens. 2020,12, 1914 3 of 26 from Sentinel-2 have a potential for mapping and detecting changes induced by bark beetle attack, particularly based on the red-edge bands or water-related indices. These changes are limitedly detectable by Landsat-8 due to the lower spectral and spatial resolution of the OLI sensor; see the comparable study [2] with 67% accuracy for Sentinel-2 and 36% accuracy for Landsat 8. Several articles have been written for the TS of the Sentinel-2 data, e.g., focused on mapping floodplain grassland plant communities, where the authors dealt with the variability in the water content using TS, especially with TS using the support vector machine and random forest classifiers [ 35 ]. Other authors [ 36 ] focused on observing the Cotton Belt using the TS of the Sentinel-2 data. The authors used random forest classification and several kinds of vegetation indices for their time survey. The temporal and spatial resolution played an important role in the observation of the TS [ 37 – 39 ]. For this reason, a fusion of the Landsat 8 and Sentinel-2, called harmonized datasets, has been designed [40–42]. Few articles have been written for the TS of the Sentinel-2 focused on forest vegetation, e.g., that dealt with the recognition of unhealthy cork oaks (Quercus suber) [ 43 ]. Based on a multitemporal comparison (using vegetation indices), researchers found that while in the wetter part of the year, the differences were small, and, during the drier part of the year (September and October), the differences between the healthy and unhealthy cork oak stands were more significant. Another study focused on forest growth using the Sentinel-2 TS [ 44 ] and attempted to determine how to distinguish bamboo forests (Phyllostachys pubescens) from coniferous or deciduous forests using the annual course of the vegetation indices (NDVI, NDMI, etc.). The results prove that May is the best month to distinguish bamboo stands from other forest types. Another study [ 45 ] showed that the classification of forest types using the random forest method based on the Sentinel-2 TS proved a high relevancy of the red-edged bands. A perspective method in the processing of RS data is cloud computing. The Sentinel Hub is one of the most used cloud-based applications, and this application allows effective and end-user friendly processing and analysis of the EO data (https://www.sentinel-hub.com). It is a user interface app for the semi-automatic satellite data analysis created by Sinergise [ 46 ]. This platform provides satellite images from various missions, e.g., Sentinel, Landsat, or MODIS. The Sentinel Hub allows the downloading, the visualization of the satellite imageries, scripting to calculate the vegetation indices, classifications, and other analyses, as well as using WMS (Web Map Service) services to extract the selected values for the given site and period [46,47]. One of the current problems in the forest ecosystems of central Europe is bark beetle calamities [ 18 ]. Their overgrowth is closely related to the climate change and the spread of non-indigenous coniferous trees in central Europe, e.g., spruce (Picea abies) [ 48 ]. This study is focused on an evaluation of the changes in the forest vegetation in selected areas of the Low Tatras National Park (Slovakia) and the Sumava National Park (Czechia) using the TS methods. The analysis of the TS is based on Sentinel-2 images (description in Table A1; with a comparison to the Landsat data). Due to the dynamic changes in the forests in Czechia and Slovakia that occurred over the last several years, it is a suitable opportunity to apply the Sentinel-2 TS for the detection of disturbing events (bark beetle) and to evaluate the forest health status during the individual stages of the disturbance (forest conditions before, during, and after the disturbance events). From this point of view, the main aim and novelty of this study is to use and test the Sentinel-2 data for the evaluation of the dynamic forest changes. The selected vegetation indices and their trajectories of the TS were interpreted and validated in relation to the in situ data investigated during field research or provided by the Administration of National Parks. In the process of collecting and evaluating the data, we cooperated with the Low Tatras National Park and the Sumava National Park. An important part of this study is the evaluation of the applicability of the EO/Sentinel-2 data and the methods in the forest management and the nature protection of national parks. For this reason, big data/cloud data methods of the Sentinel Hub were used. Remote Sens. 2020,12, 1914 4 of 26 The main objectives of this study were: • To test the Sentinel-2 data in the TS analyses for the selected case studies in Czechia and Slovakia within a three-year period of 2017–2019 and to evaluate the benefits of the Sentinel-2 data for monitoring the forest changes. • To compare the Sentinel-2 and Landsat temporal and spatial resolution for the TS analyses of the forest vegetation in mountainous areas. • To evaluate the relevancy of the vegetation indices in the study of the forest changes and the health of the forest vegetation using the Sentinel-2 data. • To perform TS analyses and compare the results in the different types of areas affected by the bark beetle invasion (the disturbed and renewing forest ecosystems). • To discuss the positives and perspectives of the Sentinel-2 data in the TS of forest changes in comparison to the traditional data used, e.g., Landsat. • To process the satellite data and perform analyses in the cloud-based tool (Sentinel Hub) and discuss the positives and negatives of cloud-based systems for end-users in forestry research and management. We would like to give responses for these research questions: • What are the main positives of the Sentinel-2 data in the evaluation of the forest vegetation affected by the disturbances? What type of change in the forest is possible to detect by Sentinel-2? What are the positives using Sentinel-2 data in the TS in comparison with the traditional satellite data, e.g., Landsat? • What are the positives of processing and analyzing the data in a cloud-based tool (Sentinel Hub)? • What vegetation indices derived from the bands of the Sentinel-2 data are useful for the detection of the forest affected by a bark beetle invasion? Which vegetation indices based on these data could detect the disturbances that occurred and individual recovery phases in the forest of mountainous areas? Are the vegetation indices, traditionally used for Landsat data, usable for the Sentinel-2 data? • How many cloud-free images of Sentinel-2 are available for the TS analysis in our mountainous areas case studies annually? What is the progress of the data availability per year in comparison with the Landsat data? Is the temporal resolution of the Sentinel-2 data sufficient to detect the changes and the health of the forest within the year? 2. Materials and Methods 2.1. Study Areas 2.1.1. In Situ and Auxiliary Data Two types of study areas, with a significant and without any significant impact of a disturbance on the forest, were selected in the Sumava National Park (Sumava NP, Czechia) and the Low Tatras National Park (NAPANT, Slovakia). The study areas were selected based on: (1) the forest archive documents of the national parks, which provided data on the status of the forests and their management, the area and category of the forests, the average age of the stock and the predominant tree species representation; (2) the in situ research and field observations focused on the height, age, and density of the trees, the plant and species composition, the characteristics of the plant undergrowth, the visual observations of the conditions of the individual trees; and (3) the aerial photo archives from the Czech Office for Surveying, Mapping and Cadastre (CUZK; orthophoto for Czechia), the office of Geodesy, Cartography and Cadastre Authority of the Slovak Republic (UGKK SR; orthophoto for Slovakia), Mapy.cz (GEODIS, TopGis for Czechia and EUROSENSE, GEODIS for Slovakia), and Google Earth Pro (for both countries). The archive data and aerial photographs were helpful for the selection of the suitable study areas, the investigation of the forest status and changes as well as for the validation Remote Sens. 2020,12, 1914 5 of 26 of the results achieved by the TS analysis of the Sentinel-2 data. The field research was carried out repeatedly in the summer and autumn of the years 2017, 2018, and 2019. The study areas were localized by a global navigation satellite system receiver (Trimble Geoexplorer 6000 Geo XT GPS). 2.1.2. Study Area Description The study areas are located in the Sumava National Park (Sumava NP, Czechia) and the Low Tatras National Park (NAPANT, Slovakia). The Low Tatras are among the most important mountainous regions of Slovakia with a maximum height of 2043 m above sea level at the summit of Dumbier. The declaration of the Low Tatras National Park occurred in 1978. The national park covers an area of approximately 728 km 2 with a buffer zone of more than 1102 km 2 . NAPANT is one of the largest protected areas, not only in Slovakia, but in the entire Carpathians. Due to its large size, the diverse soil substrate, and the variation in the forms of relief, the Low Tatras are a part of the area with the highest diversity of plant species in Slovakia. The dominant plant community in the area is forest, with the presence of extensive forest ecosystems [ 28 ]. The Sumava mountains are an extensive mountain range on the border of Czechia, Austria, and Germany. The mean altitude is 928 m above sea level. In Czechia, since 1963, the territory has been managed as a protected landscape area, and as a national park since 1991. The area of the national park is approximately 681 km 2 [ 19 , 49 ]. The dominant land cover is a coniferous forest, which is under pressure from biotic, abiotic, and anthropogenic influences [ 14 , 15 ]. In total, five study areas of interest were selected (three in the Sumava NP and two in the NAPANT), where the field research was conducted (Table 1). All the study areas are depicted in the map (Figure 1). The study areas represent different areas according to the type of the forest vegetation degradation (areas affected by a disturbance and areas without any significant disturbances). The extent of all the study areas was defined as 60 × 60 m (3600 m 2 ). This extent ensures pure pixels (unmixed) of the individual forest type and it is compatible with the 10-/20-m bands of the Sentinel-2 (6 × 6/ 3×3 pixels). Table 1. Characteristics of the study areas. Source: author’s own work (2020). Study Area Country (National Park) Latitude (WGS84) Longitude (WGS84) Altitude (m) Disturbance Description SA1dist Czechia (Sumava National park) 13.45196 49.02482 1040 YES Bark beetle calamity (red-attack 2015, gray-attack prevailed 2016–2017). From 2018 in the recovery phase. SA2recov Czechia (Sumava National park) 13.52290 48.98473 1160 YES Recovery phase after bark beetle calamity. Small spruces and shrubs as undergrowth. SA3non-dist Czechia (Sumava National park) 13.47758 49.04568 1035 NO Stable ecosystem with spruce trees, grass and moss as undergrowth. Visible impact of drought in last years. SA4recov Slovakia (Low Tatras National Park) 19.61421 48.89590 1300 YES Recovery phase after bark beetle calamity. Small trees and shrubs as undergrowth. SA5non-dist Slovakia (Low Tatras National Park) 19.66245 48.95963 1335 NO Stable ecosystem with spruce trees. SA—study area; dist —disturbance, recov —recovery, non-dist —no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Remote Sens. 2020,12, 1914 6 of 26 RemoteSens.2020,12,xFORPEERREVIEW6of26 Table1.Characteristicsofthestudyareas.Source:author’sownwork(2020). Study Area Country (NationalPark) Latitude (WGS84) Longitude (WGS84) Altitude (m)DisturbanceDescription SA1dist Czechia (Sumava Nationalpark) 13.4519649.024821040YES Barkbeetlecalamity(red‐attack2015, gray‐attackprevailed2016–2017). From2018intherecoveryphase. SA2recov Czechia (Sumava Nationalpark) 13.522948.984731160YES Recoveryphaseafterbarkbeetle calamity.Smallsprucesandshrubs asundergrowth. SA3non‐dist Czechia (Sumava Nationalpark) 13.4775849.045681035NO Stableecosystemwithsprucetrees, grassandmossasundergrowth. Visibleimpactofdroughtinlast years. SA4recov Slovakia(Low TatrasNational Park) 19.6142148.89591300YES Recoveryphaseafterbarkbeetle calamity.Smalltreesandshrubsas undergrowth. SA5non‐dist Slovakia(Low TatrasNational Park) 19.6624548.959631335NOStableecosystemwithsprucetrees. SA—studyarea;dist—disturbance,recov—recovery,non‐dist—nodisturbance;areanumber(1–3from Czechia,4–5fromSlovakia).  Figure1.StudyareasoftheSumavaNationalPark(Czechia)andtheLowTatrasNationalPark (Slovakia).Source:author’sownwork(2020)/WebMapService(WMS)OrthophotobyCzechOffice forSurveying,MappingandCadastre(CUZK)(2019)/WMSZBGIS(Digitalmapapp)byGeodesy, CartographyandCadastreAuthorityoftheSlovakRepublic(UGKKSR)(2018). Figure 1. Study areas of the Sumava National Park (Czechia) and the Low Tatras National Park (Slovakia). Source: author’s own work (2020)/Web Map Service (WMS) Orthophoto by Czech Office for Surveying, Mapping and Cadastre (CUZK) (2019)/WMS ZBGIS (Digital map app) by Geodesy, Cartography and Cadastre Authority of the Slovak Republic (UGKK SR) (2018). The study area 1 (SA1 dist ) area was selected in the Sumava NP. In this area, a bark beetle attack occurred in 2015 (the red attack). A gray-attack prevailed in 2016–2017 and the study area has been in recovery phase since 2018 (Figure 2). The average altitude is 1040 m above sea level. The study area 2 (SA2 recov ) represents a place where a massive bark beetle attack occurred in the Sumava NP with the culmination around the year 2012 (the disturbance occurred between 2009–2012). Currently, the slow renewal of the ecosystem is taking place without any forest management intervention. More than 50% of the area is still covered by non-logged dead trees (spruces). Although this SA2 recov is a wetland area, the location has been under pressure from drought over the years since 2017. New small trees (spruces) and shrubs are visible in the territory (Figure 3). The average altitude is 1160 m above sea level. The study area 3 (SA3 non-dist ) located in the Sumava NP has not incurred any significant calamity or disturbance over the past 30 years. The altitude is approximately 1035 m above sea level (Figure 4). This is a forest ecosystem covered by grass and moss in the undergrowth that maintains a stable humidity. However, in recent years, we can see the impact of the drought on the forest. A reduction in the number of needles on the trees during the long-term field survey was observed. The study area 4 (SA4 recov ) is located in the NAPANT with an average altitude of 1300 m above sea level. This study area represents a place where the bark beetle outbreak gave rise to a decline in the spruce trees on the upper edge of the forest line level between the years 2007–2012 (Figure 5). In 2007, the wind calamity Kyrill took place, which caused, after 2007, a massive expansion of the bark beetle calamity in the surrounding areas. After the disturbance, the damaged trees were not logged and, currently, the recovery phase is still in progress. Remote Sens. 2020,12, 1914 7 of 26 RemoteSens.2020,12,xFORPEERREVIEW7of26  Figure2.Studyarea1(SA1 dist ,SumavaNP,Czechia)indetail:Withabarkbeetleoutbreak(aerial photocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMSOrthophotobyCZUK (2019).Theareaofinterestisintheredpolygon.  Figure3.Studyarea2(SA2 recov ,SumavaNP,Czechia)indetail:Witharecoveryphaseafterthebark beetleoutbreak(aerialphotocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMS OrthophotobyCUZK(2019).Theareaofinterestisintheredpolygon. Figure 2. Study area 1 (SA1 dist , Sumava NP, Czechia) in detail: with a bark beetle outbreak (aerial photo created by CUZK in 2019). Source: author’s own work (2020)/WMS Orthophoto by CZUK (2019). The area of interest is in the red polygon. RemoteSens.2020,12,xFORPEERREVIEW7of26  Figure2.Studyarea1(SA1 dist ,SumavaNP,Czechia)indetail:Withabarkbeetleoutbreak(aerial photocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMSOrthophotobyCZUK (2019).Theareaofinterestisintheredpolygon.  Figure3.Studyarea2(SA2 recov ,SumavaNP,Czechia)indetail:Witharecoveryphaseafterthebark beetleoutbreak(aerialphotocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMS OrthophotobyCUZK(2019).Theareaofinterestisintheredpolygon. Figure 3. Study area 2 (SA2 recov , Sumava NP, Czechia) in detail: with a recovery phase after the bark beetle outbreak (aerial photo created by CUZK in 2019). Source: author’s own work (2020)/WMS Orthophoto by CUZK (2019). The area of interest is in the red polygon. Remote Sens. 2020,12, 1914 8 of 26 RemoteSens.2020,12,xFORPEERREVIEW8of26  Figure4.Studyarea3(SA3 non‐dist ,SumavaNP,Czechia)indetail:Withminimaldisturbances(aerial photocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMSOrthophotobyCUZK (2019).Theareaofinterestisintheredpolygon.  Figure5.Studyarea4(SA4 recov ,theLowTatras,Slovakia)indetail:Witharecoveryphaseafterthe barkbeetleoutbreak(aerialphotocreatedbyUGKKSRin2018).Source:author’sownwork (2020)/WMSZBGISbyUGKKSR(2018).Theareaofinterestisintheredpolygon. Figure 4. Study area 3 (SA3 non-dist , Sumava NP, Czechia) in detail: with minimal disturbances (aerial photo created by CUZK in 2019). Source: author’s own work (2020)/WMS Orthophoto by CUZK (2019). The area of interest is in the red polygon. RemoteSens.2020,12,xFORPEERREVIEW8of26  Figure4.Studyarea3(SA3 non‐dist ,SumavaNP,Czechia)indetail:Withminimaldisturbances(aerial photocreatedbyCUZKin2019).Source:author’sownwork(2020)/WMSOrthophotobyCUZK (2019).Theareaofinterestisintheredpolygon.  Figure5.Studyarea4(SA4 recov ,theLowTatras,Slovakia)indetail:Witharecoveryphaseafterthe barkbeetleoutbreak(aerialphotocreatedbyUGKKSRin2018).Source:author’sownwork (2020)/WMSZBGISbyUGKKSR(2018).Theareaofinterestisintheredpolygon. Figure 5. Study area 4 (SA4 recov , the Low Tatras, Slovakia) in detail: with a recovery phase after the bark beetle outbreak (aerial photo created by UGKK SR in 2018). Source: author’s own work (2020)/WMS ZBGIS by UGKK SR (2018). The area of interest is in the red polygon. Remote Sens. 2020,12, 1914 9 of 26 The study area 5 (SA5 non-dist ) is located in the NAPANT and represents a place without any disturbance, with approximately 70-year-old trees clearly in good condition (Figure 6). The average altitude is 1335 m above sea level. The predominant representation of the forest trees is spruce (Picea abies), which naturally creates a functioning ecosystem of mountain spruce forests in the NAPANT. RemoteSens.2020,12,xFORPEERREVIEW9of26  Figure6.Studyarea5(SA5 non‐dist ,theLowTatras,Slovakia)indetail:Withminimaldisturbances (aerialphotocreatedbyUGKKSRin2018).Source:author’sownwork(2020)/WMSZGBISbyUGKK SR(2018).Theareaofinterestisintheredpolygon. 2.2.OpticalData 2.2.1.Sentinel‐2 Forthisstudy,theatmosphericallycorrectedSentinel‐2L2Adata(usingtheSen2Corcorrection) withoutsnowandcloudcoverwereused.Thecloudandsnowcoverwereidentifiedbyvisual interpretationoftheimagesusingtrueandfalsecolorRGBbandcombinations.Weworkedwith Sentinel‐2L2Adata(sensedbysatellites2Aand2B)with10‐ and20‐mspatialresolution.The Sentinel‐2datawereusedforthecalculationofthevegetationindicesandTSanalysis. 2.2.2.Landsat8 ThesecondsourceofsatellitedatawastheLandsat8imageswitha30‐mspatialresolution. Landsat8datawereusedfortheevaluationofthetemporalandspatialresolutionbasedona comparisonwiththeSentinel‐2. BoththeavailableSentinel‐2andLandsat8datawereretrievedfromtheSentinelHubarchive. ThecompletenessofthearchiveoftheSentinelHubapplicationwasvalidated/checkedusingthe CopernicusOpenAccessHub(forSentinel‐2)andusingtheUnitedStatesGeologicalSurvey(USGS) EarthExplorer(forLandsat8).Theobservationperiodwasfrom28March2017(releaseofSentinel‐ 2L2Adata)to31December2019.   Figure 6. Study area 5 (SA5 non-dist , the Low Tatras, Slovakia) in detail: with minimal disturbances (aerial photo created by UGKK SR in 2018). Source: author’s own work (2020)/WMS ZGBIS by UGKK SR (2018). The area of interest is in the red polygon. 2.2. Optical Data 2.2.1. Sentinel-2 For this study, the atmospherically corrected Sentinel-2 L2A data (using the Sen2Cor correction) without snow and cloud cover were used. The cloud and snow cover were identified by visual interpretation of the images using true and false color RGB band combinations. We worked with Sentinel-2 L2A data (sensed by satellites 2A and 2B) with 10and 20-m spatial resolution. The Sentinel-2 data were used for the calculation of the vegetation indices and TS analysis. 2.2.2. Landsat 8 The second source of satellite data was the Landsat 8 images with a 30-m spatial resolution. Landsat 8 data were used for the evaluation of the temporal and spatial resolution based on a comparison with the Sentinel-2. Both the available Sentinel-2 and Landsat 8 data were retrieved from the Sentinel Hub archive. The completeness of the archive of the Sentinel Hub application was validated/checked using the Copernicus Open Access Hub (for Sentinel-2) and using the United States Geological Survey (USGS) Earth Explorer (for Landsat 8). The observation period was from 28 March 2017 (release of Sentinel-2 L2A data) to 31 December 2019. Remote Sens. 2020,12, 1914 16 of 26 of the national parks). The values of SA1 dist were the lowest from all the study areas. This area was in the onset of the recovery phase. 3.2.3. Tasseled Cap Greenness (TCG) Time Series Figure 10c represents the values of the TCG orthogonal index. The chart shows minimal differences in the values of the study areas with recovery phase or healthy forest, so it is difficult to find any specific characteristic values related to the type of the study area. The differences between the disturbed area with an initial phase of recovery (SA1 dist ) and undisturbed areas or areas in the more advanced recovery phases are visible. The values of the disturbed forest vegetation with an initial recovery phase represent a range from 0.02 to 0.09 and the undisturbed or in the recovery phase range from 0.02 to 0.20. The SA1 dist under a strong impact of disturbance had much lower variability of the values than the areas in the advanced recovery phase (SA2recov and SA4recov). 3.2.4. Tasseled Cap Wetness (TCW) Time Series The TCW orthogonal index allowed us to define the specific groups of the values (Figure 10d). Differentiating the values of the disturbed and undisturbed study areas was possible. The values of the disturbed forest or recovery forest were concentrated in a range approximately from − 0.10 to 0.01 and the undisturbed study areas ranged from 0.00 to 0.03. The values higher than 0.00 could be mostly considered as an undisturbed forest and the lower values as a forest under disturbance or in the recovery phase. Concerning the disturbed and recovery areas (SA1 dist , SA2 recov , and SA4 recov ), it is possible to recognize the more apparent recovery phase of the forest in SA4 recov (NAPANT) characterized by the higher increasing the values of the TCW index. At the end of 2019, the values of SA4recov were higher than the values of SA1dist and SA2recov. 3.2.5. Statistical Analyses Figures 11 and 12 demonstrate the values of the indices and their differences based on the mathematical statistical boxplots of all the study areas for the whole observed period. The results show that the SA3 non-dist healthy forest both in the Sumava NP Czechia) or SA5 non-dist in the NAPANT (Slovakia) had very similar values in all the observed indices. The differences in the values between the disturbed (or in the recovery mode) and undisturbed areas were possible to document on the basis of the NDVI, NDMI, and TCW (higher values in SA3non-dist and SA5non-dist and lower values in SA1dist, SA2 recov , and SA4 recov ). On the other hand, very similar values were documented in the TCG index; only the SA1 dist had slightly lower values. The disturbed areas documented a lower average value of the NDVI, NDMI, and TCW indices. The variability of the values in the disturbed areas was much higher than in the areas without disturbance. The SA1 dist with the initial phase of recovery reflected a lower value of NDVI and NDMI than the SA2 recov and SA4 recov , which were in the more advanced recovery phase. The abilities of the TCG for the detection of undisturbed and disturbed vegetation in the different phases of recovery were weak and unsuitable for this purpose. Table 3shows the standard deviation of the selected indices. Table 4shows the p values of the Shapiro–Wilk test for the individual indices and for all the study areas. The values, where the null hypothesis holds (at α =0.05), are in italics and bold. None of the indices had a normal distribution for all the study areas. Based on these results, we decided to apply the Wilcoxon nonparametric test. Remote Sens. 2020,12, 1914 17 of 26 Remote Sens. 2020, 12, 1914 17 of 27 the disturbed (or in the recovery mode) and undisturbed areas were possible to document on the basis of the NDVI, NDMI, and TCW (higher values in SA3non-dist and SA5non-dist and lower values in SA1dist, SA2recov, and SA4recov). On the other hand, very similar values were documented in the TCG index; only the SA1dist had slightly lower values. The disturbed areas documented a lower average value of the NDVI, NDMI, and TCW indices. The variability of the values in the disturbed areas was much higher than in the areas without disturbance. The SA1dist with the initial phase of recovery reflected a lower value of NDVI and NDMI than the SA2recov and SA4recov, which were in the more advanced recovery phase. The abilities of the TCG for the detection of undisturbed and disturbed vegetation in the different phases of recovery were weak and unsuitable for this purpose. Table 3 shows the standard deviation of the selected indices. Figure 11. Study areas statistics. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). Boxplots are divided into two parts by median. The cross indicates the arithmetic mean. Dots indicate outliers. TCG—tasseled cap greenness; TCW—tasseled cap wetness; SA—study area; dist—disturbance, recov—recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). SA(d)1 SA(d)2 SA(d)3 SA(d)4 SA(d)5 Index Value -0.2 0 0.2 0.4 0.6 0.8 1 NDVI NDMI TCG TCW SA3 non-dist SA4 recov SA5 non-dist SA1 dist SA2 recov Figure 11. Study areas statistics. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). Boxplots are divided into two parts by median. The cross indicates the arithmetic mean. Dots indicate outliers. TCG—tasseled cap greenness; TCW—tasseled cap wetness; SA—study area; dist —disturbance, recov —recovery, non-dist —no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Remote Sens. 2020, 12, 1914 18 of 27 Figure 12. NDVI, NDMI, TCG, and TCW statistics. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). Boxplots are divided into two parts by the median. The cross indicates the arithmetic mean. Dots indicate outliers. SA—study area; dist—disturbance, recov— recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Table 3. Comparison of the standard deviation for all the vegetation indices/study areas. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). SA1dist SA2recov SA3non-dist SA4recov SA5non-dist NDVI 0.057 0.070 0.044 0.140 0.076 NDMI 0.066 0.062 0.041 0.075 0.064 TCG 0.017 0.038 0.012 0.078 0.027 TCW 0.015 0.012 0.006 0.019 0.010 SA—study area; dist—disturbance, recov—recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Table 4 shows the p values of the Shapiro–Wilk test for the individual indices and for all the study areas. The values, where the null hypothesis holds (at α = 0.05), are in italics and bold. None of the indices had a normal distribution for all the study areas. Based on these results, we decided to apply the Wilcoxon nonparametric test. Table 4. The Shapiro–Wilk test for all study areas and vegetation indices. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). SA1dist SA2recov SA3non-dist SA4recov SA5non-dist NDVI 0.001383 0.023480 0.719900 0.000001 0.003258 NDMI 0.789300 0.044700 0.823200 0.002179 0.378700 TCG 0.773300 0.046910 0.004701 0.000001 0.000703 TCW 0.010030 0.128700 0.652100 0.131400 0.013750 Values in italics and bold have a normal distribution, the non-italic and non-bold (regular) values do not. SA—study area; dist—disturbance, recov—recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Table 5 shows the results of the Wilcoxon test for each study area combination and for each vegetation index. The values of the combinations, for which no statistically significant difference was found, are in italics and bold. The results of the Wilcoxon test show that the NDMI index seems the most appropriate for distinguishing the recovery and healthy vegetation. The Wilcoxon test indicates NDVI NDMI TCG TCW Index Value -0.2 0 0.2 0.4 0.6 0.8 1 SA1dist SA2recov SA3non-dist SA4recov SA5non-dist SA1 dist SA2 recov SA3 non-dist SA4 recov SA5 non-dist Figure 12. NDVI, NDMI, TCG, and TCW statistics. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). Boxplots are divided into two parts by the median. The cross indicates the arithmetic mean. Dots indicate outliers. SA—study area; dist —disturbance, recov—recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Table 3. Comparison of the standard deviation for all the vegetation indices/study areas. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). SA1dist SA2recov SA3non-dist SA4recov SA5non-dist NDVI 0.057 0.070 0.044 0.140 0.076 NDMI 0.066 0.062 0.041 0.075 0.064 TCG 0.017 0.038 0.012 0.078 0.027 TCW 0.015 0.012 0.006 0.019 0.010 SA—study area; dist —disturbance, recov —recovery, non-dist —no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Remote Sens. 2020,12, 1914 18 of 26 Table 4. The Shapiro–Wilk test for all study areas and vegetation indices. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). SA1dist SA2recov SA3non-dist SA4recov SA5non-dist NDVI 0.001383 0.023480 0.719900 0.000001 0.003258 NDMI 0.789300 0.044700 0.823200 0.002179 0.378700 TCG 0.773300 0.046910 0.004701 0.000001 0.000703 TCW 0.010030 0.128700 0.652100 0.131400 0.013750 Values in italics and bold have a normal distribution, the non-italic and non-bold (regular) values do not. SA—study area; dist —disturbance, recov —recovery, non-dist —no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). Table 5shows the results of the Wilcoxon test for each study area combination and for each vegetation index. The values of the combinations, for which no statistically significant difference was found, are in italics and bold. The results of the Wilcoxon test show that the NDMI index seems the most appropriate for distinguishing the recovery and healthy vegetation. The Wilcoxon test indicates similarities between the SA1 dist and SA2 recov combination or SA3 non-dist and SA5 non-dist combination for the TCW index and thus confirmed the distinction between disturbed (or recovery mode in the initial phase, not the more advanced phase like SA4 recov ) and undisturbed areas. For the NDVI index, similarity was shown only for SA3 non-dist and SA5 non-dist combination; other study areas reflected different values. The test also proved the unsuitability of the TCG index for distinguishing individual forest phases. The results may be affected by multiple testing. Table 5. The Wilcoxon test for selected vegetation indices with all combinations of study areas. Source: author’s own work (2020)/Sentinel Hub and ESA Sentinel-2 data product (2020). NDVI NDMI SA1dist SA2recov SA3non-dist SA4recov SA1dist SA2recov SA3non-dist SA4recov SA2recov 2.80 ×10−7SA2recov 2.00 ×10−9 SA3non-dist 1.60 ×10−15 4.80 ×10−15 SA3non-dist <2×10−16 6.10 ×10−16 SA4recov 1.60 ×10−90.023 3.40 ×10−8SA4recov 1.10 ×10−90.510 <2×10−16 SA5non-dist 1.30 ×10−13 3.90 ×10−13 0.377 8.50 ×10−10 SA5non-dist 6.30 ×10−14 1.70 ×10−13 0.190 8.40 ×10−15 TCG TCW SA1dist SA2recov SA3non-dist SA4recov SA1dist SA2recov SA3non-dist SA4recov SA2recov 2.20 ×10−11 SA2recov 0.371 SA3non-dist 2.30 ×10−16 0.972 SA3non-dist 2.30 ×10−16 6.90 ×10−16 SA4recov 0.002 0.182 0.346 SA4recov 3.20 ×10−51.20 ×10−41.30 ×10−15 SA5non-dist 5.30 ×10−11 0.972 0.133 0.098 SA5non-dist 6.30 ×10−14 1.70 ×10−13 0.075 1.90 ×10−14 Values in italics and bold show significant dependence between the area combination, otherwise the non-italic and non-bold (regular) values show independence. SA—study area; dist —disturbance, recov —recovery, non-dist—no disturbance; area number (1–3 from Czechia, 4–5 from Slovakia). 4. Discussion Disturbance phenomena, such bark beetle and wind calamities, are one of the most urgent problems in the forest ecosystems of Central Europe. Their overgrowth is closely related to climate changes and the spread of non-indigenous coniferous trees in Central Europe, e.g., spruce (Picea abies). In both observed national parks (the Sumava National Park in Czechia and the NAPANT in Slovakia), spruce is the dominant tree. The dominance of spruce monocultures was caused by foresters due to market-oriented forest management. Currently, under new environmental politics and under the influence of the climate changes, national parks have implemented more environmentally oriented and sustainable management practices, often without any interventions of the foresters. The primary task of that management is to achieve a natural and environmentally stable forest. Spruce monocultures, higher temperatures and abiotic calamities can cause the appropriate conditions for bark beetle calamities. On the other hand, a sustainable management should create a natural ecosystem suitable for richer communities of species compared to the spruce monoculture. Different approaches are now visible in forest management in Europe. The forests in the NAPANT national park are mostly forestry managed. In the case of the occurrence of a bark beetle outbreak, Remote Sens. 2020,12, 1914 19 of 26 the obligation of the owner/administrator to cut down the forest in the intervention zone exists. Typically, more than half of the dead trees are harvested after the calamity. However, about 15% of the national park area is currently in a non-intervention mode (SA4 recov ). In the case of a bark beetle occurrence, an obligation of the owner/administrator to cut down the forest in intervention zone exists, despite the value of the nonprofit functions in the nature–close management of forests. On the other hand, a forest management practice with very limited intervention was implemented in the Sumana NP. Hence, a natural renewal of the forest ecosystem has been taking place in the Sumava NP on a prevailing area. New small trees and shrubs are visible in the recovery territories (SA2recov). In this study, we applied and tested the Sentinel-2 data for the evaluation and monitoring of the forest ecosystems in Czechia (the Sumava National Park) and Slovakia (the NAPANT). Due to the dynamic changes in the forests in Czechia and Slovakia that occurred over the last years, there exists a suitable opportunity to apply the Sentinel-2 TS for the detection of disturbance events (bark beetle) and to evaluate the forest health (forest condition before, during, and after the disturbance events). For this reason, the study areas represent a forest under different circumstances and development (the dead trees after a bark beetle outbreak, the recovery mode after a bark beetle outbreak, and the forest vegetation without any significant influence). The selected vegetation indices and their trajectories of the TS were interpreted and validated in relation to the in situ data investigated during the field research or provided by the administration of the national parks and to the aerial photographs. The Copernicus Program with the Sentinel-2 data has brought new opportunities and perspectives in the landscape/forest monitoring. Thanks to the temporal and spatial resolution, we were able to obtain a short TS for observing the health status and changes in the forest ecosystems [ 59 ]. Due to the advanced parameters of the Sentinel mission (e.g., two satellites—Sentinel A and Sentinel B), the Sentinel-2 data allowed us to observe the status and changes of the vegetation repeatedly in a short time period (5 days or less) with multispectral characteristics in the data (13 bands). The Sentinel’s TS with the higher temporal and spatial resolution were undoubtedly one of the greatest benefits of studying regional or even global phenomena [ 60 ]. Other authors [ 61 ] claimed that the probability of acquired cloud-free pixels by the Landsat 8 during a summer period in the U.S. was 0.78. This fact often led researchers to study the changes in the forest over several years rather than during the seasons of the year (phenophases). We proved, in this study, that the Sentinel-2 data had a significantly higher density in the TS when compared with the Landsat 8 data. Clearly, many more images of the Sentinel-2 (2A and 2B), compared with the Landsat images, were available for all the areas in 2018 and 2019. Due to the launch of the Sentinel-2B, the temporal resolution rapidly increased. This resolution is a very important factor for the evaluation of the forest with dynamic changes (e.g., under a disturbance aspect) in mountainous areas where cloud cover is a relevant problem. The Sentinel-2 data brings a new dimension in the data availability and its improved temporal resolution using two satellites (2A and 2B) was more suitable for the monitoring of the forest dynamic than the traditionally used Landsat data. Concerning the processing of the data with a lower temporal resolution (such as Landsat), several preprocessing steps for preparing the datasets are necessary. To obtain comparable results, it is necessary to use normalization [ 20 , 62 , 63 ], cross-calibration methods, specialized harmonizing algorithms (like LandsatLinkr [ 29 ]), or to use other methods (like LandTrendr [ 64 , 65 ]). However, in the case of the Sentinel-2 data, normalization methods were not necessary for processing for the TS when using the entire set of measurements with a high temporal resolution (outliers caused by inhomogeneity of the environment during the year may be filtered out—or the data may be fit). The Sentinel-2 data brought much sharper images, especially due to the 10-m resolution bands. This spatial resolution allowed us to better recognize individual elements, such as damaged forests, shadows, clouds, and other covers. However, due to the low spatial resolution of Landsat, it was difficult to detect cirrus clouds. For the Landsat 8 data, it was possible to use panchromatic images with a 15-m resolution and to perform pan-sharpening methods to obtain a better spatial resolution. However, this was only possible for visible bands and this approach required complex image preprocessing. Remote Sens. 2020,12, 1914 20 of 26 The cloud-based Sentinel Hub application was used for the processing and analysis of the Sentinel-2 data. The cloud-based tools allowed for the effective and end-user friendly processing and analysis of the EO data. An advantage of the Sentinel Hub was the option to prepare and use scripts in the Configuration Utility for the semi-automatic processing of the data and the calculation of the indices from the Sentinel-2 data (NDVI, NDMI, tasseled cap greenness, and tasseled cap wetness). The Sentinel Hub also contained prepared algorithms for the calculation of the most used indices. Performing land use/land cover classifications and other image analyses was also possible. Several programming scripts using JavaScript language were written for the calculation of the vegetation indices in this study. The impact of the disturbances on the forest vegetation was evaluated using the NDVI, NDMI, TCG, and TCW indices. The trends (similarities and differences) in the values of the studied indices were analyzed using the TS methods. We found that the vegetation indices NDVI, NDMI, and TCW were able to distinguish healthy and disturbed/damaged forests and they gave most relevant results for the research purposes of this study. The reflectivity of the healthy forest vegetation was higher in the NIR band than the SWIR; however, the SWIR reflectivity was higher in the case of disturbances. This aspect played an important role as the SWIR bands responded sensitively in the case of the degradation of the forest. The NDMI index using the SWIR band was evidently able to precisely distinguish the individual types of areas affected by disturbances (without disturbance and under disturbance or after disturbance in the individual phases of recovery). The values of the disturbed and recovery areas did not overlap with the undisturbed areas and create specific groups in the TS chart. The values of the NDMI for the unaffected areas ranged approximately from 0.37 to 0.69, and the areas affected by the disturbances or in a recovery phase had evidently lower values of the NDMI with a range from − 0.12 to 0.28. The orthogonal TCW index had a similar ability; however, it was not able to precisely distinguish the disturbed and recovery phases/areas. The values of the TCW index were documented from − 0.10 to 0.00 in the disturbed and recovery areas and from 0.00 to 0.03 for the unaffected areas. The NDVI index, which is based on the visible band and the NIR ratios, was able to reflect an evolution of the plants and shrubs that quickly covered the disturbed area during the recovery phase (see Figure 10). Concerning the TCG index, the abilities of this index did allow us to detect the disturbances in the gray-attack phase. However, it was difficult to distinguish the different recovery phases of the study areas. Our results are in accordance with comparable studies [ 66 , 67 ]. The novelty of this study was to use and test the Sentinel-2 data in the TS of the dynamic forest changes; therefore, we did not find any similar studies with comparable results of the value indices in the TS analysis based on Sentinel-2 data. An important finding gained in this study was the high relevance of the data of the forest management practices in the national parks. The data contained in the forest management plans appeared to be excellent as a source for the validation of the results of the EO/TS investigation. The forest management plans contained a wide range of useful information, and this information is regularly updated. The data from the in situ research provided highly detailed data on the recovery of the forest ecosystems, especially with regards to the restoration of the herbaceous and woody undergrowth. The use of the Sentinel-2 data in the evaluation of the forest vegetation and the consequences of the differentiated approaches to conservation and forestry activities provided excellent capacities and opportunities for the monitoring of the forest and the implementation of the EO data and the TS method into the protection and management of the forest. In addition, the new Landsat 9 (with similar OLI and TIRS sensor) will be sent into orbit soon (March 2021). If the Landsat 8 is still working in the following years, the new temporal coverage of the Landsat images will be higher (similar to Sentinel-2A and Sentinel-2B). The harmonization of the Landsat 8 (or Landsat 9) and Sentinel-2 data could also increase the temporal coverage of the TS [ 68 ]. The possibilities to increase the temporal and spectral resolution can be seen in the data fusion of optical and SAR data also (e.g., Sentinel-1 and Sentinel-2). From the method point of view, it would be useful to test an application of the cloud-based Google Earth Engine (GEE), which allows for the analogical analyses of the EO data, such as the Sentinel Hub. Remote Sens. 2020,12, 1914 21 of 26 Another idea is to carry out a comparative study in other national parks in Europe and validate our achieved results in the different forest ecosystems. 5. Conclusions In our study, we evaluated the changes in forest vegetation affected by the bark beetle in selected national parks in Czechia and Slovakia (the Sumava NP and the NAPANT). The study used the open data of Sentinel-2 from the Copernicus program. The TS analysis using the selected vegetation indices was performed in the areas in both countries. Based on our achieved results, we argue that the Sentinel-2 data were able to accurately distinguish the areas that were not affected by disturbances and the areas under disturbances and to detect the individual phases of the recovery mode of the forest vegetation. The results confirm the abilities of the NDVI, NDMI, and TCW indices to distinguish disturbed and undisturbed areas. The NDMI vegetation index was useful for the detection of the disturbed forest and forest recovery after bark beetle outbreaks and provided relevant information regarding the health of the forest (the individual stages of the recovery mode). On the contrary, the TCG index demonstrated limited abilities. The TCG could distinguish the gray-attack disturbance phase; however, it was difficult to use this index for the detection of different recovery phases and to distinguish recovery phases from healthy forests. The data were processed in the cloud-based tool of the Sentinel Hub, which allowed for very flexible data processing and analyzing. A combination of a cloud-based easy operated system with the open data of Sentinel-2 with a high spatial, temporal, and spectral resolution provided a powerful tool for forest research. This study should serve as an example for the application of advanced remote sensing methods and as data if the impact of the disturbances on the forest vegetation is analyzed using vegetation indices. These results and methods should be useful and inspirational for forest managers to determine appropriate forest management practices under the circumstances of forest disasters. Author Contributions: Conceptualization, J.L. and P.S. (Premysl Stych); methodology, J.L. and P.S. (Premysl Stych); software, J.L. and N.K.; validation, D.P., R.H., J.L., P.S. (Premysl Stych) and P.S. (Pavel Svec); formal analysis, J.L.; investigation, J.L. and P.S. (Premysl Stych); resources, J.L.; data curation, N.K., J.S. and R.H.; writing—original draft preparation, J.L.; writing—review and editing, P.S. (Pavel Svec) and P.S. (Premysl Stych); visualization, J.L., D.P., R.H., N.K. and J.S.; supervision, P.S. (Pavel Svec) and P.S. (Premysl Stych); funding acquisition, J.L. and P.S. (Premysl Stych); senior author, P.S. (Premysl Stych). All authors have read and agreed to the published version of the manuscript. Funding: This study was supported by the Charles University, project GA UK No. 512217 (Charles Univ, Fac Sci, 2017–2019): “Hodnoceni vlivu disturbanci na lesni ekosystemy v Cesku a na Slovensku pomoci metod DPZ”. The UNCE “Charles University Research Center program UNCE/HUM/018”. The European Union’s Caroline Herschel Framework Partnership Agreement on Copernicus User Uptake under grant agreement No. FPA 275/G/GRO/COPE/17/10042, project FPCUP (Framework Partnership Agreement on Copernicus User Uptake), Action 2019-2-49 “Developing supports for monitoring and reporting of GHG emissions and removals from land use, land use change and forestry” (219/SI2.818795/07 (CLIMA)). The study was invited (Feature Paper) to be part of the Special Issue of Monitoring Forest Change with Remote Sensing. Acknowledgments: We would like to thank the Grant Agency of the Charles University (GA UK), project GA UK No. 512217 for their support as well as to thank the UNCE “Charles University Research Center program UNCE/HUM/018” for their support. We would like to thank the European Union’s Caroline Herschel Framework Partnership Agreement on Copernicus User Uptake under grant agreement No. FPA 275/G/GRO/COPE/17/10042, project FPCUP (Framework Partnership Agreement on Copernicus User Uptake), Action 2019-2-49 “Developing supports for monitoring and reporting of GHG emissions and removals from land use, land use change and forestry” (219/SI2.818795/07 (CLIMA)) for their support. We thank the Sentinel Hub team for using the Sentinel Hub services and for helping with this article, especially Grega Milcinski. We thank Zdenek Hlavka from the Faculty of Mathematics and Physics, Charles University for consultations on the statistical methods and interpretations of the results. We also thank the anonymous reviewers for their critical review and constructive comments. Conflicts of Interest: The authors declare no conflict of interest. Remote Sens. 2020,12, 1914 22 of 26 Appendix A Table A1. Sentinel-2 and Landsat 4–8 spectral bands and their differences. Mission Radiometric Resolution Temporal Resolution Spectral Band Coastal/Aerosol (B01) Blue (B02) Green (B03) Red (B04) NIR 2 (B05) NIR 3 (B06) NIR 4 (B07) NIR 1 (B08) NIR 5 (B8A) Atmospheric/Water Vapour (B09) Cirrus (B10) SWIR 1 (B11) SWIR 2 (B12) Panchromatic Thermal Infrared 1 Thermal Infrared 2 Date Landsat 4 and 5 8bit (16bit rescaled) 16 days Central Wavelength (nm) - 485 560 660 - - - 830 - - - 1650 2215 - 11450 - 16 July 1982 (deactivated: 15 June 2001)/ 1 March 1984 (deactivated: 5 June 2013) Bandwidth (nm) - 70 80 60 - - - 140 - - - 200 270 - 2100 - Spatial Resolution (m) - 30 30 30 - - - 30 - - - 30 30 - 120 - Landsat 7 Central Wavelength (nm) - 485 560 660 - - - 835 - - - 1650 2215 710 11450 - 15 April 1999 (active, data failure: 31 May 2003) Bandwidth (nm) - 70 80 60 - - - 130 - - - 200 270 380 2100 - Spatial Resolution (m) - 30 30 30 - - - 30 - - - 30 30 15 60 - Landsat 8 12bit (16bit rescaled) Central Wavelength (nm) 440 480 560 655 - - - 865 - - 1370 1610 2200 590 10895 12005 11 February 2013 (active) Bandwidth (nm) 20 60 60 30 - - - 30 - - 20 80 180 180 590 1010 Spatial Resolution (m) 30 30 30 30 - - - 30 - - 30 30 30 15 100 100 Sentinel-2A 12-bit 10-day (5-day with both satellites) Central Wavelength (nm) 443.9 496.6 560 664.5 703.9 740.2 782.5 835.1 864.8 945 1373.5 1613.7 2202.4 - - - 23 June 2015 (active) Bandwidth (nm) 27 98 45 38 19 18 28 145 33 26 75 143 242 - - - Spatial Resolution (m) 60 10 10 10 20 20 20 10 20 60 60 20 20 - - - Sentinel-2B Central Wavelength (nm) 442.3 492.1 559 665 703.8 739.1 779.7 833 864 943.2 1376.9 1610.4 2185.7 - - - 7 March 2017 (active) Bandwidth (nm) 45 98 46 39 20 18 28 133 32 27 76 141 238 - - - Spatial Resolution (m) 60 10 10 10 20 20 20 10 20 60 60 20 20 - - - Source: author’s own work (2020)/ESA (2015): Resolution and swath. Downloaded from: https://sentinel.esa.int/web/sentinel/missions/sentinel-2/instrument-payload/Resolution-andswath/USGS (2020): Spectral Response of the Operational Land Imager In-Band, Band-Average Relative Spectral Response. Downloaded from: https://landsat.gsfc.nasa.gov/preliminaryspectral-response-of-the-operational-land-imager-in-band-band-average-relative-spectral-response. The NIR 2-4 are the Vegetation Red Edge bands. Remote Sens. 2020,12, 1914 23 of 26 References 1. 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. [CrossRef] 2. Abdullah, H.; Skidmore, A.K.; Darvishzadeh, R.; Heurich, M. Sentinel-2 accurately maps green-attack stage of European spruce bark beetle (Ips typographus, L.) compared with Landsat-8. Remote Sens. Ecol. Conserv. 2019,5, 87–106. [CrossRef] 3. Hansen, M.C.; Stehman, S.V.; Potapov, P.V. Quantification of global gross forest cover loss. Proc. Natl. Acad. Sci. USA 2010,107, 8650–8655. [CrossRef] [PubMed] 4. Potapov, P.V.; Turubanova, S.A.; Tyukavina, A.; Krylov, A.M.; McCarty, J.L.; Radeloff, V.C.; Hansen, M.C. Eastern Europe’s forest cover dynamics from 1985 to 2012 quantified from the full Landsat archive. Remote Sens. Environ. 2015,159, 28–43. [CrossRef] 5. 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. [CrossRef] 6. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 2012,118, 83–94. [CrossRef] 7. Zhu, Z.; Woodcock, C.E. Continuous change detection and classification of land cover using all available Landsat data. Remote Sens. Environ. 2014,144, 152–171. [CrossRef] 8. Zhu, Z.; Woodcock, C.E.; Holden, C.; Yang, Z. Generating synthetic Landsat images based on all available Landsat data: Predicting Landsat surface reflectance at any given time. Remote Sens. Environ. 2015 ,162, 67–83. [CrossRef] 9. Neigh, C.S.R.; Bolton, D.K.; Diabate, M.; Williams, J.J.; Carvalhais, N. An Automated Approach to Map the History of Forest Disturbance from Insect Mortality and Harvest with Landsat Time-Series Data. Remote Sens. 2014,6, 2782–2808. [CrossRef] 10. DeRose, R.J.; Long, J.N.; Ramsey, R.D. Combining dendrochronological data and the disturbance index to assess Engelmann spruce mortality caused by a spruce beetle outbreak in southern Utah, USA. Remote Sens. Environ. 2011,115, 2342–2349. [CrossRef] 11. Hicke, J.A.; Logan, J. Mapping whitebark pine mortality caused by a mountain pine beetle outbreak with high spatial resolution satellite imagery. Int. J. Remote Sens. 2009,30, 4427–4441. [CrossRef] 12. 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. [CrossRef] 13. Griffiths, P.; Kuemmerle, T.; Baumann, M.; Radeloff, V.C.; Abrudan, I.V.; Lieskovsky, J.; Munteanu, C.; Ostapowicz, K.; Hostert, P. Forest disturbances, forest recovery, and changes in forest types across the Carpathian ecoregion from 1985 to 2010 based on Landsat image composites. Remote Sens. Environ. 2014 , 151, 72–88. [CrossRef] 14. Kozak, J.; Estreguil, C.; Troll, M. Forest cover changes in the northern Carpathians in the 20th century: A slow transition. J. Land. Use Sci. 2007,2, 127–146. [CrossRef] 15. Kuemmerle, T.; Hostert, P.; Radeloff, V.C.; Perzanowski, K.; Kruhlov, I. Post-Socialist Forest Disturbance in the Carpathian Border Region of Poland, Slovakia, and Ukraine. Ecol. Appl. 2007 ,17, 1279–1295. [CrossRef] [PubMed] 16. Butsic, V.; Munteanu, C.; Griffiths, P.; Knorn, J.; Radeloff, V.C.; Lieskovsk ý , J.; Mueller, D.; Kuemmerle, T. The effect of protected areas on forest disturbance in the Carpathian Mountains 1985–2010. Conserv. Biol. 2017,31, 570–580. [CrossRef] [PubMed] 17. Kupkova, L.; Potuckova, M.; Lhotakova, Z.; Albrechtova, J. Forest cover and disturbance changes, and their driving forces: A case study in the Ore Mountains, Czechia, heavily affected by anthropogenic acidic pollution in the second half of the 20th century. Environ. Res. Lett. 2018,13. [CrossRef] 18. Hais, M.; Wild, J.; Berec, L.; Bruna, J.; Kennedy, R.; Braaten, J.; Broz, Z. Landsat imagery spectral trajectories-important variables for spatially predicting the risks of bark beetle disturbance. Remote Sens. 2016,8, 687. [CrossRef] 19. Hais, M.; Jonasov á , M.; Langhammer, J.; Kucera, 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. [CrossRef] Remote Sens. 2020,12, 1914 24 of 26 20. Stych, P.; Lastovicka, J.; Hladky, R.; Paluba, D. 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. ISPRS Int. J. Geo-Inf. 2019,8, 71. [CrossRef] 21. Stych, P.; Jerabkova, B.; Lastovicka, J.; Riedl, M.; Paluba, D. A Comparison of Worldview-2 and Landsat 8 Images for the Classification of Forests Affected by Bark Beetle Outbreaks Using a Support Vector Machine and a Neural Network: A Case Study in the Sumava Mountains. Geosciences 2019,9, 396. [CrossRef] 22. Jensen, J.R. Remote Sensing of the Environment: An Earth Resource Perspective, 2nd ed.; Pearson Prentice Hall: Upper Saddle River, NJ, USA, 2007. 23. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the Great Plains with ERTS. In Proceedings of the 3rd ERTS Symposium, Washington, DC, USA, 10–14 December 1973; NASA SP-351 I. pp. 309–317. 24. 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. [CrossRef] 25. Gao, B. NDWI—A normalized difference water index for R. Elsevier 1996,266, 257–266. [CrossRef] 26. 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. [CrossRef] 27. Deering, D.; Rouse, J. Measuring ‘Forage Production’ of Grazing Units from Landsat MSS Data. In Proceedings of the 10th International Symposium on Remote Sensing of Environment, ERIM, Ann Arbor, MI, USA, 6–10 October 1975; pp. 1169–1178. 28. Hladky, R.; Lastovicka, J.; Holman, L.; Stych, P. Evaluation of the influence of disturbances on forest vegetation using Landsat time series; a case study of the Low Tatras National Park. Eur. J. Remote Sens. 2020 , 53, 40–66. [CrossRef] 29. Vogeler, J.C.; Braaten, J.D.; Slesak, R.A.; Falkowski, M.J. Extracting the full value of the Landsat archive: Inter-sensor harmonization for the mapping of Minnesota forest canopy cover (1973–2015). Remote Sens. Environ. 2018,209, 363–374. [CrossRef] 30. Wulder, M.A.; White, J.C.; Loveland, T.R.; Woodcock, C.E.; Belward, A.S.; Cohen, W.B.; Fosnight, E.A.; Shaw, J.; Masek, J.G.; Roy, D.P. The global Landsat archive: Status, consolidation, and direction. Remote Sens. Environ. 2016,185, 271–283. [CrossRef] 31. Franklin, S.E.; Wulder, M.A.; Skakun, R.S.; Carroll, A.L. Mountain pine beetle red-attack forest damage classification using stratified Landsat TM data in British Columbia, Canada. Photogramm. Eng. Remote Sens. 2003,69, 283–288. 32. Li, F.; Jupp, D.L.B.; Paget, M.; Briggs, P.R.; Thankappan, M.; Lewis, A.; Held, A. Improving BRDF normalisation for Landsat data using statistical relationships between MODIS BRDF shape and vegetation structure in the Australian continent. Remote Sens. Environ. 2017,195, 275–296. [CrossRef] 33. Roy, D.P.; Kovalskyy, V.; Zhang, H.K.; Vermote, E.F.; Yan, L.; Kumar, S.S.; Egorov, A. Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sens. Environ. 2016,185, 57–70. [CrossRef] 34. Immitzer, M.; Vuolo, F.; Atzberger, C. First Experience with Sentinel-2 Data for Crop and Tree Species Classifications in Central Europe. Remote Sens. 2016,8, 166. [CrossRef] 35. Rapinel, S.; Mony, C.; Lecoq, L.; Cl é ment, B.; Thomas, A.; Hubert-Moy, L. Evaluation of Sentinel-2 time-series for mapping floodplain grassland plant communities. Remote Sens. Environ. 2019,223, 115–129. [CrossRef] 36. Lambert, M.J.; Traor é , P.C.S.; Blaes, X.; Baret, P.; Defourny, P. Estimating smallholder crops production at village level from Sentinel-2 time series in Mali’s cotton belt. Remote Sens. Environ. 2018 ,216, 647–657. [CrossRef] 37. Piedelobo, L.; Hern á ndez-L ó pez, D.; Ballesteros, R.; Chakhar, A.; Del Pozo, S.; Gonz á lez-Aguilera, D.; Moreno, M.A. Scalable pixel-based crop classification combining Sentinel-2 and Landsat-8 data time series: Case study of the Duero river basin. Agric. Syst. 2019,171, 36–50. [CrossRef] 38. Ottosen, T.B.; Lommen, S.T.E.; Skjøth, C.A. Remote sensing of cropping practice in Northern Italy using time-series from Sentinel-2. Comput. Electron. Agric. 2019,157, 232–238. [CrossRef] 39. Rozenstein, O.; Haymann, N.; Kaplan, G.; Tanny, J. Estimating cotton water consumption using a time series of Sentinel-2 imagery. Agric. Water Manag. 2018,207, 44–52. [CrossRef] Remote Sens. 2020,12, 1914 25 of 26 40. Shang, R.; Zhu, Z. Harmonizing Landsat 8 and Sentinel-2: A time-series-based reflectance adjustment approach. Remote Sens. Environ. 2019,235, 111439. [CrossRef] 41. Masek, J.; Claverie, M.; Dungan, J.; Ju, J. Harmonized Landsat-8 and Sentinel-2. In User Guide; 2017. Available online: https://hls.gsfc.nasa.gov/wp-content/uploads/2017/08/HLS.v1.3.UserGuide_v2.pdf (accessed on 8 June 2020). 42. Shao, Z.; Cai, J.; Fu, P.; Hu, L.; Liu, T. Deep learning-based fusion of Landsat-8 and Sentinel-2 images for a harmonized surface reflectance product. Remote Sens. Environ. 2019,235, 111425. [CrossRef] 43. Navarro, A.; Catalao, J.; Calvao, J. Assessing the Use of Sentinel-2 Time Series Data for Monitoring Cork Oak Decline in Portugal. Remote Sens. 2019,11, 2515. [CrossRef] 44. Li, L.; Li, N.; Lu, D.; Chen, Y. Mapping Moso bamboo forest and its on-year and off-year distribution in a subtropical region using time-series Sentinel-2 and Landsat 8 data. Remote Sens. Environ. 2019 ,231, 111265. [CrossRef] 45. Grabska, E.; Hostert, P.; Pflugmacher, D.; Ostapowicz, K. Forest Stand Species Mapping Using the Sentinel-2 Time Series. Remote Sens. 2019,11, 1197. [CrossRef] 46. Gomes, V.C.F.; Queiroz, G.R.; Ferreira, K.R. An Overview of Platforms for Big Earth Observation Data Management and Analysis. Remote Sens. 2020,12, 1253. [CrossRef] 47. Sovdat, B.; Kadunc, M.; Batic, M.; Milcinski, G.; Functions, P.L. Natural color representation of Sentinel-2 data. Remote Sens. Environ. 2017,225, 392–402. [CrossRef] 48. Albrechtova, J.; Kupkova, L.; Campbell, P.K.E. Hodnoceni Stavu Smrkovych Porostu.-Pripadove Studie Sledovani Vyvoje Fyziologickeho Stavu Smrkovych Porostu v Krusnych Horach v Letech 1998–2013; Geographicaˇ Cesk á GeografickáSpoleˇcnost: NovéMˇesto, Czechia, 2017; ISBN ISBN 978-80-905642-9-9. 49. Hais, M.; Kucera, T. The influence of topography on the forest surface temperature retrieved from Landsat TM, ETM +and ASTER thermal channels. ISPRS J. Photogramm. Remote Sens. 2009 ,64, 585–591. [CrossRef] 50. Tang, D.; Fan, H.; Yang, K.; Zhang, Y. Mapping forest disturbance across the China–Laos border using annual Landsat time series. Int. J. Remote Sens. 2019,40, 2895–2915. [CrossRef] 51. Schultz, M.; Clevers, J.G.P.W.; Carter, S.; Verbesselt, J.; Avitabile, V.; Quang, H.V.; Herold, M. Performance of vegetation indices from Landsat time series in deforestation monitoring. Int. J. Appl. Earth Obs. Geoinf. 2016 , 52, 318–327. [CrossRef] 52. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation; Remote Sens. Centre, TEXAS A&M Univ.: College Station, TX, USA, 1973. 53. Xiao, C.; Li, P.; Feng, Z. Monitoring annual dynamics of mature rubber plantations in Xishuangbanna during 1987-2018 using Landsat time series data: A multiple normalization approach. Int. J. Appl. Earth Obs. Geoinf. 2019,77, 30–41. [CrossRef] 54. Purwadi, I.; Van Der Werff, H.M.A.; Lievens, C. Targeting rare earth element bearing mine tailings on Bangka Island, Indonesia, with Sentinel-2 MSI. Int. J. Appl. Earth Obs. Geoinf. 2020,88, 1–8. [CrossRef] 55. Rujoiu-Mare, M.R.; Olariu, B.; Mihai, B.A.; Nistor, C.; Savulescu, I. Land cover classification in Romanian Carpathians and Subcarpathians using multi-date Sentinel-2 remote sensing imagery. Eur. J. Remote Sens. 2017,50, 496–508. [CrossRef] 56. Kauth, R.J.; Thomas, G.S. The Tasselled Cap – A Graphic Description of the Spectral-temporal Development of Agricultural Crops as Seen by LANDSAT. LARS Symp 1976, 4B–41–4B–51. [CrossRef] 57. Crist, E.P.; Cicone, R.C. A Physically-Based Transformation of Thematic Mapper Data-The TM Tasseled Cap. IEEE Trans. Geosci. Remote Sens. 1984,22, 256–263. [CrossRef] 58. Bartsch, A.; Widhalm, B.; Leibman, M.; Ermokhina, K.; Kumpula, T.; Skarin, A.; Wilcox, E.J.; Jones, B.M.; Frost, G.V.; Höfler, A.; et al. Feasibility of tundra vegetation height retrieval from Sentinel-1 and Sentinel-2 data. Remote Sens. Environ. 2020,237, 111515. [CrossRef] 59. Zhang, W.; Brandt, M.; Wang, Q.; Prishchepov, A.V.; Tucker, C.J.; Li, Y.; Lyu, H.; Fensholt, R. From woody cover to woody canopies: How Sentinel-1 and Sentinel-2 data advance the mapping of woody plants in savannas. Remote Sens. Environ. 2019,234, 111465. [CrossRef] 60. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O’Connor, B.; Paganini, M. Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island. Remote Sens. Environ. 2018,215, 517–529. [CrossRef]