Classifying the Degree of Bark Beetle-Induced Damage on Fir ( Abies mariesii ) Forests, from UAV-Acquired RGB Images

: Bark beetle outbreaks are responsible for the loss of large areas of forests and in recent years they appear to be increasing in frequency and magnitude as a result of climate change. The aim of this study is to develop a new standardized methodology for the automatic detection of the degree of damage on single ﬁr trees caused by bark beetle attacks using a simple GIS-based model. The classiﬁcation approach is based on the degree of tree canopy defoliation observed (white pixels) in the UAV-acquired very high resolution RGB orthophotos. We deﬁned six degrees (categories) of damage (healthy, four infested levels and dead) based on the ratio of white pixel to the total number of pixels of a given tree canopy. Category 1: <2.5% (no defoliation); Category 2: 2.5–10% (very low defoliation); Category 3: 10–25% (low defoliation); Category 4: 25–50% (medium defoliation); Category 5: 50–75% (high defoliation), and ﬁnally Category 6: >75% (dead). The deﬁnition of “white pixel” is crucial, since light conditions during image acquisition drastically affect pixel values. Thus, whiteness was deﬁned as the ratio of red pixel value to the blue pixel value of every single pixel in relation to the ratio of the mean red and mean blue value of the whole orthomosaic. The results show that in an area of 4 ha, out of the 1376 trees, 277 were healthy, 948 were infested (Cat 2, 628; Cat 3, 244; Cat 4, 64; Cat 5, 12), and 151 were dead (Cat 6). The validation led to an average precision of 62%, with Cat 1 and Cat 6 reaching a precision of 73% and 94%, respectively.


Introduction
Bark beetle outbreaks are responsible for the loss of large areas of forests all over the world with irreversible changes on the landscape and significant economical loses. Bark beetle outbreaks reduce wood quality and production and influence the water/nutrient/carbon cycle, affecting ecosystem biodiversity [1,2]. In recent years, such pest outbreaks in forests appear to be increasing in frequency and magnitude all over the world as a result of climate change [2]. The probability of bark beetle attacks increases after long periods of drought, storms, windthrows [3,4], heavy snow, and other extreme weather conditions that weaken trees, making them ideal breeding grounds [5][6][7]. Once they have reached high population levels, bark beetles are even able to attack and infest healthy trees, as is the case of mountain pine beetles, which can cause the mortality of apparently healthy trees over millions of hectares [8]. In Japan, as in the rest of Asia, including Eurasia in Russia, the main common bark beetle is Polygraphus proximus [9], which is a non-aggressive phloephagous reaching a size of up to 12 mm. In 2013, a large-scale outbreak of P. proximus took place in Zao Mountains in northeastern Japan that, by 2016, had already devastated hundreds of hectares of pristine fir forests [10]. Since the area affected is part of a Quasi-National Park, the application of chemicals or any harmful material is forbidden and thus, the spread of bark beetles has been left undisturbed. This in turn, offers an excellent ground to monitor the natural process of bark beetle spread, forest resilience, and the role of abiotic factors in ameliorating or enhancing its effect. There are two common approaches that are used to evaluate the effect of bark beetles on forests, fieldwork observations and satellite images (usually hyperspectral ones). The first one is limited in space and the second one, even though it covers large areas at once, is limited in its precision since the medium resolution used ranges from 10 m to 100 m as is the case for Landsat images, although high resolution data (1-10 m) have been commonly used for further supporting detail mapping of the temporal effect of bark beetles [11], with very few of these studies conducted in Asia [12]. In order to monitor the development of damage at the individual tree level and the mapping of the whole affected forest area, it is necessary to quantify the number of infested trees and classify the degree of damage of each one of them. The most frequently used indexes for mapping bark beetle disturbances are the Enhanced Wetness Difference Index (EWDI), followed by the Normalized Burn Ratio (NBR) and Disturbance Index (DI) [11]. From these indexes, classifications, such as severely infested or almost dead, provide a qualitative but not a quantitative description of the damage and they cannot be used easily for comparison with other forest stands. It is also inaccurate to point to areas of infestation in coarse resolution images (>100 m) when the areas enclosed in large pixels could still include healthy trees [2,13]. Furthermore, the classification of the degree of damage of individual trees during field observations is subject to human criteria, experience, and accessibility to the tree canopy being classified, especially in highly dense forests [14,15]. Thus, the results are usually difficult to validate and compare with other forest stands, unless the same person and the same conditions are met, which is physically impossible, especially when the limited area the assessor can cover is considered. From fieldwork surveys as well as remote sensing studies, there is no standard criteria that can be used for comparison, first, because the patterns of damage are not discussed, and second, because there is no quantification of the degree of damage of single trees in time or space. In recent years, there have been several studies using unmanned aerial vehicles (UAV)-acquired images to evaluate and monitor forest pest infestation, mainly using hyperspectral cameras [16,17] and rarely RGB cameras [10,18]. In these latter two studies, the obtained data were also used to train deep learning models for tree health status detection. In both studies, the question was whether trees were healthy, infested, or dead. UAVs have shown an immense potential to close the gap between field observations and satellite imagery [19], since its low elevation flights produce very high-resolution images (0.02-0.03 m) and can cover areas of several tens to hundreds of hectares. The very high-resolution data have renewed interest in the use of a wide variety of data processing techniques, which include computer vision, soft computing, and deep learning [10,14,[20][21][22] among others. Furthermore, pixel counting in satellite images to estimate classes in a given area [23] can be used at a micro level for counting pixels to estimate defoliation in a given tree canopy area from very-high resolution images. Thus, the aims of this study are: (1) to identify the health condition of fir trees from UAV-acquired RGB images; and (2) quantify and classify the degree of damage caused by bark beetle infestation in single fir trees inferred from the white pixels of partially or completely defoliated tree branches.

Study Site
The study area is located in Zao Mountains, a volcano in southeastern Yamagata, Japan (38 • 09 05 N, 140 • 25 03 E). The site covers an area of 25 ha with a tree density of about 200 trees/ha. The age of fir trees in the area ranges from 41 to 103 years old. The study area is divided into five sites in which the degree of insect damage along an elevation gradient of 1332 m to 1670 m ( Figure 1) is evaluated. Site 1 (3.9 ha) is in the lower elevation and site 5 (6 ha) is in the highest elevation. Site 1 and site 2 are composed of Maries fir mixed with other deciduous tree species such as: Acer spp., Fagus crenata and Sorbaria sorbifolia. However, with increasing altitude, the number of deciduous tree species decreases and fir becomes the dominant tree species. Thus, in sites 3 (5.1 ha) and 4 (4 ha), fir is the dominant tree species. In site 5 (1508-1670 m), 95% of all fir trees are dead [11].
The study area is located in Zao Mountains, a volcano in southeastern Yama Japan (38°09′05" N, 140°25′03" E). The site covers an area of 25 ha with a tree dens about 200 trees/ha. The age of fir trees in the area ranges from 41 to 103 years old study area is divided into five sites in which the degree of insect damage along an e tion gradient of 1332 m to 1670 m ( Figure 1) is evaluated. Site 1 (3.9 ha) is in the l elevation and site 5 (6 ha) is in the highest elevation. Site 1 and site 2 are compos Maries fir mixed with other deciduous tree species such as: Acer spp., Fagus crenata Sorbaria sorbifolia. However, with increasing altitude, the number of deciduous tree sp decreases and fir becomes the dominant tree species. Thus, in sites 3 (5.1 ha) and 4 (4 fir is the dominant tree species. In site 5 (1508-1670 m), 95% of all fir trees are dead [

Data Acquisition
Sets of RGB aerial photos were collected during the growing season of the years and 2021. Several flight-missions were conducted with a DJI Mavic 2 pro Hasselblad 20c camera. The UAV acquired 20-megapixel images, following routes designed o GS pro software (DJI Inc., Shenzhen, China). The camera sensor is a 1-inch CMOS fixed focal length of 10 mm and aperture f/3.2. The UAV flew at 70 m altitude from off points, nadir view, 3 m/s and shutter interval of 2 s. The photos were acquired 90% side and front overlap. The camera was set at S priority (ISO-100) with shutter s at 1/240-1/320 s in the case of strong sun and wind and at 1/80-1/120 in more favo weather conditions. This setup kept exposure values (EV) of around 0 to +0.7 provid ground sampling distance (GSD) of 1.5-2.1 cm/pixel. Other missions used a DJI Pha 4 Quadcopter. The camera sensor in this UAV is 1/2.3-inch CMOS with fixed focal le of 24 mm. The UAV flew at speed 2 m/s, with a shutter interval of 2 s. The camera w

Data Acquisition
Sets of RGB aerial photos were collected during the growing season of the years 2020 and 2021. Several flight-missions were conducted with a DJI Mavic 2 pro Hasselblad L1D-20c camera. The UAV acquired 20-megapixel images, following routes designed on DJI GS pro software (DJI Inc., Shenzhen, China). The camera sensor is a 1-inch CMOS with fixed focal length of 10 mm and aperture f/3.2. The UAV flew at 70 m altitude from take-off points, nadir view, 3 m/s and shutter interval of 2 s. The photos were acquired with 90% side and front overlap. The camera was set at S priority (ISO-100) with shutter speed at 1/240-1/320 s in the case of strong sun and wind and at 1/80-1/120 in more favorable weather conditions. This setup kept exposure values (EV) of around 0 to +0.7 providing a ground sampling distance (GSD) of 1.5-2.1 cm/pixel. Other missions used a DJI Phantom 4 Quadcopter. The camera sensor in this UAV is 1/2.3-inch CMOS with fixed focal length of 24 mm. The UAV flew at speed 2 m/s, with a shutter interval of 2 s. The camera was set in automatic mode at a shutter speed of 1/120 s, ISO-100, EV at 0. This set up resulted in a GSD of 2.6 cm/pixel. The number of RAW images of each flight mission ranged from 150 to 390 [11].

Data and Programs
The orthomosaics had a pixel size of 0.014 m to 0.02 m and were delivered in JPGformat or tiff. They had 3-4 bands, but were generalized by a program to a 3-band data. Moreover, the sizes were standardized to a size of 0.02 m time 0.02 m. The coordinate system was Tokyo_UTM_Zone_54N. The programs that were used in this study were ArcGIS Pro version 2.7.0, python 3.7, and spyder 3.3.6. All programs were run under a Microsoft Windows 10 64-bit operating system. Nineteen different orthomosaics were collected under different environmental conditions for this study (Table 1).

Definition of Tree Health Category
Based on these images, the degree of defoliation in individual trees is evaluated. The most relevant visual symptom of tree infestation from UAV-acquired images is the whitish looking leafless branches in random areas within the canopy. The percentage of white pixels increases with the degree of tree defoliation (health decline), from full green when the tree is healthy (non-infested) to full white when the tree is dead or completely defoliated ( Figure 2). Thus, the canopy area of each tree is divided into pixels of 0.02 × 0.02 m size each, of which the white pixels represent the defoliated branches in each tree. The proportion of white pixels to the total number of pixels (white and non-white) within the tree canopy is calculated and classified into six different categories, representing the degree of defoliation of any given tree: Cat 1: <2.5% (healthy, no defoliation); Cat 2: 2.5-10% (very low defoliation); Cat 3: 10-25% (low defoliation); Cat 4: 25-50% (medium defoliation); Cat 5: Thus, the canopy area of each tree is divided into pixels of 0.02 × 0.02 m size each, of which the white pixels represent the defoliated branches in each tree. The proportion of white pixels to the total number of pixels (white and non-white) within the tree canopy is calculated and classified into six different categories, representing the degree of defoliation of any given tree: Cat 1: <2.5% (healthy, no defoliation); Cat 2: 2.5-10% (very low defoliation); Cat 3: 10-25% (low defoliation); Cat 4: 25-50% (medium defoliation); Cat 5: 50-75% (high defoliation), and finally Cat 6: >75% (dead), based on the classification of [15].

Calculation of Whiteness
The whiteness of each pixel is defined by the RGB-color system and its values. The color white consists of an equal proportion of the colors red, green and blue, with relatively high color values. As these values decrease, keeping the same ratio among them, the white changes to gray and then to black. However, depending on the light conditions during the collection of the images that were used for the different orthomosaics, the values and ratios for the white color varies considerably. Thus, in order to determine white and non-white pixels within the tree canopy, we carried out the following three steps: First, a pixel is considered white if the red (R) and blue (B) pixel values are higher than the average R and B pixel values of the whole orthomosaic. By considering the average value of the orthomosaic when defining white and non-white, different orthomosaics that were made with images collected under different light conditions can be compared. The ratio of the R and B values is used, because from successive tests, it was found that blue is the decisive value when the pixel is white. When the color value approaches white, the blue value increases and reaches similar values as the red values.
Since each orthomosaic has different light characteristics, the average R and B values are different for each of the 19 orthomosaics. In the orthomosaics used in this study, the average R and B values ranged from 100 to 180 ( Figure 3). Thus, for each orthomosaic, the threshold values were different. Since, after several trial-and-error tests, the threshold values of the averaged R and B values left some white pixels out, we decided to homogeneously decrease the threshold values 15 points in each orthomosaic. As a result, the model could identify clearly the white pixels in each orthomosaic regardless of the lightning conditions at which the images for a given orthomosaic were collected.  Second, we focus on the pixels that were considered non-white in the first step. If the pixel values of R and B, are lower than the average B pixel values divided by 3 and R pixel values divided by 2, then the pixel is considered as non-white ( Figure 4). The reason for dividing R by 2 and B by 3 is that the red value is usually much higher (pixel values) than that of blue. Therefore, the dividend must be set lower, resulting in a higher absolute minimum threshold value. It is important that as many non-white pixels are included as possible leaving white pixels out. Non-white pixels include green areas (foliated branches) of the tree, as well as dark pixels that are not part of the canopy, which might otherwise be defined as one in the next step. Second, we focus on the pixels that were considered non-white in the first step. If the pixel values of R and B, are lower than the average B pixel values divided by 3 and R pixel values divided by 2, then the pixel is considered as non-white ( Figure 4). The reason for dividing R by 2 and B by 3 is that the red value is usually much higher (pixel values) than that of blue. Therefore, the dividend must be set lower, resulting in a higher absolute minimum threshold value. It is important that as many non-white pixels are included as possible leaving white pixels out. Non-white pixels include green areas (foliated branches) of the tree, as well as dark pixels that are not part of the canopy, which might otherwise be defined as one in the next step. Second, we focus on the pixels that were considered non-white in the first step. If the pixel values of R and B, are lower than the average B pixel values divided by 3 and R pixel values divided by 2, then the pixel is considered as non-white (Figure 4). The reason for dividing R by 2 and B by 3 is that the red value is usually much higher (pixel values) than that of blue. Therefore, the dividend must be set lower, resulting in a higher absolute minimum threshold value. It is important that as many non-white pixels are included as possible leaving white pixels out. Non-white pixels include green areas (foliated branches) of the tree, as well as dark pixels that are not part of the canopy, which might otherwise be defined as one in the next step. For the final step, the ratio between the R and the B pixel value (R/B) is calculated. First, we determine the average red and green (G) color values and the ratio between them. If the ratio is higher than 1, red is the dominant color, but if green is the dominant color, the ratio is lower than 1. Thus, threshold ratios must be adjusted. When the average R is higher than the average G, and the R/B ratio value of an individual pixel is lower than 1.3, then the pixel is considered as white. However, if the ratio between R and G pixel values (R/G) is smaller than 0.9, then the pixel R/B ratio should be lower than 1.2 to be considered as white. If the average R pixel value is higher than the average G pixel value, then pixels with an R/B ratio lower than 1.1 are considered as white ( Figure 5). Under all other conditions, the pixel is considered as non-white. The R/B ratio was chosen because there are large variations between these two colors. Once all pixels have been categorized as white or non-white, then the model calculates how many white pixels are part of the area of the tree. Every white pixel is assigned the value 1, while all the non-white pixels are assigned the value 0. These steps are repeated for each individual pixel of a given orthomosaic. The result is a grid shape with 1 s and 0 s and then by using steps one, two, and three, the threshold errors are also minimized. For the final step, the ratio between the R and the B pixel value (R/B) is calculated. First, we determine the average red and green (G) color values and the ratio between them. If the ratio is higher than 1, red is the dominant color, but if green is the dominant color, the ratio is lower than 1. Thus, threshold ratios must be adjusted. When the average R is higher than the average G, and the R/B ratio value of an individual pixel is lower than 1.3, then the pixel is considered as white. However, if the ratio between R and G pixel values (R/G) is smaller than 0.9, then the pixel R/B ratio should be lower than 1.2 to be considered as white. If the average R pixel value is higher than the average G pixel value, then pixels with an R/B ratio lower than 1.1 are considered as white ( Figure 5). Under all other conditions, the pixel is considered as non-white. The R/B ratio was chosen because there are large variations between these two colors. Once all pixels have been categorized as white or non-white, then the model calculates how many white pixels are part of the area of the tree. Every white pixel is assigned the value 1, while all the non-white pixels are assigned the value 0. These steps are repeated for each individual pixel of a given orthomosaic. The result is a grid shape with 1 s and 0 s and then by using steps one, two, and three, the threshold errors are also minimized.

Polygon Definition
Polygons were drawn around every single fir tree and saved in a shape file. In this study, the polygons were drawn manually. There was a careful consideration for polygons to include as much area from the tree and as little of the non-tree area as possible to make sure that each tree canopy area was correctly represented. The number of white (1) and non-white pixels (0) are counted for each polygon using Equation (1).

Proportion of whiteness =
where ∑Pixelwhite is the sum of all white pixels of the canopy area from a single tree and ∑Pixel is the number of all pixels within the canopy area of the same tree. This equation calculates the degree of damage (a proxy for defoliation) of an individual fir tree in percentage based on the proportion of white and non-white pixels (POW).
The GIS model uses three simplified steps. The first step is the creation of polygons as well as the extraction of the RGBavg values of the orthomosaic. In the next step, the model calculates the number of white or non-white pixels and in the third step, the two previous steps are combined to calculate the POW-values and classified tree health. We used the orthomosaics from site 1, 2, 3, and 5 to calibrate the model and applied it on Site 4 (year

Polygon Definition
Polygons were drawn around every single fir tree and saved in a shape file. In this study, the polygons were drawn manually. There was a careful consideration for polygons to include as much area from the tree and as little of the non-tree area as possible to make sure that each tree canopy area was correctly represented. The number of white (1) and non-white pixels (0) are counted for each polygon using Equation (1).

Proportion of whiteness
where ∑Pixel white is the sum of all white pixels of the canopy area from a single tree and ∑Pixel is the number of all pixels within the canopy area of the same tree. This equation calculates the degree of damage (a proxy for defoliation) of an individual fir tree in percentage based on the proportion of white and non-white pixels (POW). The GIS model uses three simplified steps. The first step is the creation of polygons as well as the extraction of the RGB avg values of the orthomosaic. In the next step, the model calculates the number of white or non-white pixels and in the third step, the two previous steps are combined to calculate the POW-values and classified tree health. We used the orthomosaics from site 1, 2, 3, and 5 to calibrate the model and applied it on Site 4 (year 2021). This procedure is repeated for each tree (polygon) (Figure 6). and non-white pixels (0) are counted for each polygon using Equation (1).
where ∑Pixelwhite is the sum of all white pixels of the canopy area from a single tree and ∑Pixel is the number of all pixels within the canopy area of the same tree. This equation calculates the degree of damage (a proxy for defoliation) of an individual fir tree in percentage based on the proportion of white and non-white pixels (POW).
The GIS model uses three simplified steps. The first step is the creation of polygons as well as the extraction of the RGBavg values of the orthomosaic. In the next step, the model calculates the number of white or non-white pixels and in the third step, the two previous steps are combined to calculate the POW-values and classified tree health. We used the orthomosaics from site 1, 2, 3, and 5 to calibrate the model and applied it on Site 4 (year 2021). This procedure is repeated for each tree (polygon) (Figure 6).

Results
From the orthomosaics, we were able to visualize the proportion of defoliation of trees as an indicator of tree health status but with the results of the GIS-based model, we were able to automatically calculate the degree of tree damage (health decline) into six categories. The different levels of defoliation, shown in violet in Figure 7, represent each category: 1.1 % defoliation (Figure 7a,b) corresponded to Cat 1, while defoliations of 9.2% Figure 6. Workflow of the model for tree damage classification.

Results
From the orthomosaics, we were able to visualize the proportion of defoliation of trees as an indicator of tree health status but with the results of the GIS-based model, we were able to automatically calculate the degree of tree damage (health decline) into six categories. The different levels of defoliation, shown in violet in Figure 7, represent each category: 1.1 % defoliation (Figure 7a,b) corresponded to Cat 1, while defoliations of 9.2% (Figure 7c,d), 16.4% (Figure 7e,f), 46.9% (Figure 7g,h), 59.1% (Figure 7i,j) and 98.8% (Figure 7k,l), corresponded to Cat 2-6, respectively. The white color within the canopy of each tree, in the RGB orthomosaic, is perceived by humans but can hardly be quantified visually. In comparison to any subjective method used at present, the model provided a millimetric account of all the defoliated branches in the canopy of any given tree.
In the case of a healthy tree, the maximum number of white pixels found (below the 2.5% threshold) matched the condition of the tree, as is visualized in the images and partially in the field. The same could be said about trees in Cat 6, which represented the dead trees. In this case, the number of white pixels found (above 75%), corresponded to the characteristics of a dead tree in the orthomosaic and in the field. These figures also show that the precise definition of the polygon around the tree is of the outmost importance for the model, avoiding the inclusion of areas that are not part of the tree. Thus, the more damaged trees need special attention with regard to the polygon definition, since the canopy areas (branches included) are smaller and their shape does not follow a clear pattern in comparison to trees in other categories.
It should also be noted that if the white colors (violet in the figure) outside the tree canopy ( Figure 7) are counted, the error of calculating the POW will increase considerably.
In the orthomosaic corresponding to Site 4 in our study area (Figure 8a), the average R value was 179, while the average G and B values were 202 and 170, meaning that the whole area was predominantly green. The average R and B values were below 180 but above 160, which means that when the pixel R and B values were above 145, the pixel was declared white. When the R and B values of a pixel were below 145, the next step was followed in order to eliminate darker pixels. When the R pixel value was less than the average R value divided by 2 and the blue pixel value was lower than the average B value divided by 3, the pixel was declared as non-white. When the R pixel value was lower than 89.5 (179/2) and the B pixel value was lower than 56.6 (170/3), the pixel was declared as non-white.
trees. In this case, the number of white pixels found (above 75%), corresponded to the characteristics of a dead tree in the orthomosaic and in the field. These figures also show that the precise definition of the polygon around the tree is of the outmost importance for the model, avoiding the inclusion of areas that are not part of the tree. Thus, the more damaged trees need special attention with regard to the polygon definition, since the canopy areas (branches included) are smaller and their shape does not follow a clear pattern in comparison to trees in other categories. It should also be noted that if the white colors (violet in the figure) outside the tree canopy ( Figure 7) are counted, the error of calculating the POW will increase considerably.
In the orthomosaic corresponding to Site 4 in our study area (Figure 8a), the average R value was 179, while the average G and B values were 202 and 170, meaning that the whole area was predominantly green. The average R and B values were below 180 but above 160, which means that when the pixel R and B values were above 145, the pixel was declared white. When the R and B values of a pixel were below 145, the next step was followed in order to eliminate darker pixels. When the R pixel value was less than the average R value divided by 2 and the blue pixel value was lower than the average B value divided by 3, the pixel was declared as non-white. When the R pixel value was lower than 89.5 (179/2) and the B pixel value was lower than 56.6 (170/3), the pixel was declared as non-white. When at least one pixel was above one of these calculated values, 89.5 for R or 56.6 for B, the ratio between the R and B pixel values needed to be below 1.2 (average R divided by average G is below 0.9) to be declared white. Otherwise, the pixel was considered as When at least one pixel was above one of these calculated values, 89.5 for R or 56.6 for B, the ratio between the R and B pixel values needed to be below 1.2 (average R divided by average G is below 0.9) to be declared white. Otherwise, the pixel was considered as non-white. This was repeated for all pixels in the orthomosaic. In the last step, in each of the polygons, representing the canopy of a tree (Figure 8b), the POW value was calculated from the white and non-white pixels in the polygons (Figure 8c) and finally, each individual tree damage category was determined (Figure 8d). The results of analyzing each single tree health condition in Site 4 showed that in the 4-ha area there were 1376 trees, of which 277 were healthy (Cat 1) while 628, 244, 64, and 12 trees were found in different degrees of damage (defoliation) corresponding to Cat 2-5, respectively. Finally, the number of dead trees (Cat 6) was 151 or 11% of the total number of trees. Besides the number of trees in each of the categories, the model also provided the spatial distribution of forest decline for the year 2021.

Discussion
The first bark beetles attack in Zao Mountains was reported in 2012 and by the year 2016, almost all fir trees near the tree-line were dead (Site 5). Since then, forests in lower grounds (Sites 1-4) have not been considered as seriously affected by the spread of bark beetle infestation [24]. The lack of detailed analysis on the rate of forest decline have obscured the magnitude of infestation within the fir forest, as was pointed out by Nguyen et al. [10]. Bark beetle infestation is usually classified as green-attack stage (very early infestation), red-attack (early infestation), and grey-attack (later stage of infestation) in the wake of the infestation [25,26]. The results of our study represent the rate of forest decline as the grey-attack stage has set in and when trees start showing signs of defoliation after the severe effect of bark beetle proliferation in the phloem tissues, disrupting the translocation of water and nutrients in affected trees [27]. Our results showed that in Site 4, ranging from an altitude of 1435 m to 1495 m, only 20% of the trees of a total of 1376 trees were found healthy. Most of the trees were slightly damaged, showing a level of defoliation that extends from 2.5% to 25% (Cat 2-3) in 63% of the trees and only 6% of the trees are in the last stages of defoliation, and thus seriously damaged (Cat 4-5), while the number of dead trees (Cat 6) was 151. However, if Cat 4-5 trees are also considered in the long-term, then the number of dead trees would account for 17% of the total number of trees. Based on the distribution of tree health condition in the forest, we could observe that infested trees were distributed all over the site, regardless of the altitude or its proximity to the forest road or the ropeway line. In the northeastern areas of the study site, small pockets of heavily infested trees can be found, while most of the healthy trees were found in the southeastern area of the site. There is no direct method to validate if the more than 1000 trees that have been classified in the six tree damage categories are in agreement with the subjective evaluation in the field. Therefore, in order to test our model, the images (polygons around the canopies) of 67 trees including the six categories were selected randomly from the 4-ha site, and assessors familiar with the evaluation of insect infested forests were asked to visually classify the images for comparison with the results of our GIS-based model. The results of this comparison proved, as expected, high variability in human perception of tree infestation depending on the subjective criteria of each of the assessors ( Table 2). The comparison between both methodologies showed in general an average agreement of 62%. The highest precision was observed in Cat 1 and Cat 6 with 79% and 94%, respectively, while the lowest precision (36%) was reported for Cat 3, showing the difficulty for the assessors to distinguish the differences between small differences in tree defoliation. Similar results were found for the sensitivity and specificity values.
Our results showed that this new systematic evaluation method can be easily verified, shared, and used for comparison with other fir forests in Japan, where they extend from north to south in a transect of 700 km from center to northern Honshu Island. Thus, in comparison with other studies [10,16,19], we were able to classify different degrees of damage and provide a more thorough assessment of forest health decline.

Considerations When Using the Model
Since the critical pixels of the model are the white pixels, images collected under snow covered conditions should be avoided. Uniform light conditions in the area where pictures are being taken are ideal, avoiding as much as possible shadows covering part of the tree canopy. Preferably, cloudy conditions or late hours in the afternoon provide the best conditions for homogeneous lightning in large forested areas that are being photographed sometimes for hours. In general, the model can work relatively well under a wide range of lightning conditions. Nevertheless, care should be taken to ensure that the orthomosaic is homogeneous with respect to lighting conditions. For example, a very well-lit slope and a heavily shaded slope should not be included in the same orthomosaic. This would otherwise lead to threshold problems and misinterpretations, since an average value is used for two areas with different conditions. This is because the model uses the average values of the orthomosaic. By using the average values, different orthomosaics can be compared with each other, regardless of the different conditions, while one single large orthomosaic with many different conditions might lead to incorrect results. Different factors during the flight mission affect significantly light conditions. For example, a forest looks different at 9:00 than at 15:00 h, with seasonal light conditions and slope aspect affecting the orthomosaic too. In the area of Zao Mountains, the months with the most stable weather conditions were June and October.
Since counting pixels within the polygon (tree canopy) is the key for the precision of the model, extreme care should be taken on the precise delineation of it. If the polygon is done incorrectly, such as too small or too large, pixels will be missed or extra pixels not belonging to the tree will be included, respectively. Thus, in this study, we have chosen to establish the polygons manually for each tree, which is at present the limitation for complete automatization of the model. In the future, the automatic delineation of the polygons around the multi-shaped tree canopies in a bark beetle infested forests can be performed with deep learning. Another limitation of the GIS-based model proposed in this study is that UAV-acquired images captured tree canopies from above, which means that hidden defoliated branches cannot be detected. Accounting for those missing branches can be partially solved if we considered the different types of defoliation symptoms of infested trees. In some cases, defoliation occurs evenly on the top of the tree canopy, while in other trees defoliation is concentrated on the canopy side from top to bottom, or evenly distributed in the lower branches of the tree canopy. Sometimes, these types of infestation are also mixed. Therefore, it is crucial to consider the entire tree canopy area, and not only the tree top or partial areas of the tree, because if not we will increase the effect of 'hidden branches' and we will be ignoring the patterns of damage, especially those occurring in the lower branches.

Conclusions
The GIS-based model developed in this study estimated the defoliation of single trees from white pixels in the very high-resolution images from UAVs and used it as a proxy for the degree of damage of single trees. With this methodology, it was possible to rapidly evaluate hundreds of trees in a relatively short amount of time. The model is the first step towards the automatic classification and evaluation of pest infestation of single trees in fir forests, but it can be adapted for any type of disturbance that causes defoliation in trees. The model could also be combined with deep learning for an automatic detection of the degree of damage in individual trees. The results of the model can be ideally used as training data for deep learning models, thus, replacing the subjective classification of humans for a numerical objective calculation.
The model proposed in this study showed the potential to evaluate the spatial effect of pest infestation in fir forests covering large areas in a difficult terrain. The method that is proposed can be used by forest managers or scientists with minimum equipment capabilities.