Exploring the potential to use in-between pixel variability for early detection of bark beetle attacked trees

. The European spruce bark beetle ( Ips typographus L.) is a major disturbance agent in Norway spruce ( Picea abies (L.) Karst ) forests in Europe and it is estimated that a changing climate will result in more severe outbreaks in the future. To reduce the risk of large outbreaks it is important to have methods that enable early detection of bark beetle attacks to help forest managers to prevent population build-up, e.g by sanitary cutting. Several studies have been devoted to early detection of bark beetle attacks with Sentinel-2 data with a focus on spectral properties and vegetation indices for early detection with pixel-based methods. In this study we explore the potential to use changes in variability between pixels in windows of different sizes (3×3, 4×4 and 5×5 pixels). We compute the coefficient of variation for four vegetation indices (NDVI, NDWI, CCI and NDRS) in a time-series of Sentinel-2 data during a bark beetle outbreak in Sweden that was triggered by a drought in 2018. The results indicate that CCI is the most promising index for early detection and that the variability between pixels increase in windows with attacked trees from late July when the main swarming was the second week of May.


Introduction
The European spruce bark beetle (Ips typographus L.) is a major disturbance agent in Norway spruce (Picea abies (L.) Karst) forests in Europe (Jönsson et al. 2012) and has killed more than 150 million m 3 forest during the last 50 years (Schroeder and Cocoş 2018). Bark beetles are naturally present in spruce forest and at lower population levels they only attack weakened trees that are comparatively easy to kill; but periodically, favourable conditions enable the populations to grow rapidly, enabling the bark beetles to also attack healthy trees, resulting in large outbreaks (Hlásny et al. 2021). Such outbreaks often occur after windstorms since storm-felled trees provide the bark beetles with breeding material that facilities population growth (Schelhaas et al. 2003). Bark beetle outbreaks can also be triggered by drought events as drought stress decreases the trees' defense capacity (Jactel et al. 2012;Netherer et al. 2019). As an example, the drought that occurred during the summer of 2018 triggered the largest recorded bark beetle outbreak in Sweden.
With a warmer climate, more severe bark beetle attacks are projected, mainly due to higher temperatures and more frequent drought events (Jactel et al. 2012;Marini et al. 2013); a relationship between increasing temperatures and forest mortality due to bark beetles has already been observed (Seidl et al. 2014;Hlásny et al. 2021). In addition to abiotic factors that influence the trees' defense capacity, the bark beetles themselves are influenced by temperature. Northern regions commonly have one generation of bark beetles per year, but during longer and warmer summers the probability of bark beetles developing two generations further north increases, resulting in more rapid population growth (Jönsson et al., 2012;Bentz et al. 2019).
To prevent the beetle population build-up, e.g. by sanitary felling, it is important to detect infested spruce trees early, while the larvae are still inside the attacked trees (Jönsson et al. 2012;Hlásny et al. 2021). Several studies have been conducted to detect attacked trees in this so-called green attack stage with satellite based remote sensing. Abdullah et al. (2018) showed that using Sentinel-2 data resulted in higher detection accuracy (67%) than Landsat (36%) for the green attack stage. Fernandez-Carrillo et al. (2020) applied Sentinel-2 data to classify bark beetle damage into AGILE: GIScience Series, 4, 35, 2023. https://doi.org/10.5194/agile-giss-4-35-2023 Proceedings of the 26th AGILE Conference on Geographic Information Science, 2023. Editors: P. van Oosterom, H. Ploeger, A. Mansourian, S. Scheider, R. Lemmens, and B. van Loenen. This contribution underwent peer review based on a full paper submission. © Author(s) 2023. This work is distributed under the Creative Commons Attribution 4.0 License. severity classes and achieved high accuracy for areas with high damage but lower accuracy for areas with less severe damage. Huo et al. (2021) developed a new vegetation index, Normalized Distance Red & SWIR (NDRS), that performed better than other Sentinel-2 derived indices for early detection of stressed spruce trees. Dalponte et al. (2022) achieved an accuracy of 83% when applying a support vector machine classifier to Sentinel-2 data to detect early and late stages of bark beetle attacks.
However, these studies have focused on the spectral properties of individual pixels and not on the spatiotemporal development of bark beetle attacks over adjacent pixels. In this study we compute time-series of the coefficient of variation for windows of different sizes (3×3, 4×4 and 5×5 pixels) to explore if changes in the variability between pixels can be used for early detection of bark beetle attacks. The assumption is that the variability between pixels increases in a window with bark beetle attacked trees since the attacked trees at lower populations often are scattered or small groups of trees, and that the increase in variability might be detectable even if the change in a vegetation index is not sufficiently large for detection on pixel level.

Methods and study area 2.1 Study area
The study area was an approximately 8×8 km large area located in Västervik municipality, Sweden (Fig. 1). The area is dominated by coniferous forest but contains some deciduous forest, agricultural fields and lakes.

Sentinel-2 data and processing
Sentinel-2 L1C data was downloaded and processed to an analysis-ready data cube with the Framework for Operational Radiometric Correction for Environmental monitoring (FORCE; Frantz 2019) version 3.6.5. Wavelength bands with 20×20 m spatial resolution were resampled to 10×10 m by splitting a pixel into four without modifying the pixel values. Image-to-image registration within FORCE was utilized (Rufin et al. 2021) to handle the geometric deviations between Sentinel-2a and Sentinel-2b. where DRSmax and DRSmin are maximum and minimum DRS for coniferous forest in a Sentinel-2 image respectively (Huo et al. 2021).

Bark beetle data
Harvester data was obtained from Sveaskog, the largest forest owner in Sweden, and consisted of coordinates of attacked trees derived from harvester machines when cutting the trees. The harvester data was rasterized to create a grid with the number of attacked trees per pixel that aligned with the FORCE grid (10×10 m).
The bark beetle damage data only included attacked trees. Data over healthy pixels were derived with the aid of land cover data and forest estate data. The assumption was that in a forest estate with bark beetle attacked trees the forest manager would have marked all attacked trees for the harvester to cut since it is associated with high costs to bring the harvester. Hence, coniferous forests in the same estate as forests with harvested trees, but without harvester data, were considered healthy. The forest estates are relatively small and the trees were harvested in late stages of the attack so nearby attacked trees would have been easily detected. For more details about the damage data see Müller et al. (2022).
The date of attack was estimated from bark beetle swarming data published by the Swedish forest agency (Skogsstyrelsen, 2023). The swarming data gives the number of beetles caught in pheromone traps on a weekly basis. For 2018 the station with pheromone traps near the study area showed a distinct peak in the number of caught beetles in week 19 (second week of May), followed by weeks with low number of beetles. Hence, the trees were likely attacked during week 19 in 2018.

Calculate coefficient of variation
We identified windows of three sizes: 3×3, 4×4 and 5×5 pixels. The windows were identified by iterating over the dataset with a moving window. This means that some windows overlap, and that a single pixel could theoretically be part of up to nine unique windows among the 3×3 pixel windows, up to 16 windows in the 4×4 pixel windows, and up to 25 in the 5x5 pixel windows.
Two types of windows were identified for each window size: attacked windows, and healthy windows. Attacked windows are windows for which at least one pixel includes trees with confirmed bark beetle attack and where all pixels in the window are of the same land cover type (Fig. 2). The reason for including only windows where all pixels had the same land cover was to avoid a situation where changes in variability between pixels was caused by differences in phenology for different land cover types rather than bark beetle attacks.

Figure 2.
A 4×4 pixel window in blue (attacked window 4×4 pixels #234) with the number of attacked trees per pixel (A) and land cover (B). This window had a total of 18 damaged trees spread over 4 pixels and located on land cover class 112 (spruce dominated forest), fulfilling both conditions for attacked windows.
Healthy windows are windows for which all pixels are healthy and have the same land cover type (Fig. 3); only land cover with coniferous forests were included. The reason for calculating CV instead of working with standard deviation was to compensate for differences in phenology with generally higher VI, and hence, also higher standard deviation during the main growing season. Mean of µ and CV were calculated for attacked and healthy windows respectively, and for all window sizes, VIs and time periods. All calculations were done in Python.

Data and Software Availability
The authors are not allowed to share the bark beetle damage data. The Python scripts developed are available upon request and the land cover data are freely available.

Comparison between vegetation indices
CCI seems to be the most promising VI for early detection of bark beetle attacks with a rather stable mean CV for healthy windows during the summer and fall of 2018 except for a single peak in early Augusts that is most likely noise. The attacked windows, on the other hand, show an increase in CV during the summer (Fig. 4). The main swarming was the second week of May, i.e. the trees were most likely attacked that week, and the results indicate that the CV starts increasing for attacked pixels in the later part of July. For NDVI, the results indicate that there is a slight increase in CV during the later part of the summer but there is a sharp increase in CV also for healthy windows from mid-July. For NDWI, CV increases for both attacked and healthy windows from early summer and for NDRS, CV stays rather stable during the summer and for both attacked and healthy windows CV increase later in the summer. Hence, only CCI is included in the following results. Figure 4. Time-series of mean CV for CCI, NDVI, NDWI and NDRS for attacked and healthy windows during 2018. Window size is 4×4 pixels. The raw data (grey) are rather noisy and were hence smoothed with a 5-point rolling median to show the dominating pattern (black).

Influence of window size on CV
The time-series of mean CV for CCI behaves similarly for different window sizes with an increase from late July, and small differences between window sizes except for a period in August with high noise in the data (Fig. 5). The Pearson correlation between the attacked 3×3 and 4×4 pixel windows is 0.81, between the attacked 3×3 and 5×5 pixel windows 0.68 and between the attacked 4×4 and 5×5 pixel windows 0.67. Before computing the Pearson correlation, CV values that deviated more than 1 standard deviation from the mean of the time-series, for each window size respectively, were considered outliers and removed. The outliers included both very high and very low deviation in CV for a single time-period (Fig. 5); such deviations are most likely due to noise.

Influence of attack intensity on CV
The number of attacked trees inside a window (here referred to as attack intensity) influences the change in mean CV (Fig. 6). For all attack intensities, CV starts to increase in late July but CV is highest for windows with eight or fewer attacked trees and lower for windows with 16 or more attacked trees. The increase in CV is also lower for windows with at least 16 attacked trees. When including all windows with attacked trees the pattern is similar to windows with eight or fewer threes but with slightly lower CV.

Discussion
This study indicates that in-between pixel variability can be used to detect bark beetle attacks with CCI as the more promising VI. One reason for the better performance of CCI compared to NDRS and NDWI, which have performed well in earlier studies, might be that both NDRS and NDWI include a SWIR band that originally has a spatial resolution of 20×20 m. That means that even if data is resampled to 10×10 m the originally coarser resolution will result in lower variability.
The results indicate that the window size has little influence on the variability. With a smaller window, a single pixel with attacked trees will have a stronger influence on the variability compared to a larger window. For scattered attacked trees it is more likely to have a larger number of attacked pixels with larger window sizes. It is also important to note that a larger window size will result in fewer windows, due to the limitation of one land cover type per window. Hence, we did not test with larger window sizes. In this study, the number of attacked windows were 709 for 3×3 pixel windows, 668 for 4×4 pixel windows and 588 for 5×5 pixel windows. The number of healthy windows were larger.
The results show that a lower number of attacked trees (<=8) (less intense attacks) results in higher variability and change in CV compared to windows with a larger number (>=16) of attacked trees. This might seem unexpected but a possible explanation is that with a higher number of attacked trees it is likely that a majority of the pixels in a window are attacked, and hence, the variability between pixels will be lower compared to a situation when only a few pixels are attacked.
To conclude, the results of this study indicate that inbetween pixel variability provides information that could aid in early detection of bark beetle attacked trees. We did not compare the method with pure pixel-based methods since we see the method as a way of increasing accuracy in combination with pure pixel-based methods. When developing methods for early detection it is likely that variability can be used in combination with pure pixelbased methods to include both VI values and variability in a classifier. If forest stand data are available it would also be possible to further develop the method to work with variability within forest stands rather than regular windows of fixed sizes.