Exploring the Relationship between Burn Severity Field Data and Very High Resolution GeoEye Images: The Case of the 2011 Evros Wildfire in Greece

Monitoring post-fire vegetation response using remotely-sensed images is a top priority for post-fire management. This study investigated the potential of very-high-resolution (VHR) GeoEye images on detecting the field-measured burn severity of a forest fire that occurred in Evros (Greece) during summer 2011. To do so, we analysed the role of topographic conditions and burn severity, as measured in the field immediately after the fire (2011) and one year after (2012) using the Composite Burn Index (CBI) for explaining the post-fire vegetation response, which is measured using VHR satellite imagery. To determine this relationship, we applied redundancy analysis (RDA), which allowed us to identify which satellite variables among VHR spectral bands and Normalized Difference Vegetation Index (NDVI) can better express the post-fire vegetation response. Results demonstrated that in the first year after the fire event, variations in the post-fire vegetation dynamics can be properly detected using the GeoEye VHR data. Furthermore, results showed that remotely-sensed NDVI-based variables are able to encapsulate burn severity variability over time. Our analysis showed that, in this specific case, burn severity variations are mildly affected by the topography, while the NDVI index, as inferred from VHR data, can be successfully used to monitor the short-term post-fire dynamics of the vegetation recovery.


Introduction
Over the last decades, forest fires have become more frequent, prolonged, and severe than ever before [1,2], especially in the Mediterranean region. As such, fires have long-lasting and devastating impacts on forest ecosystems both in ecological and socio-economic terms. Therefore, it is now more important than ever to act upon enforcing sustainable fire management policies and promote forest sustainability [3].
To decide upon effective policies, national governments and stakeholders require accurate and detailed information regarding the impact of wildfires on forest ecosystems and their response after the disturbance [4][5][6]. Accurate information regarding the spatial distribution and level of damage (i.e., burn severity), as well as the post-fire vegetation recovery, is essential to quantify the impact of fire on the landscape, select and prioritize treatment measures, and plan and monitor restoration activities [7]. The overall effectiveness of post-fire management plans (restoration activities, etc.) This could be because such studies require extensive and costly field measurements for extracting the burn severity metric. Although VHR data have been mainly used for assessing burn severity, there remains a need to explore the interrelationships between remotely-sensed measurements and the physical/ecological processes that the CBI index infers [3].
So far, a limited number of studies have found statistical correlations between fire severity metrics, topography, VHR satellite data, and ecosystem responses [9,28,36,37]. This has generated the need for more studies focused on the understanding of these relationships.
The main objectives of this study are (i) to investigate the potential of the remotely-sensed post-fire vegetation response data as an expression of field-measured burn severity and topography; and (ii) to identify the best satellite variables (bands and NDVI index) able to capture the post-fire vegetation dynamics in the short-term. In this perspective this study used VHR GeoEye images of a forest fire that occurred in Evros (Greece) during summer 2011, in relationship with the CBI measured immediately after the fire (November 2011) and one year after the fire (August 2012).

Study Area
The study area is located in Evros, Greece's northernmost prefecture (40˝57 1 -41˝0 1 N and 26˝2 1 -26˝8 1 E) ( Figure 1). Evros is located in the eastern-northeastern part of Thrace and shares international borders with Turkey to the east and Bulgaria to the north. The forests in this region are mainly composed of conifers, oaks, and beeches, and they are widely known for their unique ornithological and ecological value. The climate in this region is mainly Mediterranean, although it is affected by continental conditions. Remote Sens. 2016, 8,566 3 of 15 [13,[33][34][35]. This could be because such studies require extensive and costly field measurements for extracting the burn severity metric. Although VHR data have been mainly used for assessing burn severity, there remains a need to explore the interrelationships between remotely-sensed measurements and the physical/ecological processes that the CBI index infers [3]. So far, a limited number of studies have found statistical correlations between fire severity metrics, topography, VHR satellite data, and ecosystem responses [9,28,36,37]. This has generated the need for more studies focused on the understanding of these relationships.
The main objectives of this study are (i) to investigate the potential of the remotely-sensed post-fire vegetation response data as an expression of field-measured burn severity and topography; and (ii) to identify the best satellite variables (bands and NDVI index) able to capture the post-fire vegetation dynamics in the short-term. In this perspective this study used VHR GeoEye images of a forest fire that occurred in Evros (Greece) during summer 2011, in relationship with the CBI measured immediately after the fire (November 2011) and one year after the fire (August 2012).

Study Area
The study area is located in Evros, Greece's northernmost prefecture (40°57′-41°0′N and 26°2′-26°8′E) ( Figure 1). Evros is located in the eastern-northeastern part of Thrace and shares international borders with Turkey to the east and Bulgaria to the north. The forests in this region are mainly composed of conifers, oaks, and beeches, and they are widely known for their unique ornithological and ecological value. The climate in this region is mainly Mediterranean, although it is affected by continental conditions. The area was affected by a large forest fire (5996 ha) in the summer of 2011. The fire started on 24 August 2011 and it was extinguished three days later, on 27 August 2011. Prior the the fire the area was mainly covered by coniferous (29.96%), oak (18.39%), and mixed oak-coniferous forests (27.15%). The fire is the first fire event after the reforestation interventions that took place in the mid-1970s. The predominant types of landcover at that point of time (in the mid-1970s) was makia The fire is the first fire event after the reforestation interventions that took place in the mid-1970s. The predominant types of landcover at that point of time (in the mid-1970s) was makia vegetation and sparse oak forests. The area's elevation varies from 0 m to 1650 m above sea level, whereas the elevation within the fire perimeter varies from 35 m to 433 m ( Figure 2).

Field Data
We collected two sets of data over 100 randomly-selected 60 × 60 m plots ( Figure 1). The first set of data was collected in November 2011 immediately after the fire incident, and the second one a year later, in August 2012. To select the location of our burn severity data sampling, we used a stratified random sampling scheme based on a burn severity map. Such a map was derived using the object-based image analysis (OBIA) [38] method by classifying the multispectral four-band WorldView-2 imagery into seven classes: water, healthy vegetation, bare land/other land uses (ex.agricultural areas), high severity, high-medium severity, and medium-low severity.
The whole process was fulfilled in two main phases. In the first phase, we applied a multi-resolution segmentation to the WorldView-2 image. Throughout the procedure, we segmented the image and generated image objects based on several adjustable criteria of scale, band weights, and homogeneity in color and shape. Overall, we created two segmentation levels for the WorldView-2 image. In each level, the segmentation parameters were determined through a trial-and-error procedure. The first level targeted exclusively at agricultural areas by using the land-parcel identification system (LPIS) thematic layer parameter to segment the image. The LPIS thematic dataset is a vector layer that contains information about the boundaries of all Greek agricultural areas and it was used during the segmentation process to separate the agricultural areas from its surroundings. At the second level, the area outside the LPIS boundaries was segmented again in order to derive homogeneous objects that are representative of all of the remaining classes. In the second phase, we developed a rule-based classification model to classify image objects into the burn severity classes. In our analysis, the representative features of the object classes were determined with the associated threshold value and implemented in the classification model. We regulated the membership of objects in certain classes via classification rules with fixed threshold

Field Data
We collected two sets of data over 100 randomly-selected 60ˆ60 m plots ( Figure 1). The first set of data was collected in November 2011 immediately after the fire incident, and the second one a year later, in August 2012. To select the location of our burn severity data sampling, we used a stratified random sampling scheme based on a burn severity map. Such a map was derived using the object-based image analysis (OBIA) [38] method by classifying the multispectral four-band WorldView-2 imagery into seven classes: water, healthy vegetation, bare land/other land uses (ex.agricultural areas), high severity, high-medium severity, and medium-low severity.
The whole process was fulfilled in two main phases. In the first phase, we applied a multi-resolution segmentation to the WorldView-2 image. Throughout the procedure, we segmented the image and generated image objects based on several adjustable criteria of scale, band weights, and homogeneity in color and shape. Overall, we created two segmentation levels for the WorldView-2 image. In each level, the segmentation parameters were determined through a trial-and-error procedure. The first level targeted exclusively at agricultural areas by using the land-parcel identification system (LPIS) thematic layer parameter to segment the image. The LPIS thematic dataset is a vector layer that contains information about the boundaries of all Greek agricultural areas and it was used during the segmentation process to separate the agricultural areas from its surroundings. At the second level, the area outside the LPIS boundaries was segmented again in order to derive homogeneous objects that are representative of all of the remaining classes. In the second phase, we developed a rule-based classification model to classify image objects into the burn severity classes. In our analysis, the representative features of the object classes were determined with the associated threshold value and implemented in the classification model. We regulated the membership of objects in certain classes via classification rules with fixed threshold values. Overall, the whole process resulted in the production of the previously-mentioned burn severity map.
The 100 CBI plots, which we selected based on the object-based burn-severity map, were positioned in the field by taking into account accessibility constraints. All CBI plots were located at least 100 m away from each other. Unfortunately, the Greek forest services, in an effort to aid the recovery of the ecosystem, allowed loggers to remove all of the dead trees and clear the burnt land. Due to this, we were able to use only 52 out of the initial 100 plots in our analysis.
To gauge the magnitude of ecological change caused by fire in the field, we adopted the Composite Burn Index (CBI). The CBI was measured for each plot, using an adapted version of the FIREMON landscape assessment methodology [1,8]. Burn severity was recorded using the FIREMON field form [8] and field digital photography. Five burn severity factors were visually assessed for each vegetation strata ((i) substrates; (ii) grasses, herbs, and small shrubs below 1 m; (iii) tall shrubs and trees up to 5 m; (iv) trees from 5 to 20 m; and (v) large trees taller than 20 m) within each plot, and were given a score ranging from zero to three (unburned to highest severity) [39]. In the following, a total CBI score was calculated for each plot by averaging the five rating factors measured in field. Finally, for the purposes of the statistical analysis, the CBI values were linearly scaled to the range between 0 and 1.
To provide additional information regarding the statistical attributes of the derived CBI variables (Time1CBI and Time2CBI) we plotted their histograms in Figure 3. From the figure we can see that in both examined cases (Time1 at 2011 and Time2 at 2012), the CBI scores are highly clustered around the high values (and, therefore, they are not very dispersed). This indicates that most of the plots have high CBI scores and implies that most of the areas examined in the field have been severely affected by the fire. Visual analysis of the two histograms indicates that there is a significant difference between the two distributions. In the case of the first histogram ( Figure 3a) the majority of the CBI plots have high scores that fall into the 2.75-3 class whereas, in the second histogram, the number of plots that belong to this class has decreased significantly. In the first case examined (Time1CBI), the concentration of the CBI scores around the high values can be attributed to the fact that, immediately after the fire, the effects of burn severity were very intense. This find also agrees with our field observations. In the second case (Time2CBI) ( Figure 3b) the fact that the number of CBI plots with the highest scores is decreased can be explained by the fact that one year after the fire the ecosystem has started recovering. values. Overall, the whole process resulted in the production of the previously-mentioned burn severity map. The 100 CBI plots, which we selected based on the object-based burn-severity map, were positioned in the field by taking into account accessibility constraints. All CBI plots were located at least 100 m away from each other. Unfortunately, the Greek forest services, in an effort to aid the recovery of the ecosystem, allowed loggers to remove all of the dead trees and clear the burnt land. Due to this, we were able to use only 52 out of the initial 100 plots in our analysis.
To gauge the magnitude of ecological change caused by fire in the field, we adopted the Composite Burn Index (CBI). The CBI was measured for each plot, using an adapted version of the FIREMON landscape assessment methodology [1,8]. Burn severity was recorded using the FIREMON field form [8] and field digital photography. Five burn severity factors were visually assessed for each vegetation strata ((i) substrates; (ii) grasses, herbs, and small shrubs below 1 m; (iii) tall shrubs and trees up to 5 m; (iv) trees from 5 to 20 m; and (v) large trees taller than 20 m) within each plot, and were given a score ranging from zero to three (unburned to highest severity) [39]. In the following, a total CBI score was calculated for each plot by averaging the five rating factors measured in field. Finally, for the purposes of the statistical analysis, the CBI values were linearly scaled to the range between 0 and 1.
To provide additional information regarding the statistical attributes of the derived CBI variables (Time1CBI and Time2CBI) we plotted their histograms in Figure 3. From the figure we can see that in both examined cases (Time1 at 2011 and Time2 at 2012), the CBI scores are highly clustered around the high values (and, therefore, they are not very dispersed). This indicates that most of the plots have high CBI scores and implies that most of the areas examined in the field have been severely affected by the fire. Visual analysis of the two histograms indicates that there is a significant difference between the two distributions. In the case of the first histogram ( Figure 3a) the majority of the CBI plots have high scores that fall into the 2.75-3 class whereas, in the second histogram, the number of plots that belong to this class has decreased significantly. In the first case examined (Time1CBI), the concentration of the CBI scores around the high values can be attributed to the fact that, immediately after the fire, the effects of burn severity were very intense. This find also agrees with our field observations. In the second case (Time2CBI) (Figure 3b) the fact that the number of CBI plots with the highest scores is decreased can be explained by the fact that one year after the fire the ecosystem has started recovering.

Satellite Data
Two satellite images collected by the GeoEye sensor (2 m) were used for the statistical analysis ( Table 1). The first image was captured immediately after the fire (20 November 2011) and the second image was captured one year later (23 August 2012). Henceforth, the preprocessed images will be referred, respectively, as GeoEyetime1 and GeoEyetime2.

Satellite Data
Two satellite images collected by the GeoEye sensor (2 m) were used for the statistical analysis ( Table 1). The first image was captured immediately after the fire (20 November 2011) and the second image was captured one year later (23 August 2012). Henceforth, the preprocessed images will be referred, respectively, as GeoEye time1 and GeoEye time2 .
The two GeoEye images were atmospherically corrected, using the ATCOR 2 software (Wessling, Germany), and ortho-rectified with the help of the 5 m Digital Elevation Model [40]. The total RMS errors associated with the ground control points (GCPs) did not exceed 0.5 pixels in either of the two images. Auxiliary data, such as the Land-parcel Identification system (LPIS) data for Greece and a WorldView-2 image captured soon after the fire event, were also used during the object-based classification process.
To investigate the relationship between the VHR satellite data (GeoEye) and the field-measured burn severity, for each CBI plot we identified a 60 mˆ60 m area centered on each CBI field point, where we computed five remotely-sensed variables. First, we computed the mean value of each spectral band (i.e., blue, green, red, near-infrared). Then, since various studies have shown that the variation of the vegetation cover after a severe fire event is highly correlated with the magnitude of the burn severity [12,33,41] we computed the mean Normalized Difference Vegetation Index (NDVI) [42] per sampling plot. We applied the same procedure to both post-fire images (i.e., GeoEye time1 and GeoEye time2 ). Overall, we derived 10 variables. Figure 4 presents the histograms of the derived variables. From the figure we can see that variation exists in all variables. Ultimately, all of the aforementioned variables were linearly scaled to the range (0, 1) [7]. For better understanding of the NDVI variable, Figure 5 presents the variability within the CBI plots of different levels of burn severity. The two GeoEye images were atmospherically corrected, using the ATCOR 2 software (Wessling, Germany), and ortho-rectified with the help of the 5 m Digital Elevation Model [40]. The total RMS errors associated with the ground control points (GCPs) did not exceed 0.5 pixels in either of the two images. Auxiliary data, such as the Land-parcel Identification system (LPIS) data for Greece and a WorldView-2 image captured soon after the fire event, were also used during the object-based classification process.
To investigate the relationship between the VHR satellite data (GeoEye) and the field-measured burn severity, for each CBI plot we identified a 60 m × 60 m area centered on each CBI field point, where we computed five remotely-sensed variables. First, we computed the mean value of each spectral band (i.e., blue, green, red, near-infrared). Then, since various studies have shown that the variation of the vegetation cover after a severe fire event is highly correlated with the magnitude of the burn severity [12,33,41] we computed the mean Normalized Difference Vegetation Index (NDVI) [42] per sampling plot. We applied the same procedure to both post-fire images (i.e., GeoEyetime1 and GeoEyetime2). Overall, we derived 10 variables. Figure 4 presents the histograms of the derived variables. From the figure we can see that variation exists in all variables. Ultimately, all of the aforementioned variables were linearly scaled to the range (0, 1) [7]. For better understanding of the NDVI variable, Figure 5 presents the variability within the CBI plots of different levels of burn severity.   In this study we used the NDVI index as a proxy for quantifying the magnitude of post-fire vegetation response in GeoEye images (Equation (1) NDVI is a commonly used index for post-fire vegetation response in the literature [2,33,43] and is based on the absorption and reflection characteristics of plants in the red and near-infrared spectral regions [12]. In addition to NDVI, previous work [44] has proposed the use of the Normalized Burn Ratio (NBR) [26], and the differenced Normalized Burn Ratio (dNBR) [26] for detecting and monitoring post-fire vegetation changes over time. Unfortunately, their use was not possible in this study because they require data from regions of the optical spectrum that are not available in GeoEye images (Table 1). More specifically, the NBR index requires multispectral sensors with a NIR band between 0.76-0.9 µm and a SWIR band between 2.08-2.35 µm. In this study we used the NDVI index as a proxy for quantifying the magnitude of post-fire vegetation response in GeoEye images (Equation (1)): NDVI is a commonly used index for post-fire vegetation response in the literature [2,33,43] and is based on the absorption and reflection characteristics of plants in the red and near-infrared spectral regions [12]. In addition to NDVI, previous work [44] has proposed the use of the Normalized Burn Ratio (NBR) [26], and the differenced Normalized Burn Ratio (dNBR) [26] for detecting and monitoring post-fire vegetation changes over time. Unfortunately, their use was not possible in this study because they require data from regions of the optical spectrum that are not available in GeoEye images (Table 1). More specifically, the NBR index requires multispectral sensors with a NIR band between 0.76-0.9 µm and a SWIR band between 2.08-2.35 µm.

Topographic Data
We used a five-meter (5 m) DEM [40] to derive three topographic parameters: (a) aspect; (b) slope; and (c) elevation. The produced grids of aspect and slope were calculated in radians, while the elevation in meters. Aspect is a variable measured in circular scale, with values ranging from 0˝to 360˝degrees. The fact that angles 0˝and 360˝are identical in this scale imposes the conversion of aspect values into sine and cosine values [45]. As a result, the aspect values were initially converted from degrees to radians and then to sine and cosine values using ArcMap Raster Calculator. The converted aspect sine values ranged from´1 (at due west) to 1 (at due east) [46], while cosine values ranged from´1 (at due south) to 1 (at due north) [45].
The values of all the topographic parameters (aspect cosine, aspect sine, slope, and elevation) were obtained for all CBI sampling points by computing the mean of each variable within a window of 12 pixelsˆ12 pixels (60 mˆ60 m plot area). Figure 6 presents the histograms of the derived variables. From the figure we can see that variation exists in all four variables. This implies that the CBI plots were collected from areas with diverse topographic characteristics. Finally, the values of these variables were linearly scaled between 0 and 1 in order to facilitate the statistical analysis.

Topographic Data
We used a five-meter (5 m) DEM [40] to derive three topographic parameters: (a) aspect; (b) slope; and (c) elevation. The produced grids of aspect and slope were calculated in radians, while the elevation in meters. Aspect is a variable measured in circular scale, with values ranging from 0° to 360° degrees. The fact that angles 0° and 360° are identical in this scale imposes the conversion of aspect values into sine and cosine values [45]. As a result, the aspect values were initially converted from degrees to radians and then to sine and cosine values using ArcMap Raster Calculator. The converted aspect sine values ranged from −1 (at due west) to 1 (at due east) [46], while cosine values ranged from −1 (at due south) to 1 (at due north) [45].
The values of all the topographic parameters (aspect cosine, aspect sine, slope, and elevation) were obtained for all CBI sampling points by computing the mean of each variable within a window of 12 pixels × 12 pixels (60 m × 60 m plot area). Figure 6 presents the histograms of the derived variables. From the figure we can see that variation exists in all four variables. This implies that the CBI plots were collected from areas with diverse topographic characteristics. Finally, the values of these variables were linearly scaled between 0 and 1 in order to facilitate the statistical analysis.

Methodology
To quantify the potential of VHR GeoEye satellite data to express the variations in the field-measured burn severity we used redundancy analysis (RDA) [47]. The RDA method was first developed by Rao [48,49]. RDA is a form of asymmetric canonical ordination technique widely used by ecologist and palaeoecologists [50]. The term 'asymmetric' means that the two data matrices used in the analysis do not play the same role: there is a matrix of response variables, denoted Y, which often contains community composition data, and a matrix of explanatory variables (also referred to as predictors) (e.g., environmental), denoted X, which is used to explain the variation in Y, as in regression analysis [49]. In this case, the term 'redundancy' has a specific meaning and is used

Methodology
To quantify the potential of VHR GeoEye satellite data to express the variations in the field-measured burn severity we used redundancy analysis (RDA) [47]. The RDA method was first developed by Rao [48,49]. RDA is a form of asymmetric canonical ordination technique widely used by ecologist and palaeoecologists [50]. The term "asymmetric" means that the two data matrices used in the analysis do not play the same role: there is a matrix of response variables, denoted Y, which often contains community composition data, and a matrix of explanatory variables (also referred to as predictors) (e.g., environmental), denoted X, which is used to explain the variation in Y, as in regression analysis [49]. In this case, the term "redundancy" has a specific meaning and is used interchangeably with explained variance [49,51]. RDA can also be considered as an extension of principal component analysis because the ordination axes are linear combinations of the response variables (Y), but they are also constrained to be a linear function of the predictor [49]. Accordingly, the ordination of the dependent variables is forced to be explained by the configuration of the independent variables such that the RDA axes represent the percentage of the variance of the response variables explained by the predictors [52,53]. In a RDA biplot, the correlation among variables is given by the cosine of the angle between two vectors. Vectors pointing in roughly the same direction indicate a high positive correlation, vectors crossing at right angles correspond to a near zero correlation, and vectors pointing in opposite directions show a high negative correlation [54]. Environmental variables with long vectors are the most important in the analysis [55].
In order to identify the explanatory power of field-measured burn severity and topography on the vegetation post-fire response, we carried out the RDA using the topographic and CBI data as predictor variables (X), while the remotely-sensed VHR variables represented the response variables (Y). The rationale was to see if, through time, the satellite-based data can somehow express the field-measured burn severity and to quantify to what extent, in our study case, burn severity rather than topography influences the satellite-based data. RDA was performed using XLSTAT statistical software. We tested the significance of the analysis by running 500 permutations.

Results
The results obtained by applying the RDA method are presented in Figure 7 and Tables 2 and 3. Figure 6 illustrates that the first, second. and third axis represented, respectively, 95.83%, 2.81%, and 1.10% of the total variation of the multivariate system (p-value < 0.001); the first axis is, hence, the component which accounts for nearly the totality of the system variation. interchangeably with explained variance [49,51]. RDA can also be considered as an extension of principal component analysis because the ordination axes are linear combinations of the response variables (Y), but they are also constrained to be a linear function of the predictor [49]. Accordingly, the ordination of the dependent variables is forced to be explained by the configuration of the independent variables such that the RDA axes represent the percentage of the variance of the response variables explained by the predictors [52,53]. In a RDA biplot, the correlation among variables is given by the cosine of the angle between two vectors. Vectors pointing in roughly the same direction indicate a high positive correlation, vectors crossing at right angles correspond to a near zero correlation, and vectors pointing in opposite directions show a high negative correlation [54]. Environmental variables with long vectors are the most important in the analysis [55].
In order to identify the explanatory power of field-measured burn severity and topography on the vegetation post-fire response, we carried out the RDA using the topographic and CBI data as predictor variables (X), while the remotely-sensed VHR variables represented the response variables (Y). The rationale was to see if, through time, the satellite-based data can somehow express the field-measured burn severity and to quantify to what extent, in our study case, burn severity rather than topography influences the satellite-based data. RDA was performed using XLSTAT statistical software. We tested the significance of the analysis by running 500 permutations.

Results
The results obtained by applying the RDA method are presented in Figure 7 and Tables 2 and 3. Figure 6 illustrates that the first, second. and third axis represented, respectively, 95.83%, 2.81%, and 1.10% of the total variation of the multivariate system (p-value < 0.001); the first axis is, hence, the component which accounts for nearly the totality of the system variation. The RDA ordination plot was performed using the topographic (slope, aspect ,elevation) and field-measured CBI data as explanatory variables, and the remotely-sensed variables as response variables. The percent of the total variation that is explained by the first two canonical axes (F1 and F2) are 95.83% and 2.81%, respectively.
According to Table 2, the most contributing variables to the first axis are the CBI ones at both times of measurement (Time1CBI and Time2CBI), while the topographic data are the least Figure 7. RDA biplot showing the explanatory variables (X) with lines (different color for each variable) and the response variables with orange dots.The RDA ordination plot was performed using the topographic (slope, aspect ,elevation) and field-measured CBI data as explanatory variables, and the remotely-sensed variables as response variables. The percent of the total variation that is explained by the first two canonical axes (F1 and F2) are 95.83% and 2.81%, respectively. According to Table 2, the most contributing variables to the first axis are the CBI ones at both times of measurement (Time1CBI and Time2CBI), while the topographic data are the least contributing variables. The response variables were mostly explained by the predictors NDVITime1 and NDVITime2 (Table 3); to the contrary the raw spectral data are only partially related to the considered explanatory variables.

Discussion
Burn severity is commonly defined as a measure of the magnitude of the short-to long-term ecological changes caused by fire [9]. Due to the challenges faced by resource managers in maintaining post-fire ecosystem health, there is a need for methods to assess the ecological consequences of disturbances and to monitor the vegetation changes occurring after a burning event [7].
Generally, topography plays an important role in vegetation dynamics prior and after the fire. Prior to the fire, topography influences the microclimate (temperature, precipitation, and direct solar radiation), plant productivity, and biomass accumulation that directly affects the amount of biomass available to burn during the outbreak of fire [18]. After the fire event, the topographic variations in slope, aspect, and elevation influences the microclimate at a local scale which, in turn, affects the post-fire regeneration rates [47].
Our statistical analysis results indicate that the two NDVI time variables are slightly correlated with slope, aspect, and elevation. This important finding can be attributed to two possible reasons. First, since the vegetation response to fire depends on many factors, it is likely that various factors other than the topography have an important influence on the post-fire vegetation response. Second, prior to the fire the studied area was mainly covered by Pinus reforestations, all of which were established on areas with similar altitudes. For this reason a large number of samples (CBI field measurements) were selected from areas with similar altitude. The fact that the statistical analysis is based on a dataset with samples taken from similar altitudes may explain the lack of elevation influence on the post-fire vegetation response. Even though these statements need further investigation in order to arrive to concrete conclusions, it is clear that the restoration planning can greatly benefit from insights into the vegetation response to local environmental conditions. Following, the results showed that there is a high negative correlation between the two CBI time variables and the corresponding NDVI (NDVITime1 and NDVITime2). This negative correlation is attributed to the fact that during the course of time, the ecosystem responds to disturbance which, in turn, affects the degree of burn severity. That is, the increased satellite-derived NDVI values are associated with lower degree of burn severity in the field. Several studies have shown that the variation of the vegetation cover after a severe fire event is highly correlated with the magnitude of burn severity [7,12,33]. This indicates that the NDVI data extracted from the GeoEye imagery, was able to capture the ecosystem response to fire, which was quantified using the field-based CBI measurements.
The strong, negative correlation between the CBI and NDVI variables demonstrated that GeoEye images captured soon after the fire are capable of detecting the short term (up to one year) bi-temporal variations in burn severity. A possible explanation is that the fire effects are quite homogeneous and easy to measure short after the disturbance [18]. Even though in this study we focused on the short-term relationship between the GeoEye images and the fire effects, it would be beneficial to extend our analysis by looking at the long-term relationship. To do so, we will need a much larger dataset of field-measured burn severity that is not currently available. Extracting information from imagery with higher spectral, spatial, and temporal resolutions is expected to augment existing information regarding ecosystems response after fire.
In summary, the results of this study showed that satellite-derived NDVI time variables are highly associated with the CBI field-measured variables. This observation emphasizes the role of the VHR GeoEye remotely-sensed data in monitoring burn severity through space and time. This work is a preliminary study, in which the relationship between the field-measured CBI index and the VHR GeoEye satellite data are explored. To the best of our knowledge there is no burn severity index found in literature to be specifically designed for VHR satellite data. Hence, the results of this study could possibly provide a first guide for designing a new burn severity index specifically focused on VHR satellite data. Since the results of this study have the limit to be site specific, further research is needed to fully understand the interrelationships between topography, field measured burn severity, and the post-fire vegetation response, which, in this specific case, has been quantified using VHR GeoEye satellite data. In this context, it would be very interesting to examine how the previously-mentioned relationships vary over longer time spans and different study areas.

Conclusions
Fire and burn severity studies will continue to depend on satellite data and field-based burn severity data. As generally known, field sampling is very complicated due to several factors, like survey costs, the accessibility of the locations, the limited surface coverage, etc. Finding a direct correlation between in situ measurements and remotely-sensed data also in this research field may have useful implications in terms of post-fire analysis, recovery strategies, intervention and land management work. Therefore, it is to the best interest of stakeholders and local managers to improve the quality of the burn severity information. To do so, it is critical to associate the burn severity field measurements with the appropriate satellite image. The results of this study indicate, that short term burn severity can be successfully estimated using GeoEye images. Accurate burn severity estimates will augment existing information regarding the ecosystems response after fire. Furthermore, this kind of study can open new perspectives in terms of environmental applications of new sensor products.