1 Introduction

The counting of forest fires has been increasing since that to reach a maximum of 206 fires in 2016 from a minimum of 83 fires in 1984; corresponding burned areas were 6938.8 ha and 2198.8 ha, respectively (Christopoulou et al. 2014; Kazanis and Arianoutsou 1996; Paula et al. 2009). In Greece, almost all forest lands, a total of 2,615,000 ha, have forest fire problems in the fire season. In recent years, there has been a substantial increase in the areas burned by wildfires, related to weather conditions (Koutsias et al. 2013) especially to arson, the total number of the fire incidents varying from a minimum of 968 in 2012 to a maximum of 1443 in 2018, with corresponding areas burned of 19,613 ha and 105,450 ha, respectively (Elhag et al. 2018; Turco et al. 2016).

Crete with its favorable climatic conditions in the summer (high temperature, low moisture, long drought period) and the inflammable vegetation (phryganics, thorny bushes, evergreen broad-leaf plants) has a high fire risk (Duguy et al. 2012; Riva et al. 2017). Early study of Papanastasis (1980) has gathered available information’s on fuel attributes conducted from the phrygana communities in Crete related to fire risk. Thus, the phryganic ecosystem is intended to accumulate fuel that became a self-indulgent if remains unburned. Starting from 1980 and the number of forest fires increases, many occurred due to natural causes and some were prescribed fires (Kavouras et al. 2012; Palaiologou et al. 2018).

The problem of appraising damage caused by forest fires has been of increasing interest since the appearance of forest fires. Mills et al. (1987) described the basis of appraising forest fire damages. Agee (1998) has been concerned specifically with the value of timber in the context of deciding about forest fire suppression. Shackleton et al. (2011) and Palaiologou et al. (2018) used decision analysis to evaluate the fire hazard effects of timber harvesting.

Mapping the spatial extent of burned areas is essential to evaluate ecological impacts and economical losses, to monitor Land Use and Land Cover (LULC) change, and nowadays to model climatic influences of the burning of biomasses (Elhag and Boteva 2016b; Viedma et al. 2017). Recently, a new approach based on the eXtreme Gradient Boosting (XGBoost) feature selection and Random Forest classification was conducted by Abdullah et al. (2019), which proved to have better LULC accuracies.

Optical sensors use the emission properties of materials to characterize land surface objects including forest fires, thus the forest fires were detected from space using the unique spectra characteristics of wood (Souza Jr et al. 2005; Waigl et al. 2019). In the Infrared region of the Electromagnetic radiation (EMR) spectrum, burned soils absorb the incoming energy, whereas the surrounding land features including vegetation are comparatively high reflective. Based on this unique physical principle of reflectance, several studies had been conducted to map burned areas from space (Robinson 1991; Sharma et al. 2017; Töreyin et al. 2007).

According to Coppin and Bauer (1996) and Lausch et al. (2017), Remote Sensing applications in forest fire detection and burned areas assessment has been comprehensively developed using the advantages of the high spatial resolution and multispectral bands. These were the opportunity of large-area observation in a single image, the ability to have consecutive data and data already in digital format, and the wide range of wavelength a satellite sensor can cover (Elhag 2017; Korchenko et al. 2019). Different techniques of remote sensing have been tested, including supervised and unsupervised classification (Cihlar et al. 1998; Erbek et al. 2004; Tuia et al. 2011), vegetation indices including Normalized Difference Vegetation Index (NDVI) according to Aldhebiani et al. (2018), Goward et al. (1991), Yuan and Bauer (2007), Intensity-Hue-Saturation (IHS) transformation (Carper et al. 1990; Choi 2006; Leung et al. 2013), logistic regression modeling and Principal Component Analysis (PCA) according to Rodarmel and Shan (2002), and Fauvel et al. (2009).

Karhunen–Loeve transformation and its additive models of Abbas and Fahmy (1992) and Kouassi et al. (2001) is a well-known dimensionality reduction technique, maybe leads to a description of multidimensional data in which the axis variable are uncorrelated, with the first variable (or component) containing most of the variance of the original data set (Epstein et al. 1992; Liu 1999) and the succeeding components containing decreasing proportions of data scatter (Gastpar et al. 2006; Siljestrom Ribed and Moreno López 1995; Singh and Harrison 1985). The data decorrelation produced in this process is extremely significant in change detection analysis in multi-temporal Landsat multi-spectral image data (Elhag 2016; Li et al. 2013). This specific method is known as multi-temporal PCA and discriminates the differences of the burnt areas using the pre and post-fire differences of the area of interest (Lanorte et al. 2015; Singh and Harrison 1985; Waigl et al. 2019).

There are several methods to be used to map burned areas, mostly they are conventional methods based on the extensive field visits and the use of aerial photography in Greece. Meanwhile, the method that is going to be applied in this research paper is based on the use of remote sensing indices for forest fire mapping. The Principle Component Analysis (PCA) tends to envisage the linear transformation of a single set of variables in an attempt to maximize the projected variance of the dataset under investigation (Jolliffe et al. 2003; Lanorte et al. 2015). To this end and within the temporal image’s dataset, areas with no significant changes are estimated to be highly correlated. Meanwhile, substantially changed areas were estimated to be less correlated (Fairbanks and McGwire 2004; Su et al. 2016). The natural vegetation cover as well as agricultural lands were investigated under the PCA concepts by Volpi et al. (2015) to identify the key feature affecting the behavior of the vegetation cover.The current work aims to map recently burned areas in Thasos using satellite images and to assess the accuracy of mapping through Remote Sensing. The schematic tasks are: to explore the effectiveness of multi-temporal PCA as an image enhancement technique on Landsat-8 data to map the burnt area, to define the best band combination to employ selective band multi-temporal PCA of Landsat-8 data with the purpose of mapping the burnt area and to assess the accuracy of the map produced and to compare it with the official data provided by the National Forest Service.

2 Materials and Methods

2.1 Study Area Description

The study area is located on the island of Thasos, the most northern island of Greece. It belongs to the prefecture of Kavala, Macedonia, and it extends from 24o30′ to 24o48′ E and from 40o33′ to 40o49′ N. Its area is 399 sq. km, while its perimeter is approximately 102 km (Fig. 1). It has a volcanic origin, limestone and marbles cover its crystalline base, so the island is an important source of the latter (Elhag and Alshamsi 2019). The topographic characteristic of the island is mountainous, and the maximum height is 1217 m. over sea level. The climate of Thasos is the typical Mediterranean, characterized by hot, dry, and sunny summers and cool winters (Elhag and Bahrawi 2016a). The National Forest Service, Forest Station of Thasos, holds records of fire incidents. According to these records, the biggest forest fires in this century occurred in 1928 (1500 ha), 1938 (1700 ha), 1945 (700 ha), 1984 (1669), 1985 (10,405 ha), 1989 (8401 ha) and 2000 (187 ha). However, in the period between 2010 and 2020, there was an overwhelming increase in the number of fires and surface area burned. Indeed, the largest fires in the last century occurred in this period and resulted in the loss of about 20,000 ha of Pinus brutia and Pinus nigra forests (Elhag and Boteva 2020). The area affected constituted more than half the size of Thasos Island.

Fig. 1
figure 1

The location Thasos Island at the Aegean Sea, Greece

Full size image

2.2 Remote Sensing Data

The geomorphology of the island, the large extent of the 2016 and 2018 fires, the land cover types, as well as the existence of water bodies render the location of the case study an ideal site for operational burned area mapping (Sakellariou et al. 2019). The two fires that are depicted are between 2016 and 2018. The fire of 2016 was a mixed crown and surface fire that burned for 8 days, while the fire of 2018 was a mainly crown fire that destroyed 11.870 ha of different land types (Sakellariou et al. 2019; Tampekis et al. 2015).

The data obtained for this study consisted of two satellite images with less than 5% cloud coverage, as well as a topographic map and the official fire perimeters published by the Forest Service and the ‘Roads and Coastlines’ digital map of Thasos. The two satellite images were acquired from the Landsat-8 platform, one collected 6 days after the first fire (4 August 2016) and the other after the second fire had burned out (24 October 2018).

Landsat-8 is fashioned to generate 11 spectral bands in total, 9 bands known as the Operational Land Imager (OLI), and 2 bands known as the Thermal Infrared Sensor (TIRS). The OLI bands (1–9) are registered in the Visible, Infrared, and Shortwave Infrared, the bands are composed of 30-m spatial resolution bands except band 8 (Panchromatic band, 15-m) and Band 7 (SWIR, 60-m). Additionally, the TIRS bands (bands 10 and 11) are composed of 100-m spatial resolution and it represents the Long-wavelength infrared (LWIR) bands (Table 1).

Table 1 Landsat-8 Operational Land Imager (OLI), spectra and spatial characteristics
Full size table

2.3 Conceptual Framework

Before image classification, preprocessing of remotely sensed data is required. The two major techniques used in preprocessing are radiometric and geometric corrections. Radiometric error results in atmospheric attenuation or noise, which is a result of light scattering and absorption as it travels through the Earth's atmosphere (Boers et al. 1996), while geometric distortions occur because the imagery is representing the curved surface of the Earth in two dimensions (Singh 1989). Radiometric corrections in this study were not performed since they are important only when the objective is to detect very subtle changes. Besides, most land cover-related remote sensing investigations have ignored the atmospheric correction problem, because the signals from the objects being studied are strong enough that they can be detected despite atmospheric attenuation (Boers et al. 1996).

Furthermore, the lack of time did not allow the application of such techniques. Accordingly, an assumption that differences in the reflectance are due to changes in features and land cover changes and not related to the radiometric errors was made (Elhag and Alshamsi 2019). Therefore, the only geometric correction was performed. To enable change detection to be analyzed from the satellite imagery, the data must be co-registered and preferably matched to a map projection system (Vogelmann et al. 2001).

Image classification for the identification of different classes related to the Land Cover of the study area, a supervised classification was performed, based on the two satellite images, and previous knowledge of the terrain. For the Landsat-8 images, all the seven bands in a color composite were used, the band combination was 4-3-2, which corresponds respectively to the Near IR, Red, and Green bands of the sensor. This combination gives the best visual contrast between the different land cover categories. Training areas were located for each category using the Area Of Interest (AOI) Tools, and their spectral signatures were defined (Elhag 2017). Eight discretional classes conducted from the early image, in addition to the ninth class (burnt areas) conducted from the late image were defined based on the pixel spectral signature according to Chaudhry et al. (2006) Elhag and Boteva, (2016a) using the Maximum likelihood classifier of Ritchie et al. (2018); Robert and Gene Hwang (1996).

This parametric algorithm has been the most popular for the classification of remote sensing imagery (Xu et al. 2005). It is based on the ranges of values within the training data to define regions within a multidimensional data space, and the unclassified pixels that fall within the regions defined by the training data are assigned to the appropriate categories (Skidmore 1989). Casasent and Neiberg (1995) mentioned that the other classification classifiers are generally comparable to the Maximum Likelihood classification if the only objective is to classify targets at the "macro" level. The Likelihood classification is performed according to the following equation:

$$ g_{i} \left( x \right) = 1np\left( {\omega_{i} } \right) - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$2$}} 1n\left| {\mathop \sum \nolimits i} \right| - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$2$}} \left( {x - m_{i} } \right)^{T } \mathop \sum \limits_{i}^{ - 1} \left( {x - \mu_{i} } \right), $$
(1)

where i is the class, x is the n-dimensional data (where n is the number of bands), p(ωi)  is the probability that class ωi occurs in the image and is assumed the same for all classes, |Σi|= determinant of the covariance matrix of the data in class ωi, Σi−1 is the its inverse matrix, mi is the mean vector.

Correspondingly, the accuracy of the two resulting classifications was assessed by using a contingency error matrix and a Kappa coefficient. Error matrices compare, on a category-by-category basis, the relationship between known reference data (76 ground truth data points collected on the field) and the corresponding results of automated classification (Lillesand et al. 2014). The error or confusion matrix is a very effective way to represent map accuracy as well as other accuracy measures, such as overall accuracy, producer's and user's accuracy (Congalton 1991; Congalton and Mead 1983). The overall accuracy is computed by dividing the total number of correctly classified pixels (the sum of the elements along the diagonal) by the total number of reference pixels (Lillesand et al. 2014).

The user's accuracy tells the user of the map if a given class on the map corresponds to the same class in the ground, while the producer's accuracy informs the analyst how each class was correctly classified (Casasent and Neiberg 1995). However, these procedures only indicate how the classification strategy being employed works well in the training areas, and nothing more (Lillesand et al. 2014). For this reason, proposed applications of other techniques as a mean of improving the interpretation of the error matrix among them, the Kappa coefficient were developed. The computation of this coefficient is used to determine whether the results presented in the error matrix are significantly better than a random result (Green and Congalton 2004).

Producer’s accuracy is calculated as follows:

$$ {\text{Producer accuracy }} = \frac{{C_{aa} }}{{C_{*a} }} \times 100\% , $$
(2)

where, \(C_{aa}\) is an element at position ath row and ath column, \(C_{*a}\) is column sums, User’s accuracy is calculated as follows:

$$ {\text{User accuracy }} = { }\frac{{C_{ii} }}{{C_{i*} }} \times 100\% . $$
(3)

\(C_{ii}\) is an element at position ath row and ath column, \(C_{i*}\) is row sums.

Overall accuracy is calculated as follows:

$$ {\text{User accuracy }} = { }\frac{{\mathop \sum \nolimits_{a = 1}^{U} C_{aa} }}{Q} \times 100\% , $$
(4)

where Q and U is the total number of pixels and classes, respectively.

Matching of user’s and producer’s accuracies delivers accurateness to the classification and assures a robust liability of the implemented accuracy assessment (Congalton 1991).

Khat statistics is a second measure accuracy agreement. This measure of agreement is based on Congalton et al. (1983) findings. Khat was calculated using the following equation:

$$ {\rm K}_{{{\text{hat}}}} = \frac{{N \cdot \sum\limits_{i = 1}^{r} {x_{ii} } - \sum\limits_{i = 1}^{r} {\left( {x_{ij} \cdot x_{ji} } \right)} }}{{N^{2} - \sum\limits_{i = 1}^{r} {\left( {x_{ij} \cdot x_{ji} } \right)} }}, $$
(5)

where r is the number of rows in the error matrix, xii is the number of observations in row i and column i (the diagonal cells), xi+ is the total observations of row i, x+I is the total observations of column i, N is the total of observations in the matrix.

PCA is a multivariate statistical method in which dataset axes are transformed into principal axes, or components, that maximize dataset inconsistencies (Lanorte et al. 2015). The first step was to derive a PCA transformation from the corrected 2018 image. After combining the different components, or bands, of the outcome of this transformation, band #5 is identified as the one with the more contrast relevant to the burned areas of 2018. The problem though with this product is that it also includes the area that was burned through the 2016 fire. To apply the multi-temporal image analysis from the two initial images, the layer stack technique was used and the three bands from each image were selected and comprised in the new image. This analysis comes from a band selection that has been applied in many scholarly works and it is well documented (Chuvieco et al. 1997). Following Monahan (2000), PCA fundamental equations are:

First vector w(1) should be answered as follows:

$$ w_{\left( 1 \right)} = \arg \max_{w = 1} \left\{ {\mathop \sum \limits_{i} \left( {t_{1} } \right)_{\left( i \right)}^{2} } \right\} = \arg \max_{w = 1} \left\{ {\mathop \sum \limits_{i} \left( {x_{i} . w} \right)^{2} } \right\} . $$
(6)

The matrix form of the above equation gives the following:

$$ w_{\left( 1 \right)} = \arg \max_{w = 1} \left\{ {Xw^{2} } \right\} = \arg \max_{w = 1} \left\{ {w^{{\text{T}}} X^{{\text{T}}} Xw} \right\}. $$
(7)

w(1) should be answered as follows:

$$ w_{\left( 1 \right)} = \arg \max \left\{ {\frac{{w^{{\text{T}}} X^{{\text{T}}} Xw}}{{w^{{\text{T}}} w}} } \right\}. $$
(8)

Originated w(1) suggests that first component of a data vector x(i) can then be expressed as a score of t1(i) = x(i) ⋅ w(1) in the transformed coordinates, or as the corresponding vector in the original variables, (x(i) ⋅ w(1)) w(1).

NDVI image differencing change detection is the result of subtracting the NDVI ratios for two dates of imagery. This technique has been used to describe vegetation dynamics and changes in vegetation cover (Singh 1989). In this study, the normalized difference vegetation index images were calculated from Landsat-8 images. For the extraction and subtraction of the NDVI values, an assumption was made: the spectral and temporal differences between the two sensors produce minor differences. This assumption was considered reasonable in a similar study of Aldhebiani et al. (2018) since the objective was not to perform an accurate assessment but a relative evaluation of the difference that occurred between the two dates (2016 and 2018). From the two images, NDVI was performed according to the formula:

$$ {\text{NDVI = }}\left( {{\text{NIR}} - R} \right)/({\text{NIR}} + R), $$
(9)

where NIR, R is the Near Infrared and the Red bands of Landsat-8.

The two resulting NDVI images were subtracted from the Landsat-8 temporal data set, and the resulting image was scaled to assign to the study area classes of change. The threshold was based on visual interpretation and NDVI values. Lillesand et al. (2014), mentioned that immediate visual feedback on the suitability of a given threshold can be obtained from most image analysis system. The flowchart of the adopted methodology is illustrated in Fig. 2.

Fig. 2
figure 2

Schematic flowchart of the adopted methodology, different shapes mean different procedural steps, inputs (hexagonal), outputs (square), and conditional processors (circle)

Full size image

The current methodology is based on the integration of the PCA and the use of the NDVI values to assess and to map the burned areas in Thasos Island. The PCA was implemented as a data transformation technique to comprehensively envisage areas with multi-temporal changes. The multi-temporal images main segments are associated with continual Land cover types, areas of substantial changes will be addressed in the PCA components. Specifically, each successive component encompasses less of the total multi-temporal dataset variance. Consequently, the integration of the NDVI shall emphasize the variances of small areas (Lanorte et al. 2015; Lasaponara 2006).

3 Results and Discussion

Seventeen GCP was created to register the image to the topographic map. The topographic map is projected according to Transverse Mercator / Spheroid GRS 1980 / Datum EGSA87. The accuracy achieved was lower than one pixel (RMS error = 0.39) and the resampling method applied was the nearest neighbor. The second step is to do an image-to-image rectification using the 2016 image as a reference for the 2018 image. The accuracy, in this case, was calculated in meters (RMS error = 11.44), which is lower than a pixel (1 pixel ~ 30 m).

From the two satellite images Landsat-8, two supervised classification maps were produced. Eight classes were derived after the supervised classification of each image. Tables 2 and 3 show the confusion matrices, which summarize the agreement and confusion of the classified images with the reference data from the ground truth. Pixels that were used for testing accuracy are located along the diagonal of the error matrix, while pixels misclassified are represented along with the non-diagonal elements of the error matrix. Table 4 shows the total accuracy and Kappa coefficient of all the two supervised classifications. The result showed that the Kappa coefficient was 0.92 and 0.94 for the temporal image analysis, respectively.

Table 2 Error matrix of the early Landsat supervised classification
Full size table
Table 3 Error matrix of the late Landsat supervised classification
Full size table
Table 4 Total accuracy and Kappa Coefficient of supervised classification maps
Full size table

Where OG is the Olive Groves, NG is Natural Grassland, SV is Sparsely Vegetated Areas, MF is Mixed Forest, ES is Mineral Extraction Sites, UA is Urban Areas and SM is Salt Marshes. The eighth class is the surrounding water bodies (sea surface) and it did not interfere or overlap with any other conducted classes due to its spectral behavior (Elhag 2017). The night class (Burnt Area) appeared only on the late acquisition of the Landsat-8 image (2018) because of the forest fire occurrence (see Fig. 3).

Fig. 3
figure 3

The Maximum Likelihood classification image of the study area acquired in 2016 (up) and 2018 (down)

Full size image

The classification results generate many limitations in vegetation classification and mapping. The problem of the effect of the heterogeneity and fragmentation of Mediterranean landscapes was pointed out by Modugno et al. (2016). They had mentioned that mixing the different life forms and vegetation patterns require an understanding of the scale dependency of the vegetation classes particularly in Mediterranean countries, since there are substantial mixing classes of vegetation. The typical pattern of Mediterranean vegetation is the difference in species composition, size, and life form between north- and south-facing slopes (Rishmawi and Gitas 2001).

The following task was to apply PCA on a multi-temporal image from the two initial images. The first step was to derive a PCA transformation from the corrected 2018 image (Fig. 4). After combining the different components, or bands, of the outcome of this transformation, band four (Table 5) is identified as the one with the more contrast relevant to the burned area of 2018. The problem though with this product is that it also includes the area that was burned through the 2016 fire.

Fig. 4
figure 4

PCA transformation of the late acquired Landsat image (2018)

Full size image
Table 5 PCA transformation contrast values of the 2018 Landsat image
Full size table

The non-burned pixels represented water, urban fabric, bare surfaces, vegetation, and shade. Burned areas were assigned the value zero, whereas non-burned areas were assigned the value one to create the dependent variable in the modelling process. The PCA stacking was structured using the radiometrically corrected bands along with transformed values arising from multispectral transformations. Following the construction of the equations, the performance of each was evaluated by calculating the percentage of correct classified observations (burned, unburned and overall). The models with the better performance were applied to the entire satellite data to map the burned areas.

From the previous PCA, it was known that the bands which explain better the burned areas are three, four, five, and seven, so these bands were preferred for the new investigation. The results showed the burned area with black tones in the PC5 (Table 6). In this case, there is no confusion with the old burned area of 2016, which remains in white tones, but with the clouds. The loadings showed that the most important bands for PC5 are bands 5 and 6 from both images. Once again, the NIR and the MIR seem to be the most decisive regions of the spectrum for discriminating burned areas.

Table 6 PCA transformation contrast values of the temporal Landsat images
Full size table

Following another concept was examined. NDVI was derived for both images and included as a different layer through layer stack (Table 7). The grayscale interpretation of band five can be identified in black tone as the burned areas of 2016 (Fig. 5) and on the contrary, the burned area of 2018 is easily recognized because of its white color by the interpretation of band six (Fig. 6). The main threshold to differentiate between the burned and the non-burned areas is the estimated variances obtained from the PCA (Lanorte et al. 2015).

Table 7 The PCA of the temporal Landsat images with the NDVI understanding
Full size table
Fig. 5
figure 5

Multi-temporal PCA of Landsat-8 image Band-5 acquired in 2016

Full size image
Fig. 6
figure 6

Multi-temporal PCA of Landsat-8 image Band-6 acquired in 2018

Full size image

Stacking the PC-4, PC-5, and PC-6 from Table 6 to visualize the burned areas in color (Fig. 7), burned areas appear precisely in white concerning NDVI changing and PCA analysis. The burned area estimated in this last analysis is more helpful since it not only discriminates the difference between the two burned areas but also using a combination of the sixth and the fifth band the clouds over the area of interest can be discriminated against.

Fig. 7
figure 7

The burned areas appear in green of the temporal Landsat images

Full size image

The NDVI differencing image resulted from the subtraction of the late image (2018) from the early image (2016) demonstrated in a grayscale single band image (Fig. 8), and based on the values of the NDVI differencing image and a visual interpretation, three categories of change were derived: positive change, moderate and negative change. According to Lillesand et al. (2014), the analyst can obtain immediate visual feedback on the suitability of a given threshold. The numeric results of Landsat-8 temporal images, as well as the values of NDVI differencing, are reported in Table 8. Furthermore, areas of change in terms of vegetation (indicated by the biomass), where negative values with dark color correspond to vegetation regression (negative change), and positive values with bright colors correspond to vegetation evolution (positive change).

Fig. 8
figure 8

The burned area mapping using the NDVI differencing of Landsat-8 images

Full size image
Table 8 Statistics of Landsat-8 temporal NDVIs differencing
Full size table

Spatial variations in habitat conditions and the effect of disturbance caused by forest fires can be detected by comparing NDVI maps and may provide the basis for an early warning of degradation (Lanfredi et al. 2003). According to Hall et al. (2016), burned areas can be best detected during wet periods as the differences between burned and non-burned areas stand out most clearly. This justifies the choice of the date of image acquisition in this study. Furthermore, the choice of these two dates was favorable since it reduces the problem related to vegetation phenology differences (Van Leeuwen 2008).

The comparison of the classification with the burned area perimeter using 256 points is reported in Table 9. The classification error matrix shows the commission and emission errors in the classification. Of the 128 pixels predicted by the hue component to be burned, 11 pixels were erroneously classified as burned (commission error). The landcover types confused with the burned areas were bare/low vegetated areas and coastal areas. The emission error was higher than the commission error where 20 pixels were erroneously classified as unburned. It was noticed that slightly or moderately burned areas were classified as unburned. The overall classification accuracy was estimated to be 87.9%. It is interesting to note that topographically shaded areas were not confused with burned areas.

Table 9 Burned area accuracy assessment
Full size table

From the results of the NDVI subtraction of late from early acquisitions, the mean value was found—0.036. This means a priori a regression in the vegetation cover. Jimenez-Gonzalez et al. (2016) confirmed that the brighter the pixel is, the greater the amount of vegetation matter. From the image of NDVI differencing (excluding the water body), dark pixels with low NDVI values belong to areas where high regression of vegetation occurred. On another hand, bright pixels represent areas of vegetation evolution which unfortunately does not concern semi-natural vegetation, as a result, the soil becomes very sensitive to degradation (Bahrawi et al. 2016).

Using the image from 2016, the conclusion can be drawn as the burned areas were not included in the forest fire estimation because of the limestone quarries coverage and not vegetation coverage even before the fire of 2018 (Christopoulou et al. 2014). Having this in mind, the accuracy might even be greater than calculated. Another reason for this difference can be that the spatial resolution of the image is not high enough to differentiate vegetation when it is over on areas with high reflection like bare rock surfaces (Koutsias et al. 2013; Turco et al. 2016).

To better understand the changes that occurred between the two dates of image acquisition (2016 and 2018), especially in terms of land degradation caused by forest fires and its impact on vegetation, change detection techniques are very helpful. Vegetation change due to forest fires will be followed by a high decrease in primary production and biomass and the best example is Xilokastron in Greece, where 20% of the vegetation cover was mainly associated with the destruction of permanent shrublands and woodlands, which were subject to land degradation (Palaiologou et al. 2018).

4 Conclusions

Forest fires are considered as an integral mechanism of the ecology of the Mediterranean ecosystems and have defined, through thousands of years, the main characteristics of the Mediterranean basin forests. However, throughout the last decades the occurrence of these fires has increased, making their extent to be crucial and causing a big problem mainly because of the intense human influence, and the severe urbanization mainly on islands and remote rural areas. These have a large effect on erosion and are related to global climate change.

High spatial resolution images, provided by satellites as Landsat-8, can be used in cooperation with powerful digital image processing techniques, to present accurate results about burned areas on the Earth’s surface. The results of this research confirm the capability to precisely delineate fire scars using remote sensing science. Images can be greatly enhanced by statistic-based techniques like PCA and in the case of multi-temporal PCA; the outcome can be interpreted to discriminate burned areas at different moments in time. The bands that have been proven most useful in this work are bands four and five. NDVI on the other hand has proven to influence the image only in a way to detect the clouds and the shadows and distinguish them from the burnt area.

Finally, the approximation of the burned area using the adopted methodology reaches a high accuracy percentage of 84.61%. Multi-temporal PCA separates the two different burned areas of 2016 and 2018 although it underestimates the burned area in the middle of the island which is covered by bare rock. It would be very interesting to examine the prospect of enhancing the results by combining other methods as HIS or even to try to include band seven in the incoming research.