Combining aerial LiDAR and multi-spectral imagery to assess post-fire 1 regeneration types in a Mediterranean forest

Abstract


Introduction
Wildfires have been the most important natural disturbance in Mediterranean ecosystems since at least the late Quaternary (Carrión et al. 2010;Pausas et al. 2008).In the Mediterranean region, humans have lived with fire and used it in their agricultural and rural activities during millennia.However, in the last mid century an increase of ignition sources has taken place, causing increasing fire risk and frequency of uncontrolled fires (Gonzalez-Olabarria and Pukkala 2007; San- Miguel-Ayanz et al. 2013;Schelhaas et al. 2003).This increment has been attributed to a combination of factors related, among others, to land abandonment, the increase in the number of days with extreme-fire-hazard weather, and the increasing number of humanrelated ignition sources (Loepfe et al. 2010).In the Iberian Peninsula, large wildfires (>500 ha) have hit almost every type of forest ecosystem over the last decades, representing almost half of the total burnt area during this period (Cubo María et al. 2012).This scenario of increased fire impacts may be further magnified in the future, as climate forecasts point to prolongation of droughts and hot spells likely to further aggravate forest fire risk (Keenan et al. 2011;Lindner et al. 2010;Piñol et al. 1998;Resco de Dios et al. 2006).
In addition to changes in fire frequency and extent, some areas have suffered an increase in the occurrence of high-intensity crown fires affecting forest types that had not historically been subject to them, such as the montane (sub-Mediterranean) forests dominated by Mediterranean black pine (Pinus nigra Arn ssp.salzmannii) or Scots pine (Pinus sylvestris L.) (Ordoñez and Retana 2004;Pausas et al. 2008;Vilà-Cabrera et al. 2012).These pine species lack direct post-fire regeneration mechanisms and usually show almost nil regeneration after crown fires (Pausas et al. 2008).In these areas substitution of the pines by resprouting hardwoods (mostly Quercus species) are generally observed when the latter were present in the understory of the burnt stands (Puerta-Piñero et al. 2011;Rodrigo et al. 2004).In the particular case of P. nigra forests, the thick bark and high self-pruning ability of this species allow the persistence of some surviving trees in the form of small islands interspersed across the burned landscape (Ordoñez et al., 2005).The existence of such islands, together with the moderate shade-tolerant character of this species (Niinemets and Valladares, 2006), may lead to greater opportunities for its mid-to long-term colonization of the burned area.In contrast, in the areas without presence of sprouting species or remaining trees, important problems of soil recovery by woody vegetation are likely to arise, leading to increasing soil erosion and forest degradation, and causing long-term environmental and economic damage (Selkimäki et al. 2012).
Several factors inherent to forests in the Mediterranean basin, such as the small size of forest ownership, slow growth rate due to limited water availability, or the dominance of mountainous terrain, raise the cost of forestry operations and hamper the development of management actions (EFI 2010;Gonzalez-Olabarria et al. 2008), including post-fire restoration measures (Espelta et al. 2003a).These impediments may be in part tempered if detailed and spatially continuous assessments directed to classify the burned area by regeneration success, vegetation types and/or management are available (Vallejo et al. 2012).However, the extent of these types of disturbances makes these detailed and spatially continuous assessments often unattainable through field data gathering.Thus, the general approach used to run post-fire assessments usually involves implementing a statistical sampling design to decide the location of inventory plots (Pausas et al. 2004;Proença et al. 2010;Puerta-Piñero et al. 2012;Shatford et al. 2007) and then applying modeling (i.e.extrapolation) techniques.The use of remote sensing (RS) data has recently emerged as an efficient alternative to provide adequate regeneration assessments over large areas affected by forest disturbances.In this line, some recent studies have used indicators like the Normalized Difference Vegetation Index (NDVI) (Tucker 1979) extracted from satellite imagery to evaluate short-term regeneration success and recovery rates in terms of the total plant cover (Belda and Meliá 2000;Diaz-Delgado et al. 2003;Gouveia et al. 2010;Leeuwen et al. 2010;Viedma et al. 1997;Vila and Barbosa 2010) and even distinguishing the cover of woody shrubs or tree-species regeneration (Riaño et al. 2002;Vicente-Serrano et al. 2011).
From an operational point of view, the value of these post-disturbance regeneration assessments might significantly increase if complemented with data describing the structure of the vegetation (e.g.relative cover of the different strata or functional groups, etc.).In relation to this, Light Detection And Ranging (LiDAR) has recently emerged as a powerful tool for characterizing such structural attributes (Lefsky et al. 2002;Wulder et al. 2012) and it holds great potential for evaluating medium-to long-term post-disturbance regeneration (Debouk et al. 2013;Johnstone et al. 2004).Combining LiDAR data with multi-spectral images could further improve the assessment of post-fire vegetation by making it possible to accurately characterize low vegetation attributes (Erdody and Moskal 2010;Riaño et al. 2007) and to identify individual tree species (Hill et al. 2010;Holmgren et al. 2008;Waser et al. 2011).It also enables to identify vegetation typologies (Bork and Su 2007;Goetz et al. 2010;Mutlu et al. 2008), which may be especially appropriate in the case of low resolution data.Typological characterization of forest stands for management purposes have been widely used on adult forest stands (e.g.Aubury et al. 1990;Herbert and Rebeirot 1985;Martin-Alcon et al. 2012;Reque and Bravo 2008), providing detailed and objective classifications of the stands according to their structural and floristic attributes.In the case of regenerating forest stands, such attributes are expected to largely influence the long-term dynamics of the forest and its accompanying ecosystem services (Elmqvist et al. 2003).
The aim of this study was to develop and test a simple and cost-effective methodology for conducting a detailed and spatially continuous characterization of post-fire regeneration (i.e.recovery of vegetation cover) from widely available remote sensing data, using a large forest fire in the Mediterranean basin as study case.To this purpose, we defined and mapped the main post-fire regeneration types in the study area using a combination of low-resolution airborne LiDAR data (useful for characterizing the three-dimensional structure of plant canopies) with NDVI data computed from multi-spectral data from aerial images, potentially suitable for differentiating between species groups with different spectral signatures.We applied this methodology to the particular case of a large wildfire in a Mediterranean zone, but we anticipate that it should be easily adapted to assess vegetation responses to other forest disturbances and a wide range of forest ecosystems, such as large-scale disease outbreaks or windthrow events.

Material and methods
The methodological approach adopted to assess post-fire regeneration types based on remote sensing data can be summarized in the following steps: (1) remote sensing data processing (raw LiDAR data and raw imagery data), (2) establishment and measurement of a set of field inventory plots located in two subsets of the study area (training areas), (3) selection of the most informative remote sensing (RS) variables to be used as candidate predictors, (4) definition of alternative post-fire regeneration typologies based on training field data, (5) selection of a regeneration typology, on the basis of the classification accuracy of a supervised classification model using RS variables, (6) extrapolation of the regeneration typologies to the whole study area through the application of the corresponding supervised classification model, and (7) validation of the resulting classification using a set of observation points randomly located across the study area (validation stands) (Figure 1).

Study area
The study was conducted in an area affected by a large forest fire in the central region of Catalonia (NE Spain).This wildfire burned nearly 24,000 ha in 1998 (Figure 2), leading to almost complete loss of forest cover (Rodrigo et al. 2004).Prior to the wildfire, the area was composed of the typical Mediterranean mosaic-like landscape, with cultivated lands and scrublands (7,700 and 2,000 ha, respectively) interspersed with forest areas (14,000 ha).
According to data in the Forest Ecological Inventory of Catalonia (Burriel et al. 2004), the dominant species before the wildfire was black pine (Pinus nigra Arn.ssp salzmanii), covering 75% of the total forest surface.Aleppo pine (Pinus halepensis Mill.) was the second most abundant species, covering 15% of the forest surface.Both black and Aleppo pine forests appeared as pure stands or as two-layered stands, the latter with pine dominating the overstory and resprouting hardwoods the understory (Figure 2).The remaining forest surface was mainly covered by mixed hardwoods (mainly Quercus ilex L. and Quercus cerrioides Wilk) and Pinus sylvestris L. stands.The area presents a gentle relief with low hills ranging in elevation from 480-910 m.a.s.l. and a dry-sub-humid to sub-humid Mediterranean climate according to the Thornthwaite index.

#Figure 2#
From the total forest area affected by the 1998 wildfire, we selected a portion of 3,973 ha that burned at high severity (i.e. with no remaining adult trees) (Figure 2).Areas affected by roads, agricultural fields and other non-forest uses were excluded from the analysis, leaving a total of 3,192 ha to be used for the post-fire regeneration assessment.Additional criteria to determine the study area were the characteristics of the available remote sensing data, as a minimum seasonal similarity in both LiDAR data and aerial images was required in order to avoid potential mismatches due to the use of data from different seasons (e.g.before and after budbreak in broadleaves).In our particular case, the LiDAR data were taken in the early summer of 2009 and the aerial photographs were taken in the early spring of 2011.

LiDAR data acquisition and processing
LiDAR data were acquired via the LidarCAT project (Cartographic and Geological Institute of Catalonia).The entire study area was covered by seven flight lines captured with an ALS50 II LiDAR sensor mounted on a Cessna Caravan 208B aircraft.The flight date ranged from May 30, 2009 to June 3, 2009.The average first-return point density of these LiDAR data was 0.5 pulses•m -2 and each pulse captured up to 4 returns.LiDAR point coordinates were adjusted according to the methodology proposed by Kornus and Ruiz (2003).After filtering the clouds, low intensity returns and air points, LiDAR returns were automatically classified with the TerraScan ground classification routine (Terrasolid 2012) in order to differentiate ground from non-ground returns.Ground routine classified ground points by iteratively building a triangulated surface model.The routine started by selecting some local low points that are confident hits on the ground.It was assumed that any 40 by 40 meter area would have at least one hit on the ground and that the lowest point would be a ground hit.The routine built an initial model and started molding the model upwards by iteratively adding new laser points to it.
Each added point made the model following the ground surface more closely.Iteration parameters of angle (4º) and distance (1 m) determined how close a point must be to a triangle plane for being accepted as a ground point and added to the model.The remaining non-ground points were visually inspected and manually classified in order to extract wires and towers (since the low point density of the LiDAR data advised against the implementation of automatic algorithms for this purpose).Building roofs were extracted automatically using a new routine which classified non-ground points which formed a planar surface of at least 40 m 2 .Finally, the remaining non-ground returns were classified as vegetation, and LiDAR height was replaced by its vertical distance to a triangular irregular network (TIN) generated from ground returns.
The LiDAR vegetation point cloud was analyzed in a 10x10 meter regular grid with the FUSION software system (McGaughey and Carson 2003) to obtain structured statistical information on the laser returns.A height of 0.15 meters, predefined from the vegetation field inventory, was used as a threshold to separate woody and shrubby vegetation from herbaceous vegetation and bare soil.Intensity data was normalized by the range normalization to a userdefined standard range (Donoghue et al. 2007;García et al. 2010): where I' is the normalized intensity, I is the raw intensity value, R is the range (LiDAR sensortarget distance), and R S is the standard range (in this case 2,225 m, which corresponds to the mean range in the study area).This method eliminated the effect of path length variations on the intensity recorded by the system, providing values equivalent to the intensity that would have been recorded if all points were at the same range.The exact value of the range was not available for each LiDAR point, therefore it was approximated by the use of the difference between the average altitude of the flight and the elevation of each point, and the projected distance between the point and the flight line.This approach should cause smaller errors than the option of using the altitude difference instead of the actual range (Coren and Sterzai 2006;Kukko et al. 2008).Only the first echoes returned from the aboveground vegetation were used for computing a set of LiDAR height and intensity metrics suitable to characterize forest vegetation structure and floristic composition (e.g.Donoghue et al. 2007;García et al. 2010;Holmgren et al. 2008;Latifi et al. 2012;Morsdorf et al. 2006).These metrics were the following: (i) percentiles of pulse height (H_P05, H_P10, H_P20, ..., H_P80, H_P90, H_P95) and pulse intensity (Int_P05, Int_P10, Int_P20, ..., Int_P80, Int_P90, Int_P95); (ii) the mean (H_MEAN, Int_MEAN), mode (H_MODE, Int_MODE), standard deviation (H_SD, Int_SD), coefficient of variation (H_CV, Int_CV), interquartile range (H_IQ, Int_IQ), skewness (H_Sk, Int_Sk) and kurtosis (H_Kur, Int_Kur) of both pulse height and pulse intensity values.In addition, the following three ratios were generated: percentage of vegetation points (height > 0.15 m) in relation to all the first returns (FR_VEG), percentage of vegetation points in relation to all returns (AR_VEG), and percentage of the first returns above the mean height (FR_abMEAN) and the mode (FR_abMODE).

Multi-spectral image acquisition and processing
NDVI imagery were obtained from the Cartographic and Geological Institute of Catalonia annual coverage flights at 22 cm GSD (Ground Sample Distance), and were generated from aerial photos with RGB and near infrared bands.The aerial photos were taken between April 2 and April 10 2011, with DMC (Digital Mapping Camera) cameras DMC-26 and DMC-14, and were later ortho-rectified without stitching.The implemented procedure to produce NDVI from DMC imagery used the original DMC LR4 files (low resolution multispectral bands) with absolute radiometric calibration and was based on manufacturer's calibration of the camera, following the methodology proposed by Martínez et al. (2012).Radiance values from the red and near infrared bands, and its respective reflectances (i.e.ratio between incoming energy from the sun and reflected energy modulated by some geometric factors such as location, data and time of exposure during acquisition) were obtained for each pixel using the following equations: where L R and L NIR are radiance values from the red and the near infrared bands, µ is a geometric factor and E 0 is the extraterrestrial solar radiance.Then, NDVI was calculated using the following equation: where R and NIR correspond to the reflectance values of each pixel from the red and the near infrared band, respectively.NDVI values were then aggregated to match the same 10×10 meter grid as that use to aggregate the LiDAR data and a set of NDVI statistics potentially suitable to describe regeneration structure and composition was calculated.The quartiles of NDVI distribution (NDVI_Q1, NDVI_Q2 and NDVI_Q3), the interquartile range (NDVI_IQ), mean (NDVI_MEAN), coefficient of variation (NDVI_CV), and standard deviation of the NDVI values (NDVI_SD) were obtained.NDVI values in the regenerating area generally ranged between 0 and 0.6.On this basis, the portion of each 10x10 meter cell presenting NDVI values within a variety of sub-ranges were computed, using 0.05 increments (0.15-0.2, …, 0.55-0.6),0.10 increments (0.15-0.25, ..., 0.45-0.55),and 0.15 increments (0.15-0.3, 0.3-0.45 and 0.45-0.6), in order to maximize the chance of differentiating species or groups of them based on their different NDVI values.

Field data gathering
A field inventory was performed in order to define post-fire regeneration types.The inventory consisted of a set of forty-four 10×10 meter plots placed along two sub-areas considered representative of the overall study area (named training areas) (Figure 2).Stratified random sampling was applied to capture the range of regeneration patterns in the training areas.For this purpose, LiDAR variables which were expected to describe plant height distribution and the distance to unburned patches were used to ensure that we covered both the structural (relative abundance of trees) and compositional gradients (since pines are expected to be more abundant near the unburned patches) characterizing the post-fire vegetation.For each plot, a set of variables describing the floristic composition and aboveground biomass structure were measured.These were the percentage of the plot area covered by soil or herbs (%SOIL), low shrub species (0.15-0.5 m tall) (%LOW SHRUBS), high shrub species (>0.5 m) (%HIGH SHRUBS), pine regeneration (%PINES) and tree species regeneration (i.e.referred to all trees including pines) (%TREES).Low and high shrub covers were later summed as %SHRUBS.We also measured mean height (cm) of all tree species (HTREES), pines (HPINES) and hardwoods (HHW).To improve the accuracy of the measured vegetation variables each 10×10 meter plot was subdivided into 2.5×2.5 m subplots.The variables were measured at the subplot level and converted into plot level values by averaging the values of the measured sub-plots within each plot.

Remote sensing data reduction
A first pre-selection of the sixty-four initially considered LiDAR-and NDVI-derived metrics was executed in order to avoid using redundant variables as candidate predictors.Principal Component Analysis (PCA) with a Varimax rotation was executed on the whole set of RS data in order to identify groups of highly collinear variables.In parallel, we computed Pearson correlation coefficients to measure the associations between RS variables and the list of forest variables measured in the field plots.A total of fourteen variables from the original sixty-four were finally preselected.These variables were the ones showing the highest correlation with field variables among each group of highly collinear variables (i.e. variables with very similar PCA factor loadings).They comprised five height metrics (coefficient of variation, kurtosis, 10th and 95th percentile height, percentage of first returns of vegetation, and percentage of first returns above the mode of the height distribution), three LiDAR intensity metrics (standard deviation, 50th and 90th percentile intensity), and five NDVI-derived variables (1st quartile NDVI, percentage of plot with NDVI between 0.45-0.5,0.25-0.35,0.15-0.3,and 0.45-0.6).

Definition of post-fire regeneration types
A subset of the variables measured in the field (%SOIL, %PINES and %TREES) were deemed sufficient to define post-fire regeneration types, as they described the distribution of the plot cover among the main strata (soil/herbs, trees and, indirectly, shrubs), and the main floristic groups of the tree layer (pines and hardwoods).The forty-four training plots were clustered according to their squared Euclidean distances using Ward's hierarchical method (Ward 1963), based on these variables.Several alternative partitions of the resulting dendrogram (i.e., those exhibiting high between-groups distance (Hair et al. 2009) and strong ecological rationale) were saved.

Development of the classification model
We used pre-selected RS variables from the field plots and the Random Forest classification algorithm (Breiman 2001) to evaluate the different partition alternatives and to obtain a model for the classification of the whole study area.Random Forest (RF) is a nonparametric supervised classification technique that has shown good performance in classifying remotely sensed data (e.g.Falkowski et al. 2009;Hudak et al. 2012).The RF technique uses a bootstrap approach for achieving higher accuracies while simultaneously addressing over-fitting problems associated with traditional classification tree models.A large number of classification trees are produced from a random subset of training data (approximately 63% random subset), permutations are introduced at each node, and the most common classification result is selected.
We ran each RF model with 5,000 bootstrap replicates (i.e., individual classification trees).
With the aim of avoiding bias in the prediction caused by imbalanced classes, the number plots per class in bootstrap samples was equal to the number of plots of the less frequent class (Evans and Cushman 2009).Out-of-bag (OOB) error estimates were calculated for each tree by classifying the portion of training data not selected in the bootstrap sample, and overall accuracy was calculated by averaging error rates across all trees in the model; this is analogous to crossvalidated accuracy estimates (Cutler et al. 2007).
A model selection procedure (i.e., variable reduction) was employed to select the optimal RS variables to use in the classification of post-fire regeneration types.The procedure was formulated to develop the most parsimonious classification model, while retaining the highest possible classification accuracy.We ran a Random Forest model selection function that uses Model Improvement Ratio (MIR) standardized importance values (Evans et al. 2011;Evans and Cushman 2009) to objectively choose the most important RS variables for predicting the regeneration type of each plot.The MIR uses the permuted variable importance, represented by the mean decrease in OOB error, standardized from zero to one.The variables were subset using 0.10 threshold increments on the original model's variable importance, with all variables above the threshold retained for each model.Each subset model was compared and the model that exhibited the lowest total OOB and lowest maximum within-class error was selected.After selecting the best model for each one of the different alternative partitions, we determined the final regeneration typology as that whose classification model had obtained the highest classification accuracy.For this study, the RF algorithm was implemented using the RandomForest package (Liaw and Wiener 2002) in the R statistical program (R Development Core Team 2007).

Mapping of regeneration types and field-based validation
Each 10 x 10 m plot in the study area was classified into one of the regeneration types by applying the final classification model on the RS variables.Regeneration types were then mapped and a set of landscape metrics were generated for explaining the spatial configuration of the landscape according to the distribution of regeneration patches.These size-and shaperelated metrics were the area and percentage covered by each regeneration type, the percentage of the type area by patch size classes, the mean patch size (PSm) and its Coefficient of Variation (PScv), and the Mean Shape Index (MSI), as defined by Rempel et al. (2012).
The validity of the final classification for the whole study area was evaluated by comparing the regeneration type assignations made from supervised classification of remotesensing data (using the classification model generated by RF) against visual assessments obtained in the field.For this purpose, 15 sampling points per regeneration type were distributed along the study area using stratified random sampling (Figure 2).An independent user was asked to assign each of these stands to one of the regeneration types, based on a visual examination of the stand characteristics, using the type descriptions in terms of forest variables (see Table 1 and Figure 3) as a guide.Measures of height and cover of the different vegetation strata were taken when visual estimation was not obvious.Finally, a confusion matrix was generated to estimate the accuracy of the remote sensing-based assignations with respect to the visual ones.

Post-fire regeneration typologies
Dendrogram resulting from the application of Ward's clustering method on the variables measured in the 44 field plots 1 led to three possible choices for the placement of the cut-off point of the hierarchical tree (i.e. three partition alternatives), creating a 4-type, a 5-type and a 6type partitions.By comparing the classification accuracy of the models developed using Random Forest for the 4, 5, and 6 type partitions, we found the 5-type one to be the best in terms of overall accuracy, with 79.55% of the plots being correctly classified (i.e.OBB estimate of error rate was 20,45%), with the use of 6 RS variables as predictors.The 4-type 2 and 6 -type 3 partitions presented accuracies of 75% and 71,45% respectively, and none of them selected less than 6 RS variables as predictors.Therefore, type-5 partition was selected as the basis for representing the regeneration types (Table 1).The characteristics of the five regeneration types were as follows: the first type was characterized by a very low vegetation cover (TYPE 1); a second type had very low tree regeneration but remarkable shrubs abundance (TYPE 2); the third type included hardwoods in low to moderate cover, mixed with shrubs (TYPE 3); the fourth type presented high cover of hardwoods regeneration (TYPE 4); and the fifth type  of stands dominated by pine regeneration in moderate to high cover (TYPE 5) (Table 1, Figure 3).

#Table 1# #Figure 3#
The six RS variables selected to identify the five regeneration types were: the 10th and 95th percentile height (H_P10 and H_P95), the percentage of first returns above the mode of the height distribution (H_abMODE), the coefficient of variation of the height distribution (H_CV), the 1st quartile NDVI (NDVI_Q1) and the percentage of plot with NDVI between 0.15-0.3(NDVI_15_30) (Figure 4).OOB classification accuracies of the RF model by regeneration type were 80% for TYPE 1, 100% for TYPE 2, 81.8% for TYPE 3, 70% for TYPE 4, and 76.9% for TYPE 5.

Mapping post-fire regeneration types and validation with visual estimates
The whole study area (3191.8ha) was classified into the five regeneration types (Figure 5).
More than a decade after the wildfire, almost a half of the area (42.6%) showed high tree cover with clear dominance of resprouting hardwoods (TYPE 4).About 17% of the study area was still dominated by resprouting hardwoods, but with low to moderate tree cover (TYPE 3).Pine regeneration appeared in moderate to high cover in approximately 11% of the study area.
Finally, the rest of the area showed sparse or no tree regeneration (TYPE 2 and 1) (Table 2).
Patches of TYPE 1 and TYPE 3 were usually very small (i.e. with an area smaller than 0.1 ha), with TYPE 1 showing the most regular patch shape (i.e. closer to a circular shape).In contrast, TYPE 4 was the one presenting the highest area covered by large patches (0.1-1 ha, 1-5 ha and larger than 5 ha) and showed the most irregular shapes.TYPE 2 and TYPE 5 were also present in a wide range of patch sizes, but in them small patches were more common than in TYPE 4.

#Figure 5# #Table 2#
The evaluation of classification accuracy compared to visual estimation resulted in 82.7% well-classified plots (between 74.1% and 91.2% well classified with a 95% confidence interval), being close to the OOB classification accuracy obtained for the training dataset.

Discussion
This study presents a methodological approach that combines airborne LiDAR data and singledate NDVI data computed from aerial images to provide precise geo-referenced and spatially continuous information on post-fire regeneration over a large area.The presented methodology was not intended for short-term assessments-for which analysis of NDVI values extracted from satellite imagery may be enough (Diaz-Delgado et al. 2003;Gouveia et al. 2010;Leeuwen et al. 2010;Vicente-Serrano et al. 2011;Vila and Barbosa 2010), but for mid-to-long-term assessments which are crucial for appropriately designing and planning post-disturbance silvicultural treatments (Alloza and Vallejo 2006;Bauhus et al. 2013;Stephens et al. 2010).The information generated by this methodology provides forest stand structure and composition information that can be used to predict its potential evolution patterns, detect of areas with persistent regeneration problems, and delimit forest stands for planning management operations.

This ultimately enables management interventions to be prioritized and allocated.
There are some differences between our RS-based assessment and evaluations as traditionally made from field-based inventories.On one hand, the level of detail in describing the post-fire vegetation communities is in general higher in field-based inventories where measurements are conducted at individual plant level (e.g.Curt et al. 2009;Proença et al. 2010).However, RSbased assessments entail less time and fewer spatial constraints than field inventories, as they are much less time-consuming and less exposed to the uncertainties associated with inference and estimation from sample surveys in areas showing fine-scale heterogeneity.Assessments based on RS data also allow estimating quantitative tree and forest attributes across large areas by using regression methods (e.g.Andersen et al. 2005;García et al. 2010;Latifi et al. 2012;Wulder et al. 2009).However, that approach requires working with high resolution LiDAR data, especially when forest attributes are related to low vegetation.In our case, we aggregated remote sensing metrics to the plot-level (10x10 meter) to make a community-based assessment instead of a plant-based one, and we estimated forest types instead of quantitative forest attributes, in concordance with the potential and limitations of the available low resolution LiDAR data.Our approach permits discerning between regeneration types by using height and intensity variables from LiDAR data and NDVI metrics computed from single-date multispectral aerial imagery, thus providing forest managers with most of the information they need for mid-term operational planning of restoration activities and silvicultural operations (Vallejo et al. 2012).In addition, the provision of continuous spatially-explicit information at the landscape scale allows easy implementation of relevant forestry applications such as the characterization and mapping of forest canopy fuels and fire risk (Erdody and Moskal 2010;García et al. 2011;González-Olabarria et al. 2012;Mutlu et al. 2008;Pierce et al. 2012;Riaño et al. 2007).The resulting maps of regeneration types also offer a starting point for the analysis of post-disturbance vegetation dynamics at large scales (Debouk et al. 2013;Goetz et al. 2010;Holmgren et al. 2008).Our methodology is appropriate for large-scale assessment, since regional and national LiDAR coverages are often collected at low-resolution.The NDVI data we used were derived from 4-band aerial imagery, which are also more and more used for national surveys.They have higher spatial resolution that the ones derived from satellite imagery, although they are more limited for multi-temporal analyses due to the generally low temporal resolution.
Our approach requires a prior categorization of the vegetation that is to be extrapolated using RS data.Such characterization can be achieved by generating vegetation typologies from the analysis of quantitative field data (as presented here), but also by using existing classifications when types of vegetation are previously known (Falkowski et al. 2009).Other areas may have different pre-fire forest composition and structure, or be more or less severely impacted by fire, and present different post-fire regeneration trajectories.Thus the number and characteristics of the regeneration types may change from one area to the other, and the RS variables selected for the classification model will be different in each particular case, since other types of vegetation may be better explained by different combination of RS variables.In our case, the best fit was attained just with a combination of LiDAR height and NDVI variables.
LiDAR normalized intensity variables were not found to be essential in discerning our regeneration types, but they could be more important in other cases (Bork and Su 2007;Donoghue et al. 2007;García et al. 2010).
Our methodology is applicable as long as the input data (LiDAR and NDVI data) is available and the chosen typologies can be accurately reproduced using supervised classification of RS data.Even so, our approach calls for a certain period of time to pass between the occurrence of the disturbance and the acquisition of the RS data to allow height differentiation between species and groups of species, in order to exploit the full potential of combining LiDAR and NDVI data for vegetation characterization.Other aspects which must be considered when applying this methodology are the characteristics of the RS data available, particularly in terms of collection dates (year and season) and spatial resolution.In this regard, increasing intensity of field sampling would be necessary for a disturbed area where date and season of collection of the RS data are not comparable.Meanwhile, the spatial resolution of the RS data available may dictate how the level of detail in characterizing regeneration types needs to be adapted.Finally, the use of this methodology also requires good knowledge of the vegetation in the specific area under study and a decision about the pretended level of detail when categorizing it, in order to adequately design the field inventory, with more heterogeneous environments typically requiring a denser network of field plots.
The five regeneration types obtained in our study area adequately represent the wide range of post-fire conditions described in previous field studies (Retana et al. 2012;Retana et al. 2002;Rodrigo et al. 2004).Furthermore, the results of the field validation revealed high (>80%) overall classification accuracy levels.For three of five types, the field validation matched the classification implemented through RS variables in over 80% of cases.The TYPEs dominated by shrubs with very low tree cover (i.e.TYPE 2) and by low to moderate tree cover dominated by hardwoods (i.e.TYPE 3) showed sensibly lower accuracy levels (of around 70-75%).Most of the misclassified TYPE 2 plots were assigned to TYPE 3, meaning that the cover of hardwoods was overestimated in the classification model based on RS data.Misclassified TYPE 3 plots were assigned to TYPE 4 (very similar to TYPE 3, as they are dominated by hardwoods regeneration, but with higher hardwoods cover).Overall, misclassification problems could be considered relatively small and occurred among closely related regeneration types.Interestingly, two of the regeneration types of major interest for management purposes (TYPEs 1, with defective cover of woody vegetation, and 5, with moderate to high cover of pine regeneration), were identified with high accuracy.In the case of TYPE 5, this may indicate very good performance of NDVI-derived variables for species differentiation, as previously shown in other studies (Hill et al. 2010;Key et al. 2001).
Regarding the spatial distribution of the forest regeneration, we found hardwoods resprouts to be presently dominating near 60% of the area formerly dominated by pines.These results demonstrate the important species dominance shifts that crown fires may induce in Mediterranean black pine forests, as reported by Rodrigo et al. (2004).However, it should be stressed that the level of detail achieved in our study did not allow proper identification of the presence of small pines on forest areas dominated by high shrubs or hardwoods, leading to certain underestimation of the pine regeneration.Interestingly, we also found less than 10% of the burned area to be in high risk of soil erosion due to a defective cover of woody vegetation.This lack of regeneration appeared in very small patches, which sometimes could be associated to the presence of small rocky outcrops.
The information generated through our methodological approach should be useful for the definition of vegetation management and restoration strategies over large areas affected by disturbances, as it helps to identify areas where woody species are struggling to regenerate (thus potentially needing restoral actions) and to provide information on areas showing regeneration success and their current vegetation structure.For example, restoration actions designed to facilitate tree cover recovery could be envisaged areas of our case study where the predominant post-fire regeneration pattern was identified as TYPE 1 and TYPE 2, especially in the larger patches (more frequently classed as TYPE 2).In parallel, early thinning interventions designed to reduce competition to increase tree growth and vigor and to reduce the amount and continuity of fuels could be planned for patches classified as TYPE 5 (Gonzalez-Olabarria et al. 2008;Moya et al. 2008;Verkaik and Espelta 2006), while coppice thinning could be proposed for stands classified as TYPE 4 (Cotillas et al. 2009;Espelta et al. 2003b;Sanchez-Humanes and Espelta 2011).

Figure 2 :
Figure 2: Location and boundaries of: a) the wildfire location in the Iberian Peninsula and Catalonia region; b) the study area and the calibration areas within the wildfire zone; c) the dominant forest types in the study area before the 1998 wildfire (CREAF 1993; DGCN 2001); d) the field plots and the validation stands within the study area.

Figure 4 :
Figure 4: Variable importance plots for the for the 5-type classification model of post-fire regeneration types.Higher value of Mean Decrease Accuracy represents higher variable importance in the model.

Figure 5 :
Figure 5: Spatial distribution of regeneration types in the study area.In black, the area covered by each of the regeneration typologies.

1
Figure 1S in Supplementary material 2 Table 1S in Supplementary material 3 Table 2S in Supplementary material