SATELLITE-BASED HYPERSPECTRAL IMAGING AND CARTOGRAPHIC VISUALIZATION OF BARK BEETLE FOREST DAMAGE FOR THE CITY OF ČABAR

After enormous amount of ice rain during 2014 huge damage was done in forests in Croatia, especially in the city of Čabar area. Damage of forests is reflected in wide spread of bark beetle. Bark beetle damaged forest have different spectral range from healthy forest. Copernicus satellite land monitoring imagery enables distinguishing healthy from unhealthy forest. In this paper, the width of bark beetle infection spread in forests in the city of Čabar area using satellite images and semi-automatic classification will be determined.


INTRODUCTION
In the first days of February 2014 a huge amount of ice rain fell in the area of Gorski kotar.The consequence of such weather is enormous damage on forest ecosystem with instant damage reflected in broken and pulled out trees.The aftermath of broken trees is the progressed infestation of bark beetles [10].To prevent further infestation, the damaged forest should be cut down.This is a serious task because of the wide area infected.The first task is to determine which trees are infected.To do so, the easiest way is using satellite imagery.In this paper, the classification on training site located in the city of Čabar area will be presented.Municipality of Čabar is located on 28,205 hectares of which 94% are forests.[11].

Characteristics of Spruce Bark Beetles
The Republic of Croatia can be separated in three geographic regions: lowland, mountain area, and the Mediterranean.In mountain area we can find spruce, beech and fir trees.Čabar is located on 650 to 1,200 meters above sea level and therefore is placed in mountain region covered with spruce and fir forests [11].Fig. 1 shows spruce distribution in the city of Čabar area.As shown in Fig. 1, the whole area of the city of Čabar is covered with mainly spruce forest and this area is convenient for research of bark beetle damage.According to [1], in spruce forest spruce bark beetles have a great role and we can distinguish two different types of spruce bark beetles.Main characteristic of bark beetle damaged spruce forest is first yellow treetop, then red treetop, and eventually the infected tree is dead [1].Fig. 2 represents spruce forest infected with bark beetle.There are other symptoms, such as dispersal of green needles and dripping resin out of tree bark.However, in this paper the main focus will be only on color changing symptom because it is very easy to distinguish different colors on satellite images.

SEMI-AUTOMATIC CLASSIFICATION
Classification is analysis technique for remotely sensed image processing.It contains three types of methods: supervised, unsupervised and hybrid.The common classification algorithms contain the K-means, the parallelepiped, ISODATA, maximum likelihood classifier and minimum distance to mean [3].In QGIS there is available plugin called Semi-Automatic Classification Plugin for supervised classification.In this plugin there are three different classification algorithms: • minimum distance, • maximum likelihood, • spectral angle mapping.

Minimum Distance
Minimum distance classification is probably the oldest and simplest approach to pattern recognition, namely template matching.In a template matching we choose class or pattern to be recognized, such as healthy vegetation.Unknown pattern is then classified into the pattern class whose template fits best the unknown pattern.Unknown distribution is classified into the class whose distribution function is nearest (minimum distance) to the unknown distribution in terms of some predetermined distance measure [4].

Maximum Likelihood
Maximum likelihood algorithm is considered the most accurate classification scheme with high precision and accuracy [3], and because of that it is widely used for classifying remotely sensed data.Maximum likelihood classification is method for determining a known class of distributions as the maximum for a given statistic.An assumption of normality is made for the training samples.During classifications all unclassified pixels are assigned to each class based on relative probability (likelihood) of that pixel occurring within each category's probability density function [5].

Spectral Angle Mapping
Spectral image mapper is a spectral classifier that is able to determine spectral similarity between image spectra and reference spectra by calculating the angle between the spectra, treating them as vectors in a space with dimensionality equal to the number of bands used each time [6].Small angles between the two spectrums indicate high similarity, and high angles indicate low similarity [6].

DATASETS AND METHODOLOGY
The study area is located in the Republic of Croatia, in the city of Čabar area.It is located in the north part of Gorski kotar near Slovenian border.The whole area mostly consists of spruce and fir forests.According to [1], bark beetle infestation is best seen in spring and summer, where healthy spruces and fir have green treetops and the infected ones have orange/red treetops.Therefore, from the official Copernicus site [19], Sentinel-2 images were obtained for late July, precisely on July 22, 2016.It is important to choose imagery which has little or zero cloud coverage.Since only bands 2, 3, 4, 8 are needed for analysis, the other bands were dismissed.Bands 2, 3, 4 are needed for true image representation as shown in Fig. 3. Bands 3, 4 and 8 are necessary for calculating RGI and NDVI as shown in Eq. ( 1) and Eq. ( 2).The result after merging all bands into one raster image is true color image of our training area as shown in Fig. 3.  Second, the buildings (city, villages, roads) reflect differently than vegetation.Likewise, it shows that the healthy vegetation is easily detected by plain sight as it reflects in green, while infected vegetation is brown/orange/red.To confirm this theory, the red-green index is calculated.It is sensitive to conifer tree mortality according to [7], and this is shown in Fig. 4.Moreover, according to [7], normalized difference vegetation index is calculated and it is sensitive to green (healthy) vegetation as shown in Fig. 5. Red-green index is calculated from Eq. ( 1) and normalized difference vegetation index from Eq. (2). 4 3 Fig. 4 shows green areas, which represent infested forest, and yellow -red pixels represent healthy vegetation.In Fig. 5, white/very bright pixels which represent buildings (houses, roads, etc.) and green pixels can be separated.There are dark green pixels and they represent healthy forest, while light green pixels represent unhealthy forest.Since it has been determined that there are infested spruce trees in the training area, supervised classification can be done.As mentioned earlier, the maximum likelihood algorithm is currently most accurate classification method and it is going to be used on training site.First step is to manually choose reference pixels throughout the training site.Having features in training site in mind, four different macro classes have been proposed.First class contains buildings (city, village, roads, etc.).Since some of the objects are smaller than the size of one pixel, QGIS plugin joins pixels in class based on average value of pixel and on predefined threshold.Second class contains vegetation such as meadows and lawns.Third class represents healthy forest, and the forth class represents bark beetle infested/damaged forests.Fig. 6 shows spectral signature of obtained classes and shows the spectral plot and reflectance of each class.If looked at closely, it can be seen that first class differs from the other three classes and therefore classification for the first class can be expected.As it regards other three classes, they are slightly different one from another but classification should also give satisfying results.The availability of a large number of remote sensing variables has caused the difficulty to select useful variables [8].Thus, it becomes very important how to select the remote sensing variables that significantly contribute to increasing accuracy of distinguishing land use and land cover classification types [8].When the number of remote sensing variables is relatively small, simple and traditional methods, such as bar graph spectral plots and feature space plots, are usually utilized [8].If selected classes have spectral plot such as shown in Figure 6, the classification can start and should give good results.If selected classes have spectral plots that are very near each other or overlapping, classification can start but the results may not be as expected.Overlapping spectral plot of classes may occur if one of the class samples has been joined to the other class and vice versa so one should be very careful when assigning samples to classes.After it has been determined that obtained classes have satisfying spectral plot, the classification starts.Maximum likelihood algorithm in QGIS is time consuming so the final result is available in couple of hours.Later, the classification training site looks as shown in Fig. 7.

RESULTS
After successful classification a lot of training area is populated with damaged forest.Results of classification can be presented in numbers as shown in Tab. 2. Official data available online [11] shows the difference between total area of the city of Čabar in 2,631.57hectares.This could be the result of imperfect vector data acquired from web page [13].Tab. 2 and Fig. 7 show the correspondence between class number in Tab. 2 and legend entry in Fig. 7. Therefore, class number 4 represents unhealthy forest or, as shown in Fig. 7, bark beetle damage.Based on the conducted analysis, over 31% of the city of Čabar area is populated with unhealthy forest.It cannot be concluded from satellite images if the damage is from bark beetle or for some other reason.But official records of bark beetle infestation [20] reflect that around 4,000 hectares of public forests are infected with bark beetle.Official data for municipality Čabar says that 4,187 hectares or forested area is in private property [9] and it can be assumed that 10 -25% is infested with bark beetle.If these two values are added up, the result is around 5,000 hectares of infested forest.This value is different from the value given in Table 2 because not all of the damage comes from bark beetle.However, it can give a good overview on damaged forest and their location.

CONCLUSION
After huge amount of ice rain in Gorski kotar in 2014 there was an outburst of bark beetles which attack primary spruce and fir forest.To prevent further development of disease, infected trees must be cut down.However, the detection of such trees is difficult due to the large area infected.With Copernicus program being developed and its Sentinel satellite missions, it has become a lot easier to track down unhealthy forests.Using satellite images and proper tools to classify these images useful information about damaged forests can be obtained.However, it should be taken into consideration that such classification cannot always produce helpful information so one should be very careful when using such classification results.Nevertheless, if classification results are combined with results obtained on the field, quality map, which allows us to track down infested forest with higher accuracy, can be produced.

Figure 3
Figure 3 Training area -the city of Čabar area

Fig. 3
Fig.3represents several things.First, the cloud coverage is minimum and should not affect analysis.Second, the buildings (city, villages, roads) reflect differently than vegetation.Likewise, it shows that the healthy vegetation is easily detected by plain sight as it reflects in green, while infected vegetation is brown/orange/red.To confirm this theory, the red-green index is calculated.It is sensitive to conifer tree mortality according to[7], and this is shown in Fig.4.Moreover, according to[7], normalized difference vegetation index is calculated and it is sensitive to green (healthy) vegetation as shown in Fig.5.Red-green index is calculated from Eq. (1) and normalized difference vegetation index from Eq. (2).

Figure 4 Figure 5
Figure 4 Red-green index

Figure 6 Figure 7
Figure 6 Spectral signature plot

Table 1
Name, spectral range (nm) and spatial resolution (m) of corresponding

Table 2
Number of pixels per class and corresponding area