Forest Monitoring Based on Sentinel-2 Satellite Imagery, Google Earth Engine Cloud Computing, and Machine Learning

1. Introduction

The increasing occurrence of combined biotic and abiotic disturbances in recent years is putting forests at greater risk of severe damage all around Europe [1]. The frequency and severity of large-scale droughts, fires and insect outbreaks have reached unprecedented levels [2]. In addition, climate change exacerbates both types of disturbances [3], thus studying forest ecosystems on the European level is necessary and could reveal negative responses to forest disturbances.

As is the case in Europe, in Hungary, the occurrence of windfalls, snow breaks, droughts, floods (abiotic), and insect outbreaks, viruses, or fungi (biotic) have all increased in number and intensity in recent decades, which is visible in the damage trends as well [4]. When tree stands are damaged, the yield significantly decreases, and mortality increases. This results in lower photosynthetic activity, and tree vitality [4]. Damage chains can appear as well, since the probability of biotic damage made by insects is higher after severe abiotic damage such as a drought or windbreak, when the defence mechanism of trees is weakened and more vulnerable.

Both biotic and abiotic damage can be detected by remote sensing on which complex forest monitoring systems can be built [5], [6]. By remote sensing (RS) we can monitor the physical characteristics of forests from a distance by utilizing emitted and reflected solar radiation [7]. Satellite-based RS methods offer great opportunities to collect comprehensive data about the current state of health and temporal dynamics of forests, with a high temporal (from 2-5 days up to intraday) and spatial (from 10 m to submeter) resolution [8]. In addition to abiotic and biotic damage, satellite images also detect anthropogenic activities that affect forest health, such as forest management with clearcuts, selective cutting, and reforestation.

The satellites described in this article are equipped with optical sensors that provide data in near-real time, which is ideal for large-area forest mapping and monitoring all disturbance types [9].

The utilization of optical satellite image series data provides a powerful, precise, and cost-efficient means of detecting changes in forest health [10], [11]. In 2008, the National Aeronautics and Space Administration (NASA) made the moderate, 15x15 - 30x30 m spatial and 18-day temporal resolution Landsat 8 imagery available for free [12]. This, together with the European Space Agency’s (ESA) Sentinel-2 (S-2) satellite imagery, available since 2015, with 10x10 m spatial resolution and 2–5-day revisit time, has multiplied the available amount of geospatial data, making it possible to develop forest monitoring approaches and systems, even with harmonized, long-term datasets.

Forest monitoring approaches can be divided into two groups: systems and services. Systems are created for entire countries, utilizing a national database, developed by professionals from institutes with specific forestry knowledge. The results are published on interactive web pages, and the systems’ services are updated regularly. A forest monitoring system was developed in Hungary, called ‘TEMRE’, capable of detecting different types of damage based on vegetation indices derived from MODIS [13] and Sentinel-2 and 3 imageries [11]. The change in forest health is monitored in wall-to-wall systems utilizing harmonized Landsat and Sentinel-2 image series in the Slovak Republic [14], and Sentinel-2 and Planet in the Czech Republic [15]. The Czech system focuses on spruce bark beetle (Ips typhographus) damage which has become a severe issue lately in several European countries [16]. Buras et al. [2] observed and quantified drought impacts with the German Forest Condition Monitor utilizing MODIS images on forests in Germany and Europe, which is also a frequent problem all over the continent.

The second group of monitoring approaches is services. They are usually developed and tested for study areas and specific reasons, e.g., a natural disaster. The results are published in articles and the maps are not available online for access. The scale of these services ranges from nationwide ones to smaller study areas.

Landsats were utilised to investigate how the spatial and temporal characteristics of the Landsat archive can support forest monitoring in Finland [17]. Kern et al. [18] focused on biotic damage in certain parts of Hungary and Croatia caused by the oak lace bug (Corythucha arcuata), while other Sentinel-based development projects have begun in several countries (Poland, Czech Republic, Slovak Republic) for several reasons, including forest health change [19], [20]., treefall gap [21], defoliation [22], and forest succession detection [23].

In both types of monitoring approaches Machine Learning (ML) can be applied, which is an application of artificial intelligence that can be used in the Google Earth Engine (GEE) cloud computing platform [24], to expand the possibilities of finding forest disturbances automatically [25], as well as to classify tree species, which can support the proper interpretation of forest maps. Several ML algorithms such as Random Forest (RF), Minimum Distance Estimation, Support-Vector Machine, k-Nearest Neighbour Regression, and Gradient Boost Regression Tree can be utilized [26] to achieve these goals. This study applied the RF algorithm, which is suitable for vegetation [27] and land cover monitoring [28], but also tree species classification [29]. The above-mentioned ML methods could potentially provide valuable information about forest health. Hence, they could be used to predict future states, which could further enhance monitoring efficiency.

The main objective of this study was to test the novel GEE cloud computing methods on Sentinel-2 imagery to reveal spatial and temporal changes in forest cover and health. The GEE is capable of executing the monitoring tasks in the online cloud: storing, processing, analysing, and visualizing big data enables fast and convenient monitoring of forests [30], [31]. The other objective was to try to see whether validation of satellite maps is possible using ground-based datasets of the Hungarian Forestry Database (FD) and the Hungarian National Forest Damage Registration System where data is collected and stored regularly for the whole country in a uniform way. The ultimate goal is to automatize the monitoring process using both the remotely sensed and field-surveyed datasets in a wall-to-wall monitoring system.

2. Materials and Methods

2.1. Study Site

The Nagyerdő (Great Forest) of Debrecen is situated north of the city of Debrecen (21.63° N, 47.57° E) in eastern Hungary (Fig. 1) and covers a 1092 ha, contiguous, protected area [32]. In the past, the loess and sandy soils in the area had a favourable water supply that maintained large, continuous forests. Today, a drier climate, river regulations, and forest cutting have reduced the forest to smaller patches. The Nagyerdő is a good example of the effects of these unfortunate anthropogenic changes.

Figure 1. Location of Nagyerdő of Debrecen in Eastern Hungary.

Figure 1. Location of Nagyerdő of Debrecen in Eastern Hungary.

The most typical forest communities in the Nagyerdő are oak with lily of the valley (Convallario - Quercetum roboris), but oak-hornbeam (Querco robori - Carpinetum) and open oak forest on sand (Festuco rupicolae - Quercetum roboris), which were common in previous times, have disappeared.

The forest tree composition includes several species, but the most important, largest (20-30 m high), and oldest (100-300 years) pedunculate oaks (Quercus robur) can be found in the old forest. Domestic compound species are silver poplar (Populus alba), wild cherry (Prunus avium), Tatar maple (Acer tataricum), field maple (Acer campestre), field elm (Ulmus minor), wych elm (Ulmus glabra), European crab apple (Malus sylvestris), European wild pear (Pyrus pyraster), large-leaved lime (Tilia platyphyllos) and silver lime (Tilia tormentosa). The foreign species are black locust (Robinia pseudoacacia), Scots pine (Pinus sylvestris), and northern red oak (Quercus rubra).

The forest has been suffering for decades from increasing abiotic damage as well as from more and more frequent drought stress during the spring and summer due to decreased precipitation, and a lower groundwater level, resulting in decreasing yield and lower vitality [33], while ice damage takes place during springs when seedlings and young trees freeze, and game damage is also observed in sites of reforestation.

The forest is divided into four quarters; drought and frost damage occur mostly in the NW and NE quarters, but it has also occurred in the SW and SE quarters, to a lesser degree. One hundred thirty-one forest compartments out of 394 (33%) were damaged.

2.2. Satellite-Based Dataset

Sentinel-2 imagery from the European Space Agency has been available since 2015 with a 2-5 days revisit time (with twin satellites), high 10x10-60x60 m spatial resolution, and a free and open data policy, which makes it ideal for forest monitoring (ESA 2022). On board, Sentinel-2 twin satellites, Multispectral Instruments (MSI) are installed to measure the radiance of the Earth’s surface in 13 spectral bands (Tab. 1).

Table 1. Spectral bands of Sentinel-2 MSI. .

Band Number	Bands	Central Wavelength (µm)	Resolution (m)
Band 1	Coastal aerosol	0.443	60
Band 2	Blue	0.490	10
Band 3	Green	0.560	10
Band 4	Red	0.665	10
Band 5	Vegetation Red Edge	0.705	20
Band 6	Vegetation Red Edge	0.740	20
Band 7	Vegetation Red Edge	0.783	20
Band 8	Near-infrared	0.842	10
Band 8A	Vegetation Red Edge	0.865	20
Band 9	Water vapour	0.945	60
Band 10	Short-wave infrared Cirrus	1.375	60
Band 11	Short-wave infrared	1.610	20
Band 12	Short-wave infrared	2.190	20

Table 1. Spectral bands of Sentinel-2 MSI. .

Band Number	Bands	Central Wavelength (µm)	Resolution (m)
Band 1	Coastal aerosol	0.443	60
Band 2	Blue	0.490	10
Band 3	Green	0.560	10
Band 4	Red	0.665	10
Band 5	Vegetation Red Edge	0.705	20
Band 6	Vegetation Red Edge	0.740	20
Band 7	Vegetation Red Edge	0.783	20
Band 8	Near-infrared	0.842	10
Band 8A	Vegetation Red Edge	0.865	20
Band 9	Water vapour	0.945	60
Band 10	Short-wave infrared Cirrus	1.375	60
Band 11	Short-wave infrared	1.610	20
Band 12	Short-wave infrared	2.190	20

Sentinel-2 satellite data processing, analysis, and visualization took place online in the cloud using the Application Interface of GEE [34]. GEE uses both JavaScript and Python programming languages, and, in this study, JavaScript was chosen. The method used in GEE consists of several steps, shown in a flowchart (Fig. 2).

Figure 2. Methodological flowchart of processing, visualizing, and exporting Sentinel-2 images in GEE for forest damage monitoring.

Figure 2. Methodological flowchart of processing, visualizing, and exporting Sentinel-2 images in GEE for forest damage monitoring.

The ESA Sentinel-2 MSI L2A satellite image mosaics [8] provided the basis for the remote-sensing survey at 10x10 m spatial resolution as the first step in the process. The bottom of atmosphere reflectance-based (level 2, Sentinel-2_SR collection in GEE) annual composites were created with GEE for each year between 2017 and 2020 using spatial, temporal, and quality filtering: region of interest (ROI), time window, and cloud coverage. The annual composites are made from the S-2_SR collection by applying the pixel-wise filter and reducer functions of GEE, which are essential to creating cloud-free images for vegetation periods of each year in the study area.

Spatial filtering by filterBounds was calibrated to filter pixels outside of the study area (the same as ROI), which was given in a shapefile polygon and a WGS84 (EPSG 4326) coordinate reference system. The polygons are from the Hungarian Forest Database, and they were filtered to include only forest-covered areas and exclude forest compartments with water surfaces, roads, buildings, openings, nursery gardens, game areas, and shrubbery. This step tends to increase the monitoring accuracy.

Similarly, we defined a time window by filterDate based on the vegetation period of each year (e.g., filterDate('2017-04-15', '2017-10-15')) to keep pixels from the active state of the vegetation. Time filtering is also important for reducing the dataset since GEE is not capable of processing high-resolution data for larger areas and longer periods.

Cloud filtering was completed in two steps. The first step involved pre-filtering cloudy pixels (ee. Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 5)) that were less than or equal to 5% of coverage, which was applied on the S2 collection. The maximum 5% of cloud coverage limit was necessary to keep enough images for mosaicking. After this, the second step utilized the QA60 quality band as a bitmask, where 0 values of 10 and 11 bits refer to cloudless pixels. The two-step filtering resulted in a nearly cloud-free image collection for the selected study areas and period.

The second step involved the calculation of vegetation and water indices. Normalized Difference Vegetation Index (NDVI) (1), Normalized Difference Water Index (NDWI) (2), and Enhanced Vegetation Index (EVI) (3) maps were created from the filtered composites. Change maps were calculated from state maps based on differences between years. NDVI change (4) is based on the subtraction of two following years while standardised Z NDVI (5) is calculated by NDVI standardization.

The Normalized Difference Vegetation Index (NDVI) (1) is the most widely used index for indirectly measuring the photosynthetic activity of vegetation. It is calculated according to the following formula [35]:

NDVI = (NIR - RED) / (NIR + RED)

(1)

where RED is band 4 and NIR (Near-Infrared) is band 8 in the case of Sentinel-2. High NDVI values refer to strong photosynthetic activity, high vegetation cover, and healthy vegetation, while low index values refer to unhealthy vegetation or lower vegetation cover. The range of NDVI is between -1 and +1, where 1 is theoretically healthy vegetation, 0 is a vegetation-free area, and below zero is the water surface. In practice, a healthy forest reaches 0.8-0.9 depending on geographical location, tree species, and vegetation period. The ee.image.normalizedDifference (['B8', 'B4']) function was used to calculate NDVI in GEE.

The Normalized Difference Water Index (NDWI) (2) is used for detecting the water content of the forest canopy and can be calculated by this formula [36]:

NDWI = (NIR-SWIR) / (NIR+SWIR)

(2)

where NIR (Near-Infrared) is band 8 and SWIR (Short-Wave Infrared) is band 11 in the case of S-2. High NDWI values correspond to high vegetation water content and cover, while low index values correspond to low vegetation water content and cover. NDWI decreases during a period of drought, water stress, or after logging, which makes it ideal for canopy water content monitoring [37], [38]. The range of NDWI is similar to NDVI, where -1 refers to the driest state and +1 to the wettest state.

Compared to NDVI, the Enhanced Vegetation Index (EVI) (3) has several correlation factors in its formula to correct known issues [39]:

EVI = 2.5 * ((NIR - RED)) / ((NIR + C1 * RED - C2 * BLUE + L))

(3)

where NIR (band 8), RED (band 4), and BLUE (band 2) have atmospherically corrected surface reflectance bands, L is the canopy background adjustment for differential NIR and red radiant through a canopy, and C1 and C2 are the resistance coefficients, which utilize the blue band to correct aerosol influences in the red band. The adapted coefficients are L = 1, C1 = 6, C2 = 7.5, and G (gain factor) = 2.5. Healthy vegetation has an index value between 0.2 and 0.8 depending on the forest type. It is important to note that although EVI is usable for Sentinel-2, it was created for MODIS and MODIS-based coefficients.

The change indices (4) are calculated from all three indices (NDVI, NDWI, EVI) above by subtracting a given year’s value (second) from the previous year’s value (first) according to this formula:

NVDI_change = NDVI_{second year} - NDVI_{first year}

(4)

NDWI and EVI changes were calculated similarly; the range is, in theory, from -1 to +1. However, these change indices only show the difference in forest state between two exact dates and might not be sufficient to detect forest damage due to interannual variation of index values.

To reduce this variation, standardized NDVI values (Z NDVI) (5) [40] were calculated by this formula:

Z NDVI = (NDVI - NDVI_mean) / NDVI_std

(3)

where NDVI is the NDVI of the actual year, NDVI_mean the multiple-year average of NDVI (2017-2020 in this case), and 𝑁𝐷𝑉𝐼_std is the standard deviation of NDVI values (also the 2017-2020 period). The range of Z NDVI could have large variations between the end values, most often between -3 and 3 in the case of the Nagyerdő, where the range of -3- -1 refers to serious disturbance, -1 - +1 to a slightly changed health state in a negative or positive direction, respectively, and + 1-3 stands for stronger regeneration.

The third step was the aggregation of the dataset of VI with reducers, which is needed to visualize and export the dataset due to the limitations of GEE. Different filters, such as median (.median()), mean (.mean()), and standard deviation (STD) (.stdDev()), were used for reduction. Mean and STD are used in the formula of Z NDVI, while the median is for displaying and exporting maps. The creation of true colour composites (RGB) was the fourth step of the process, using red (band 4), green (band 3), and blue bands (band 2) and a median filter as well. Without these reducers, the GEE is not capable of calculating mosaics due to the limitations and the applicability of certain functions.

Setting the displaying parameters and properties was the next step: the map center (with ROI polygon as center and zoom level), colour palettes (green-red transitional scale for NDVIs and EVI, and white-blue for NDWI with 5-5 classes), and borders (colour, width) were described for visualization. These parameters were used in the sixth step when the actual visualization occurred on the GEE platform: RGB maps, VI maps, and charts as well. The charts were made with the ui.Chart.image.series function and based on the complete five-year dataset, where every index map was reduced with a median reducer.

The seventh and final step involved the export of the maps with Export.image.toDrive and charts with print(ui.Chart.image.series) to Google Drive and PC in georeferenced TIFF, CSV, and PNG formats, respectively. This offers the possibility of further offline analysis and more detailed image visualization. For the latter, QGIS 3.22.3. (“Białowieża”) open-source desktop GIS software was used [41].

The tree species classification was made with ee.Classifier.smileRandomForest classifier. A median composite of 2021 was used as the base of the dominant species classification: pedunculate oak, black locust, Scots pine, and red oak, plus clearcuts. Every S-2 band was used in the classification and combined into a single median raster as an input.

The sampling approach was based on the manually selected training points on the RGB image and these points were used as training input for tree species classification. Five classes were defined by the 123 points (with the geometry imports function) which were based on the tree species dataset of the Hungarian Forestry Database, and the RF classification was completed with the training samples (15 black locust, 22 red oak, 27 Scots pine, 40 pedunculate oak, 19 clearcut) and 100 decision trees. The different tree species layers were also merged into one. The classified map was exported with image.toDrive function to Google Drive and PC.

2.3. Ground-Based Dataset

Ground-based forest damage reports were applied for validation. According to the forest protection damage reports of the Hungarian National Forest Damage Registration System (HNFDRS) of the Hungarian National Land Centre, late frost damage occurred in young stands in 2017, while in 2018, 2019, and 2020 [42,43] the Nagyerdő suffered from drought, frost, and game damage. The damage type was registered in the reports and is available on the HNFDRS website. The data contains information about damage frequency and intensity (or severity), which is given for each forest compartment expressed as the percentage of damaged trees for each tree species (0-100%) (Fig. 3), and the intensity of damage (0-100%), respectively [42]. The damaged area is given in hectares. Damage ratio (6) was introduced according to the following formula, measured on a 0-100% scale:

Damage ratio = (damaged area)/(total area)*100

(6)

Figure 3. Aggregated damage frequency map of Nagyerdő of Debrecen based on field surveys from 2017 to 2020. The northern half of the forest shows the most considerable severe damage, but overall, at least 33% of the area of Nagyerdő was damaged in one year.

Figure 3. Aggregated damage frequency map of Nagyerdő of Debrecen based on field surveys from 2017 to 2020. The northern half of the forest shows the most considerable severe damage, but overall, at least 33% of the area of Nagyerdő was damaged in one year.

Damage registered in the forest damage reports was compared to Z NDVI values on the pixel level in QGIS. For this, the polygons of reports were rasterized to a 10x10 m resolution grid that is identical to the Sentinel-2 images in size. The field reports were reclassified according to damage ratio: when this was below 30% it was assigned a value of 0, while above 30% the value was 1. In the case of Z NDVI values, -0.5 was selected as a threshold; every pixel was considered forest damage if it fell below this, being assigned a value of 0, while all pixels exceeding the threshold were assigned a value of 1. The result maps were analysed with the accuracy assessment postprocessing function of the Semi-Automatic Classification Plugin (SCP) of QGIS [44].

To validate RF-based tree species classification results, the majority value was calculated for all forest compartments from RF pixels in QGIS (raster analysis, zonal statistics) where the most frequent pixel values (dominant species) within the compartments were compared to the Forestry Database attributes. The attribute of the dominant species is given in the FD and calculated from RF-based pixels as well. The forest compartment-level comparison was based on the number of correctly classified RF and ground-based compartments. The accuracy of tree species classification is calculated in MS Excel according to this formula (7):

Tree species accuracy = (RF compartments)/(FD compartments )*100

(7)

2.4. Statistical Analysis

We performed the Shapiro-Wilks test for normality on 100 randomly selected points from all four Z NDVI rasters. This method was described in the original paper of the Z NDVI formula [40] and the authors there achieved 80% normality. The points were created with random points in the layer bounds function of QGIS and analysed in Past 4.10 statistics software [45]. The distribution was normal (p<5%) in the case of at least 95% of the points in the years 2017, 2018 and 2019 as well, with p values of 0.37, 0.3, and 0.11, respectively. In 2020 this was not the case (0.01) but this was not due to the outliers; the distribution of bins is rather even. The W values were between 0.97 and 0.99.

3. Results

State and change index maps were derived from S-2 composites for each year: 2017, 2018, 2019, and 2020. The state index maps provided information on the actual condition of the forests in the given year, while the change index showed the difference between the sequential years (2017-2018, 2018-2019, 2019-2020). Standardized NDVI values were also calculated utilizing the entire time series, showing deviation from mean of the four years. These Z NDVI maps seemed to be the most promising for further analysis because they exhibited greater inter-annual and spatial changes.

Analysis of Z NDVI maps (Fig. 4) showed that the composite of the base year 2017 contained large healthy, dark green coloured areas with high Z NDVI values but also damaged, orange or red areas with clearcuts and other sparsely vegetated forest compartments. In 2018, a general Z NDVI decrement was detected (orange colour) in most of the forest with new red areas representing clearcuts following the borders of forest compartments. A positive anomaly caused by artificial plantings after clearcuts, and the regeneration of less-forested areas appeared on the maps in 2019 and 2020. Despite the general positive changes, negative changes were also apparent. New clearcuts occurred every year. Moreover, the oak-dominated quarter II of the forest became orange coloured on the map in 2020. This negative change is due to drought stress caused by the lack of precipitation to which pedunculate oak is vulnerable [46].

Figure 4. Annual Z NDVI composites of Nagyerdő of Debrecen in 2017 (a); 2018 (b); 2019 (c); and 2020 (d). Clearcuts (with dark red), drought, and frost damage (orange or yellow) are visible on every map marking a significant photosynthetic activity drop in the given years compared to the long-term mean. The difference between VI values is also due to different dominant tree species. Regeneration was detected as well with green colours.

Figure 4. Annual Z NDVI composites of Nagyerdő of Debrecen in 2017 (a); 2018 (b); 2019 (c); and 2020 (d). Clearcuts (with dark red), drought, and frost damage (orange or yellow) are visible on every map marking a significant photosynthetic activity drop in the given years compared to the long-term mean. The difference between VI values is also due to different dominant tree species. Regeneration was detected as well with green colours.

Natural differences have been observed between the studied years. For example, both the ground-based and the RS-based methods detected that several forest compartments were damaged in 2017. However, these differences are less evident in 2018. The year 2018 was generally drier than 2017, and the whole forest appeared to be in worse condition according to the Z NDVI values, while 2019 was a better year with more positive values and very few damage reports. The year 2020 was worse again, with more damage reports, lower VI values, and a visible phenomenon in quarter II, which is pedunculate oak-dominated and suffers from a decades-long history of heat and drought stress. However, there is no perfect match between the ground-based survey of disturbances and RS data, which could be due to weaker, no-disturbance changes in forests.

Confusion matrices were created for each year to compare classified Z NDVI values to reference ground-based damage reports (Tab. 2) at the pixel level. Thresholds were selected to focus on more severe damage; thus, the damage ratio was > 30% and Z NDVI was < – 1. According to these rules, every pixel was reclassified and was given a value of 0 in the case of damage and 1 in the case of no damage. User’s (UA), Producers’ (PA), and Total Accuracy (TA) were calculated for all years (Congalton 2021). The mean TA for all four years was 70.85% (Tab. 3). PA was high (99-100%) but UA varied between 62.7 and 74% and TA between 63.24 and 74.51%.

Table 2. Confusion matrix of forest damage datasets for years 2017-2020. Values are given in ratio, in % format.

		2017		2018
		Reference
		Damaged	Non-damaged	Damaged	Non-damaged
Classified	Damaged	72	25	55	32
	Non-damaged	3	0	1	12
		2019		2020
		Reference
		Damaged	Non-damaged	Damaged	Non-damaged
Classified	Damaged	75	25	72	28
	Non-damaged	0	0	0	0

Table 2. Confusion matrix of forest damage datasets for years 2017-2020. Values are given in ratio, in % format.

		2017		2018
		Reference
		Damaged	Non-damaged	Damaged	Non-damaged
Classified	Damaged	72	25	55	32
	Non-damaged	3	0	1	12
		2019		2020
		Reference
		Damaged	Non-damaged	Damaged	Non-damaged
Classified	Damaged	75	25	72	28
	Non-damaged	0	0	0	0

Table 3. Accuracy assessment of forest damage datasets for years 2017-2020. .

	2017	2018	2019
Producer's Accuracy (%)	99.19	99.89	100
User's Accuracy (%)	74.01	62.69	74.62
Total Accuracy (%)	73.70	63.24	74.51

Table 3. Accuracy assessment of forest damage datasets for years 2017-2020. .

	2017	2018	2019
Producer's Accuracy (%)	99.19	99.89	100
User's Accuracy (%)	74.01	62.69	74.62
Total Accuracy (%)	73.70	63.24	74.51

Differences between VI values within forest compartments are not exclusively due to damage but are also due to different tree species, which were detected both on the RF-based classified and the VI maps.

RF classification showed the presence of coherent patches of black locust, red oak, and Scots pine in quarters I and III (often corresponding to forest compartment borders), while quarters II and IV were dominated homogeneously by pedunculate oak (Fig. 5), with smaller patches of other species.

Figure 5. Main tree species of Nagyerdő made with Random Forest classification of Sentinel-2 images from 2020.

Figure 5. Main tree species of Nagyerdő made with Random Forest classification of Sentinel-2 images from 2020.

The majority-based species comparison on the forest compartment level indicated 76.1% accuracy in the case of five classes (four dominant tree species plus clearcut). Of the 339 compartments, in 258 we found agreement on the species. However, the method applied here aimed to classify the dominant tree species; there are several compartments with secondary, tertiary, quaternary, or even more mixed species, which made the result less accurate. It is also important to note that there is ongoing tree species change in several compartments due to nature conservation when black locust and pine are exchanged with pedunculate oak. This process resulted in compartment-sized clearcuts and lower classification accuracy.

The comparison of VI and classified tree species maps indicated a connection between them, and different species reacted differently to forest damage, which resulted in different VI values. The varying VI values of different species were visible on the Z NDVI maps of 2017 and 2018 (Fig. 4a, b), in the bottom part of quarter I, where black locust, red oak, Scots pine, and pedunculate oak-dominated stands were situated in a row. The stands with black locust and the pedunculate oak majority had lower vitality and VI values, while pine and red oak were in better condition compared to them. Scots pine and pedunculate oak showed higher Z NDVI values in 2017 (with dark green) compared to the other species mentioned above; however, in the following years, vitality oak stands decreased, and larger orange areas appeared in certain forest compartments in 2019 and 2020 due to drought stress, especially in quarter II (Fig. 4 c, d). Frost damage in young stands can be seen as well in 2017, 2018, and 2020 in several compartments in quarters I and III.

Besides the maps, time series can also be analysed in charts; thus these were made from the Sentinel-2 imagery showing median NDVI values of the Nagyerdő between 2017 and 2020 (Fig. 6). Deviations from the ideal state of the vegetation period can be detected when the actual values are compared to the ideal NDVI curve. If forest damage takes place, the shape of the NDVI curve will change according to the severity of the damage. Drought affected the NDVI curve in 2018, 2019, and 2020, which is visible on the graph in the form of a decline.

Figure 6. NDVI median graph of Nagyerdő between 2017 and 2020 in the GEE monitoring system. Deviations from the ideal curve are marked with a black frame, which could refer to forest damage.

Figure 6. NDVI median graph of Nagyerdő between 2017 and 2020 in the GEE monitoring system. Deviations from the ideal curve are marked with a black frame, which could refer to forest damage.

4. Discussion

The application of Sentinel-2 imagery in GEE has several advantages that made the forest monitoring successful; however, it had some limitations as well. Utilizing simultaneously the remotely sensed and field-surveyed datasets in a wall-to-wall monitoring system is possible, also in an automatized way but improvements are needed to fix the problems before running an operational system.

One problem was high cloud cover that caused gaps in the dataset, especially in April, May, and June, months during which cloudless images of the area of interest were virtually unavailable for some years. April, May, and June are particularly important months in the vegetation period because they mark the start of the growing season (leaf unfolding) with increasing photosynthetic activity. According to Swedish studies [47,48], these months provide valuable information concerning the spring phenological stage and the RF classification is more accurate when these months are used. Still, cloud masking remains an existing problem in the S-2 dataset; however, an increasing number of refined datasets and algorithms such as the Cloud Probability dataset, the s2cloudless, or the Fmask algorithm were written in R, Python, or JavaScript languages, and are available in GEE [49].

The other problem was the validation of the RS results with the ground-based dataset results. The disadvantage of the ground-based dataset is the large sampling area (several hectares) compared to satellite pixels (10x10 m), which results in an uneven distribution of data. This different sampling method is visible when scattered pixels are compared to whole forest compartments. To solve this problem, mean or median values could be calculated on a zonal basis from the pixels for whole compartments, but this would significantly reduce the diversity of VI values and make it harder to detect any changes in the health state of the forest.

The advantage of the compartment-based comparison method is its easier application for practical use since the VI maps made for forest compartments of the Forestry Database are used by forest managers. It could reduce the workload of fieldwork, thus using damage reports made by foresters could represent an alternative to repetitive fieldwork triggered by each new research study.

Another problem was the date of the ground-based registration. If data is collected after satellite image acquisition, the overlap will not be visible on the maps dedicated to showing the same year. The RS could show the damage which had occurred before the ground survey, which could be a useful tool to survey the damaged area before going out into the field. Drought effects can appear in the coming year rather than the given year; thus, it is not always completely evident which datasets will be compared to each other. The difference in sampling size is visible between the Producer’s and User’s Accuracy as well, where ~99% and 70% showed a remarkable difference in the accuracy of damage detection.

Studies from other regions of the world which deal with ML-based forest monitoring showed varied results. A temperate forest degradation study was made by Chen et al. [50] for Georgia using GEE and Landsat imagery, where the UA of the forest degradation class was 69% and the PA was 83%, while tropical forest disturbances were studied by Chen et al. [51] using harmonized Landsat and Sentinel-2 data, showing 84.5 and 95.5% TA for Tanzania and Brazil respectively. Later, in the boreal region, Yang [52] investigated biotic forest damage caused by the spruce bark beetle in Sweden on multitemporal Sentinel-2 collection, with maximum likelihood and EVI, GNDVI (Green Normalized Difference Vegetation Index) indices, and 89% TA was achieved.

Machine learning was successfully applied to Sentinel-2 data aimed at tree species classification to support VI maps. The combination of cloud systems and ML enabled the visualization of the distribution of different forest types and different forest damage types in all the investigated years. However, some previous studies reported better results using multi-temporal S-2 images with pedunculate oak, Scots pine, silver birch (Betula pendula), dunkeld larch (Larix × marschlinsii), and Norway spruce in Sweden with RF resulting in 88.2% Total Accuracy [47]. Fourteen S-2 images were used in the multi-temporal imagery resulting in 88.2% TA, but for a single image in July, it was 75.9%, which is similar to our result. Also in Sweden, 87% TA was obtained for the same species [48] with the Bayesian inference method using 23 pieces of S-2 images. From Latvia, 92-94% TA was reported with three S-2 images and K-means clustering and DynLand methods on Scots pine, Norway spruce, silver birch, and black alder (Alnus glutinosa) [53]. In Italy, Puletti et al. [29] obtained 86.2% TA with the RF method on four mixed forest types, emphasizing that the collection of multitemporal images at different phenological periods is required. Based on these studies, it is evident that multi-temporal imagery provided better results than single images.

It would be beneficial to test other ML methods on our study site as well, and several supervised and unsupervised algorithms are available in GEE besides RF. These algorithms can replace the process of downloading and processing satellite images as well as the need for external programs since the whole process can be made in the cloud.

On the other hand, ground-based datasets like forest damage reports cannot be uploaded directly to public GEE servers due to security issues and policies. This could make it impossible to carry out a fully online analysis, which would be desirable, especially on a larger scale (i.e. for the whole of Hungary or even Europe). Thus, the method has a potential for large-scale usage; however, its performance in different climatic regions needs to be validated first and data security has to be taken into consideration. The ML can be used to find forest damage as well, besides the tree species classification, which could be an expansion of this study in the future.

Utilized datasets could be expanded too. Since drought has become an increasing problem both worldwide and in Hungary, climatic and meteorological datasets (precipitation, temperature, soil moisture, etc.) could be compared to VI values to show correspondence(s) between them, as Birinyi et al. [54] suggested with corn yield comparison to NDVI and EVI in Hungary with GEE.

Ultimately a wall-to-wall system will be created for the entire of Hungary based on cloud-free Sentinel-2 composites, cloud computing powered by GEE and supported by ML. The application of ML could ensure both more accurate damage threshold determination and damage detection as well. A fully automatized, regularly updated country-wide map would be a great asset for both researchers and practical foresters in forest surveying and damage monitoring.

5. Conclusions

Results indicated that the combined dataset of VIs derived from satellite imagery and ground-based reports provided suitable input for forest damage monitoring conducted with GEE. The applied method successfully identified different types of forest damage (drought and frost) on Z NDVI maps in the surveyed period with 70.85% Total Accuracy, while tree species classification (black locust, red oak, Scots pine, pedunculate oak) with RF showed 76.1 % TA. These values are satisfactory but could be higher with the dynamic damage threshold computed by algorithms and removing the clearcut class and only focusing on forest-covered compartments. This latter resulted in 81.79 % accuracy instead of 76.1 %.

This GEE and RF-based method is of great importance when forest damage is more and more frequent all-around Europe due to climate change, and the fact that several tree species are unable to adapt to warmer climatic conditions, thus constant and accurate monitoring is needed to obtain data about forest health change.

The great advantage of GEE is resilience and flexibility: study sites can be modified easily, and the creation of maps and charts can be achieved rapidly for other areas covered by the same RS dataset, and Sentinels and Landsats cover the mainland surfaces of Earth in 5-8 days.

Cloud filtering is also an easily modifiable part of the GEE code which is an essential part of forest monitoring. The cloudy pixel percentage filtering and the QA 60 bitmasking used in this study proved satisfactory results since high-quality mosaics were made for each studied period but there were gaps in the dataset during cloudy springs or summers. However, novel cloud-filtering algorithms like Fmask or s2cloudless can be applied to S-2 images and easily integrated into the codes, and with adequate cloudless images created with novel cloud filtering methods, it will be possible to analyse forest health on a weekly or monthly basis instead of on annual composites. On these denser series, seasonality in vegetation indices can also be studied more accurately. Although there could be still areas (high-mountains, islands) with frequent cloud cover where it is not possible to achieve cloudless images for some periods, on mainland areas the fusion of Sentinel and Landsat could solve the problem of cloud coverage [55].

The RS-based approach supported by GEE could reduce the workload of the forestry labour force as well and provide forest data more quickly and more accurately than the traditional methods used in fieldwork. The usage of field reports and forestry data from national forest damage databases could enhance the possibility of expanding the method into other regions of Hungary, since the datasets are available for the entire country in a uniform way, and no extra field measurements would be needed. That is one of the reasons why we used this database in this study instead of conducting our measurements in the field. As we hypothesized, the field reports can be used for validation despite the different scales; however, the accuracy can be improved. The result of our study can also help to carry out field measurements and reports more effectively. To enhance accuracy the thresholds of damage could be determined by ML instead of a manual setting. ML could detect trends in anomalies of time series which sit outside these patterns and identify the thresholds.

The suitability of the proposed approach would be desirable to test on a larger scale. The expansion of the novel GEE-based monitoring method to a wall-to-wall monitoring system would be possible, for Hungary or even Europe. The similar conditions and forest damage types in Europe could be monitored uniformly regarding data collection, processing, analysis, and visualization, which would be a novel solution instead of national approaches. Cloud computing and machine learning both have great relevance in this plan; with them, the entire monitoring process can be made online.

Besides spatial expansion, other types of forest damage could be surveyed with the GEE-based method, namely the spruce bark beetle and oak lace bug gradations, which have both caused massive damage in several European countries. This monitoring on the European level could show the outbreak and spread of the species from one country to another. The oak lace bug was studied in Hungary and Croatia in these terms but on MODIS imagery, not on Sentinel-2. The bark beetle-related studies also focused on study areas or single countries; thus, it would be interesting to try the method for a larger-scale damage survey.