*Article* **A Framework for Multi-Dimensional Assessment of Wildfire Disturbance Severity from Remotely Sensed Ecosystem Functioning Attributes**

**Bruno Marcos 1,2,\*, João Gonçalves 1,3, Domingo Alcaraz-Segura 4,5,6, Mário Cunha 2,7 and João P. Honrado 1,2**


**Abstract:** Wildfire disturbances can cause modifications in different dimensions of ecosystem functioning, i.e., the flows of matter and energy. There is an increasing need for methods to assess such changes, as functional approaches offer advantages over those focused solely on structural or compositional attributes. In this regard, remote sensing can support indicators for estimating a wide variety of effects of fire on ecosystem functioning, beyond burn severity assessment. These indicators can be described using intra-annual metrics of quantity, seasonality, and timing, called Ecosystem Functioning Attributes (EFAs). Here, we propose a satellite-based framework to evaluate the impacts, at short to medium term (i.e., from the year of fire to the second year after), of wildfires on four dimensions of ecosystem functioning: (i) primary productivity, (ii) vegetation water content, (iii) albedo, and (iv) sensible heat. We illustrated our approach by comparing inter-annual anomalies in satellite-based EFAs in the northwest of the Iberian Peninsula, from 2000 to 2018. Random Forest models were used to assess the ability of EFAs to discriminate burned vs. unburned areas and to rank the predictive importance of EFAs. Together with effect sizes, this ranking was used to select a parsimonious set of indicators for analyzing the main effects of wildfire disturbances on ecosystem functioning, for both the whole study area (i.e., regional scale), as well as for four selected burned patches with different environmental conditions (i.e., local scale). With both high accuracies (area under the receiver operating characteristic curve (AUC) > 0.98) and effect sizes (Cohen's |d| > 0.8), we found important effects on all four dimensions, especially on primary productivity and sensible heat, with the best performance for quantity metrics. Different spatiotemporal patterns of wildfire severity across the selected burned patches for different dimensions further highlighted the importance of considering the multi-dimensional effects of wildfire disturbances on key aspects of ecosystem functioning at different timeframes, which allowed us to diagnose both abrupt and lagged effects. Finally, we discuss the applicability as well as the potential advantages of the proposed approach for more comprehensive assessments of fire severity.

**Keywords:** ecological disturbance; ecosystem functioning; EFAs; fire severity; satellite image time-series; wildfires

**Citation:** Marcos, B.; Gonçalves, J.; Alcaraz-Segura, D.; Cunha, M.; Honrado, J.P. A Framework for Multi-Dimensional Assessment of Wildfire Disturbance Severity from Remotely Sensed Ecosystem Functioning Attributes. *Remote Sens.* **2021**, *13*, 780. https://doi.org/ 10.3390/rs13040780

Received: 21 January 2021 Accepted: 18 February 2021 Published: 20 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Wildfire disturbances are considered an integral part of the natural dynamics of ecosystems in several biomes [1,2]. Notwithstanding, wildfire events can modify the composition, structure, and functioning of ecosystems and, consequently, the provision of ecosystem services to humankind [3]. Among those alterations, changes in ecosystem functioning are of particular interest, as fire can cause rapid modifications in the matter and energy budgets of ecosystems [4]. This is because attributes of ecosystem functioning have a shorter response time to disturbances than structural or compositional attributes of landscapes and are more directly connected to ecosystem services [5].

Wildfires can have profound impacts on many key aspects of the flows of matter and energy. These include dimensions of ecosystem functioning related to the biogeochemical cycles of carbon (e.g., primary productivity, biomass), and water (e.g., vegetation water content, soil moisture), as well as energy balances (e.g., albedo, latent heat, sensible heat). Namely, fires play an important role in the terrestrial biosphere carbon cycle [6], as they can cause substantial losses in carbon storage [7], including both aboveground biomass [8] as well as soil organic matter [9]. Moreover, fire-mediated changes in nutrient concentrations can ultimately limit productivity [10] as well as induce phenological changes [11]. Wildfires can also affect both water supply [12] and quality [13,14], as well as induce vegetation mortality due to prolonged loss of foliar moisture [15]. Additionally, wildfires can have negative impacts on evapotranspiration [16,17], which, in turn, also influences radiation exchanges, as evapotranspiration is related to latent heat flux [18]. Changes in wildfire frequency and/or severity may also have strong impacts on Earth's surface radiative budget [19,20]. Burn severity influences the magnitude of changes in land surface albedo [21,22], which can undergo a darkening and/or a brightening that can persist in time [23,24]. Such alterations cause an increase in land surface temperature immediately following fires [25], which influences both burned area and the duration of fires [26], since vegetation is a regulator of land surface energy fluxes [27]. Although the effects of wildfires span across multiple dimensions of ecosystem functioning, most studies evaluating those effects seldom address those multiple dimensions at the same time, thus failing to fully depict spatiotemporal changes caused by these disturbances and limiting the assessment of environmental impacts.

As wildfire events have been increasing in recent decades, both in terms of frequency and intensity [28,29], exacerbated by global climate change and by shifts in land use and forest management [30], there is an increasing need for methods to assess and monitor the ecological consequences of such disturbances [31]. Over the last decades, remote sensing (RS) has revolutionized the way environmental changes are monitored, with new satellite sensors (e.g., Sentinel-2 MSI, Landsat-8 OLI/TIRS) providing valuable data at increasingly higher temporal, spatial, and spectral resolutions. Other satellites provide long-term archives of images that are useful to support temporal studies (e.g., Terra/Aqua MODIS, Landsat 5/7 missions). This growing variety and quality of remotely sensed data enables a wide range of fire-related applications, from the detection of fire occurrences [32] to an evaluation of fire severity [33] and monitoring of post-fire recovery [31,34,35] and resilience [36]. The ability to assess and map the heterogeneous (both spatially and temporally) effects of wildfires on ecosystems with remotely sensed data is a major asset not only for risk assessment and governance but also for post-fire management and restoration [30,37–39]. In this regard, spectral indices and metrics derived from satellite multi-spectral images have been used to assess different aspects of fire severity, usually through the comparison of pre-fire versus post-fire values. For instance, spectral indices sensitive to photosynthetically active radiation and/or water content, such as the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR) have been used, in combination with data collected in-field (e.g., through the Composite Burn Index; CBI), to assess "burn severity" (e.g., [39–41]), which is a post-fire measure of organic matter loss and soil alteration [30,37,42]. Furthermore, both satellite-derived land surface temperature (LST; [43]) and albedo [22] have been used for the same end. Nevertheless, to fully characterize and

understand the effects of wildfire disturbances on multiple dimensions of ecosystems, a better translation of spectral indices into informative ecosystem variables is needed. More diverse sets of indicators are thus required to extend "burn severity" assessments to a wider variety of fire effects, i.e., towards more comprehensive approaches in "fire severity" assessment.

In recent decades, Ecosystem Functioning Attributes (EFAs) derived from satellite image time-series (SITS) have been increasingly used in a wide range of ecological applications due to their strong relation to biophysical properties and processes of ecosystems [5]. These include wildfire-related applications, such as improving the detection of wildfire disturbances in space and time [44], but also a wide range of other ecological applications, such as describing major ecological patterns [45–47], predicting and projecting species distributions [48–50], and supporting the definition of conservation priorities [51]. More specifically, EFAs extracted from SITS can provide information on inter-annual quantity (e.g., mean, minimum, maximum), seasonality (i.e., seasonal range or variability), and timing (e.g., dates of specific moments) components over multiple dimensions of ecosystem functioning. These include frequently used satellite-derived descriptors of the intra-annual dynamics of carbon gains, such as primary productivity, vegetation seasonality, and phenology (e.g., [46,52,53]). However, analogous measures can also be extracted from RS-derived proxy descriptors related to other dimensions of ecosystem functioning, such as water content, albedo, and sensible heat (e.g., [50,54]). By combining these four complementary dimensions of ecosystem functioning across each component (i.e., quantity, seasonality, and timing), remotely sensed EFAs can enable in-depth and integrative characterizations of fire severity [3,44].

The main objective of this study is to propose, describe, and showcase a framework to support enhanced assessments of fire severity, encompassing multiple dimensions of ecosystem functioning. Our approach is based on the effects of wildfire disturbances, at short (i.e., the year of the fire) to medium term (i.e., up to the second year after the fire), on descriptors of the intra-annual dynamics (i.e., EFAs) of four essential dimensions of the flows of matter and energy in ecosystems: (i) primary productivity; (ii) vegetation water content; (iii) albedo; and (iv) sensible heat. To illustrate our approach, we assessed the potential of inter-annual anomalies (i.e., deviations from the normal inter-annual variability) of a large set of EFAs extracted from MODIS data, for the period between 2000 and 2018 in the NW Iberian Peninsula, as estimators of the effects of wildfire disturbances. To this end, we compared, ranked, and analyzed the main patterns of those inter-annual anomalies at the regional scale. We then showcased the translation of inter-annual anomalies of the selected EFAs into indicators of wildfire disturbance severity for four selected burned patches within the study area and compared and analyzed their main spatial and temporal patterns at the local scale. Finally, we discussed the potential and added value of the proposed approach to improve RS-based assessment, mapping, and monitoring of wildfire disturbance severity on multiple dimensions of ecosystem functioning.

#### **2. Materials and Methods**

*2.1. Generic Framework*

#### 2.1.1. General Workflow

The proposed approach to evaluate the effects of wildfire disturbances on ecosystem functioning consists of using indicators derived from Ecosystem Functioning Attributes (EFAs) for each of the four key dimensions of ecosystem functioning (i.e., primary productivity, vegetation water content, albedo, and sensible heat). Figure 1 presents the general workflow of the proposed approach, which is composed of five steps, described in the following sub-sections.

**Figure 1.** General workflow of the proposed approach to assess wildfire disturbance severity using indicators based on Ecosystem Functioning Attributes (EFAs) derived from satellite image time-series (SITS). TCTB: Tasseled Cap Transformation Brightness feature; TCTG: Tasseled Cap Transformation Greenness feature; TCTW: Tasseled Cap Transformation Wetness feature; and LST: land surface temperature.

#### 2.1.2. Step 1—Satellite Time-Series

First, satellite image time-series (SITS) for each of the four dimensions of ecosystem functioning have to be collected from available satellite products and preprocessed (Figure 1—"Step 1"). From the large number of spectral band combinations and vegetation indices available that have been successfully used for fire applications (e.g., [43,44,55–57]), we propose, for that end, the use of land surface temperature (LST) for the "heat" dimension, as well as the Tasseled Cap Transformation (TCT) features of "Brightness" (TCTB), "Greenness" (TCTG), and "Wetness" (TCTW) for the "albedo", "primary productivity", and "vegetation water content" dimensions, respectively. LST is a measure of the thermal energy emitted by different surfaces, calibrated from thermal emissivity measured by the satellite instrument [58]. On the other hand, TCT features are well-known, sensorspecific, linear combinations of bands in the visible, near-infrared, and short-wave infrared regions of the electromagnetic spectrum [59]. These three TCT features of "Brightness", "Greenness", and "Wetness" result from a rotation of principal component axes, derived from a global sample, to be aligned with the biophysical parameters of albedo, amount of photosynthetically active vegetation, and soil moisture, respectively [55].

#### 2.1.3. Step 2—Extraction of EFAs

Next, metrics describing aspects of the intra-annual dynamics of ecosystem functioning (i.e., EFAs) have to be extracted from the SITS for each target year (Figure 1—"Step 2"). These types of metrics are commonly used to describe intra-annual aspects of ecosystem functioning related to quantity, seasonality, and timing (e.g., [5,46,48]). Examples of "quantity" metrics are the annual mean or median, the maximum or minimum values, or an integrated cumulative sum during a particular part of the year (e.g., the growing season). For "seasonality", metrics such as the intra-annual standard deviation or range can be used. Examples of "timing" metrics are the dates of the annual maximum, minimum, or other important moments (e.g., the start/end of the growing season).

#### 2.1.4. Step 3—Computation of EFA Anomalies

To obtain inter-annual EFA anomalies, it is necessary to first compute reference values for each pixel (Figure 1—"Step 3a"). To this end, EFAs are extracted for a multi-annual reference profile, which is obtained by averaging all years, after excluding values for years with any record of wildfire occurrences (this procedure is similar to the one used in [34,60], except for the exclusion of values in years with fire occurrences). Then, the anomalies are calculated by subtracting the reference values to those of each year, for each EFA (Figure 1— "Step 3b"). More precisely, in the case of linear (i.e., non-circular) metrics, the anomalies are calculated using the following expression:

$$A\_{\rm m} = X\_{\rm m} - R\_{\rm m\nu} \tag{1}$$

where *Am* represents the annual anomalies of EFA *m*; *Xm* is the annual values of EFA *m*, and *Rm* represents the values of the reference year for EFA *m*. In the case of circular metrics—e.g., a day of the year (DOY)—anomalies are calculated through the smallest differences between angles, in the following manner:

$$A\_{\mathfrak{m}} = \left(\frac{365}{2\pi} \times \theta\right) + 35,\tag{2}$$

where *Am* represents the annual anomalies of EFA *m*, and theta is defined as follows:

$$\theta = \begin{cases} \arctan\left(\frac{y}{x}\right), & \text{if } x > 0 \\ \arctan\left(\frac{y}{x}\right) + \pi, & \text{if } x < 0 \text{ and } y \gg 0 \\ \arctan\left(\frac{y}{x}\right) - \pi, & \text{if } x < 0 \text{ and } y < 0 \\ \frac{\pi}{2}, & \text{if } x = 0 \text{ and } y > 0 \\ \frac{-\pi}{2}, & \text{if } x = 0 \text{ and } y < 0 \\ \text{undefined}, & \text{if } x = 0 \text{ and } y = 0 \end{cases} \tag{3}$$

with:

$$\infty = \cos(X\_{\mathfrak{m}} - R\_{\mathfrak{m}}),\tag{4}$$

and:

$$y = \sin(X\_{\text{ll}} - R\_{\text{ll}}),\tag{5}$$

where *Xm* and *Rm* represent the same as in Equation (1).

#### 2.1.5. Step 4—EFA Ranking and Selection Procedures

Since multiple aspects of the intra-annual dynamics of ecosystem functioning can be obtained through the computation of EFAs, a selection process is necessary to reduce the initial set of metrics to an essential (i.e., relevant and less correlated) parsimonious set. In this selection process, classification models are used to rank all the EFAs considered in terms of their potential to measure the remotely sensed effects (as proxies of severity) of wildfire disturbances on the different EFAs (Figure 1—"Step 4a"). To this end, burned-area maps considered the response variable, while the anomalies of all EFAs are considered the predictive variables. The resulting ranking of EFA anomalies is used to compare all EFAs across (i) base SITS, with each corresponding to a different dimension of ecosystem functioning (i.e., "primary productivity", "vegetation water content", "albedo", and "sensible heat"); and (ii) the specific metric (e.g., mean vs. median) and (iii) component of the inter-annual dynamics of ecosystem functioning (i.e., "quantity", "seasonality", and "timing") used, each corresponding to an EFA. Furthermore, potential temporally lagged effects can also be evaluated by including values of EFA anomalies of different years concerning the year of the wildfire disturbance event (e.g., the anomalies for the year of the fire event, plus the first year after, the second year after, etc.).

The resulting ranking can be coupled with additional selection criteria, such as pairwise correlations, to minimize potential redundancy between EFAs so that a compact set of twelve EFAs (i.e., one for each dimension and each component of the inter-annual dynamics of ecosystem functioning) are selected (Figure 1—"Step 4b"). Then, both the magnitude and the direction of the pairwise relationship between the response variable (e.g., burned areas) and each of the selected EFAs are assessed (e.g., using Cohen's d effect size statistics [61]). This allows for the inter-comparison of selected EFAs across all the dimensions and components of ecosystem functioning, as well as other factors considered (e.g., the year of the early post-fire period). Based on these procedures, a final selection of EFAs is carried out to support the analysis of the main patterns of wildfire disturbance severity on the dimensions of ecosystem functioning considered.

#### 2.1.6. Step 5—Translation into Indicators of Wildfire Disturbance Severity

To support a multi-dimensional assessment of the effects of wildfire disturbances on ecosystem functioning, the final set of selected EFAs are then converted into indicators of wildfire disturbance severity (Figure 1—"Step 5"). For obtaining these indicators, the absolute values of the corresponding EFA anomalies are used, as these constitute measures of the magnitude of deviations from reference, thus allowing easier comparison across different dimensions. The resulting compact set of indicators is then used to support more detailed analyses.

#### *2.2. Test Case*

#### 2.2.1. Study Area

To illustrate our proposed approach, we chose the northwest of the Iberian Peninsula (NW-IP) as the overall study area (Figure 2), since this is a hotspot in terms of wildfire occurrences within the European context, being a highly fire-prone landscape [3]. This area extends to the north and the west to the Atlantic Ocean, while its southern boundary was delimited by the edge of MODIS tile h17v04. Its eastern boundary, on the other hand, was defined to include the northern half of mainland Portugal, as well as the Spanish autonomous community of Galicia.

The NW-IP features strong environmental gradients and includes a major biogeographic transition, from an Atlantic climate with temperate deciduous broadleaf and mixed forests in the north and northwest, to a Mediterranean climate with evergreen sclerophyllous vegetation towards the southeast. Land cover and uses are well diversified and highly heterogeneous (Table 1), including cropland areas, urban areas, plantation forests (mostly maritime pine, eucalyptus, and mixed stands), and shrublands in different successional stages, along with historical use of fire for agrosilvopastoral purposes. Wildfire occurrences are more frequent during the months between July and October—but especially in August—with burned areas over 10,000 ha being increasingly common. Wildfire events have been increasing in recent decades because of changes in land use, which mainly resulted from extensive abandonment of farming and husbandry, increasing the fuel load and its continuity in the landscape [30,62]. The frequency of extreme fire events, also known as mega-fires (i.e., complex catastrophic events in terms of size, intensity, resistance to control, severity, etc. [30]), has also been increasing, despite the huge investments in fire suppression [28], making this one of the regions with both the highest annual values of burnt area in Europe and the highest density of ignitions among southern European countries [63].

**Figure 2.** The study area for the illustrative test case, with the (**a**) location of the overall study area for regional-scale analyses within the Iberian Peninsula; (**b**) size distribution of all fires; (**c**) distribution of the area burned, by day of the year (DOY); and (**d**) location of the four individual burned patches selected for local-scale analyses (labels A–D), overlapped with the number of fires that occurred. For these figures, all wildfires with area > 100 ha which occurred between 2000 and 2018 (according to the MODIS MCD64A1 Burned Areas product) were considered.

To showcase our proposed framework, the full study area in NW-IP, including both burned and unburned areas of NW-IP between 2000 and 2018, was used for regional-scale analyses. Then, for local-scale analyses, four individual burned patches (see boxes labelled A–D in Figure 2) were selected from different geographical locations within the study area, featuring diverse baseline environmental characteristics (e.g., altitude, distance to coast, climate, proportions of major vegetation types), and with the focal wildfire event having occurred at different years within the considered period (A: 2003; B and C: 2005; D: 2013). Table 1 presents a summary of the characteristics of both the full study area, used for regional-scale analyses (i.e., NW-IP), as well as the four individual burned patches used for local-scale analyses.

#### 2.2.2. Satellite Data Preprocessing

In our test case, we used SITS extracted from MODIS products. Most tasks were undertaken within the R statistical programming environment [67], using mainly the "raster" package [68], with additional R packages complemented by other software for more particular tasks, as specified whenever relevant.

**Table 1.** Summary of baseline environmental characteristics of the overall study area for regionalscale analyses (northwest Iberian Peninsula, NW-IP), as well as the four selected burned patches selected for local-scale analyses (A–D). Elevation was extracted from MERIT DEM [64]; climate variables were extracted from the CHELSA Bioclim dataset [65]; the percentages of land-cover classes were extracted from CORINE Land Cover (CLC) 2000 [66].


We started by downloading and preprocessing all available images encompassing NW-IP (i.e., MODIS tile h17v04), between 2000 and 2018, from the following MODIS products, using the "MODIStsp" R package [69]:


Preprocessing consisted of converting to the GeoTIFF file format, re-projecting to a common reference system (WGS84/UTM29N), and resampling to 500-m pixel size using the nearest neighbor interpolation algorithm. Additionally, a filter based on the Hampel outlier identifier [73,74] with a window size of 5 was applied to the resulting SITS to reduce residual noise that hinders time-series data.

Next, LST and TCTB, TCTG, and TCTW were extracted from the MOD11A2 and MOD01A1 products, respectively, using "GDAL" [75] and the "rasterio" Python package [76]. For LST, only the day temperatures were used, which were then rescaled to degrees Celsius. The three TCT features of "Brightness", "Greenness", and "Wetness" were obtained by combining the 7 bands available in the MOD09A1 product, using the coefficients presented in Table 2 [59].

Burned areas extracted from the MCD64A1 product were filtered by size, eliminating all burned patches smaller than 100 ha, following a set of criteria detailed in previous work (see [44]). To further reduce the remaining noise in the data, a Whittaker–Henderson smoother [77,78] was applied to the four extracted SITS, with a lambda parameter equal to 2 (a rather low value, to conserve the characteristics of the time-series while still reducing extreme noise).


**Table 2.** MODIS-specific coefficients to calculate the Tasseled Cap Transformation (TCT) features "Brightness", "Greenness", and "Wetness". These features result from a rigid rotation of principal component axes, which are aligned with the biophysical parameters of albedo, the amount of photosynthetically active vegetation, and soil moisture, respectively [59].

NIR: near-infrared; SWIR: short-wavelength infrared.

#### 2.2.3. EFA Anomalies Computation

A total of 60 EFAs were extracted, using 15 different intra-annual metrics (Table 3), from each of the four SITS—LST, TCTB, TCTG, and TCTW—for each year between 2000 and 2018 for NW-IP. Among these, well-known measures were computed for "quantity"—mean, median, maximum, and minimum annual values—and "seasonality"—standard deviation, median absolute deviation, and absolute range. In addition to those, the relative range and a non-parametric relative range were also computed using the following expressions:

$$rrl = \log\_{10} \left( \left| \frac{r \text{ng}}{a v \text{g}} \right| \right) \tag{6}$$

and

$$mpr = \log\_{10}\left(\left|\frac{mg}{mdn}\right|\right),\tag{7}$$

where *rrl* is the relative range metric, while *rng* is the absolute range, *avg* is the mean, *npr* is the non-parametric relative range, and *mdn* is the median.

As for "timing", the times (of the year) of the maximum and minimum values were used. However, as variables of time of the year are circular by definition, we also tested two linearized versions for each, called "winterness" and "springness", using the following expressions:

$$\text{winterness} = \cos\left(\frac{2\pi}{365} \times (x - 35)\right) \tag{8}$$

and:

$$\text{springness} = \sin\left(\frac{2\pi}{365} \times (x - 35)\right),\tag{9}$$

where *x* is the day of the year (DOY) of either the maximum or minimum value (tmx or tmn, respectively). Note that these last two expressions (i.e., Equations (8) and (9)) include the conversion from DOY to radians, assuming that a full year has 365 days; as well as a rotational realignment of -35 days so that the maximum and minimum values of winterness and springness (i.e., 1 and −1, respectively) correspond to the mid-point of each season (in the case of NW-IP, as it is located in the northern hemisphere, the mid-point of winter is in early February). These metrics can be interpreted as corresponding to the temporal proximity, within the year, to either the start of winter or spring, respectively.

Finally, anomalies were then obtained for all EFAs for each of the four SITS by calculating either linear differences (Equation (1)), for all metrics except "tmx" and "tmn", or the smallest differences between angles, for those two circular metrics (Equations (2)–(5)).


**Table 3.** List of metrics extracted from each of the four satellite image time-series used in this study (i.e., land surface temperature and the Tasseled Cap Transformation features of "Brightness", "Greenness", and "Wetness"); as well as the number of predictive variables for models extracted from each one.

#### 2.2.4. Ranking and Selection of EFAs

A ranking of EFA anomalies was obtained for NW-IP through a classification modeling procedure, in which the response variable was the binary burned class (i.e., burned vs. unburned) extracted from the MODIS burned area product. As predictive variables for these models, we used the inter-annual anomalies in each one of the 60 EFAs (i.e., each of the 15 intra-annual metrics described in Table 3, for each of the four SITS), for each one of three years since the fire occurrence: the year of the fire, as well as the first and second years after (i.e., years 0, +1, and +2). This was to allow for the evaluation of potential temporally lagged effects (for which only the years from 2003 to 2016 were used). In total, the number of predictive variables used in these models was 180 (see Supplementary Information S1 for additional information on the predictive variables used in these models). To this end, the Random Forest (RF) "ensemble learning" algorithm was used, since this technique scales well on large and/or high-dimensional data [79], thus allowing for the full set of variables to be tested without the need for a selection before modeling.

As the number of pixels identified as "burned" represented only ca. 1.3% of the total overall study area (i.e., NW-IP), data balancing through down-sampling of the majority class (i.e., "unburned") was applied. To that end, 100 random sample sets were defined, consisting of all the pixels identified as "burned", as well as an equal number of randomly selected "unburned" pixels for each year, in a total of 122,658 observations (i.e., n = 61,329, per class). Then, an RF classification model was calibrated for each one of the 100 sample sets, using the "caret" R package [80], with the hyperparameters set to the default values (including a fixed value for "mtry" equal to the default squared root of the number of predictive variables), and under a leave-group-out cross-validation (LGOCV) scheme with 50:50 train/test split in each one of 10 repetitions.

Model performance was then evaluated using measures specific to two-class problems sensitivity (or true positive rate), specificity (or true negative rate), and the area under the receiver operating characteristic curve (AUC)—which were then averaged across the 100 resulting models. Finally, variable importance scores—a common approach to extract interpretable information on the contribution of different variables from RF models [79] based on the mean decrease in accuracy, were extracted and also aggregated to the median values across the 100 sample sets to support the EFA ranking. Based on the ranking of variable importance scores obtained from the RF models, one EFA was selected for each

dimension and each component, giving a total of 12, to support the analysis of the main patterns of wildfire disturbance severity on multiple dimensions of ecosystem functioning in NW-IP.

#### 2.2.5. Analysis of Indicators of Wildfire Disturbance Severity

To support the analysis of both the directionality (i.e., increasing or decreasing) and magnitude of the pairwise relationships between the regional-scale burned areas and each of the selected EFAs, Cohen's d effect size statistics [61] and Spearman pairwise rank correlations were calculated for the 100 sample sets. From this, as well as from the distributions of the selected EFA anomalies, the main overall patterns of the effects of wildfire disturbances on ecosystem functioning were analyzed for NW-IP and compared across the four dimensions and three components of ecosystem functioning, as well as across years of the early post-fire period.

It should be noted that effect sizes were applied here to avoid the use of statistical hypothesis testing based on significance (i.e., *p*-values), which have important known limitations such as not measuring the effect size or the importance of a result [81] or the strong influence of sample size, which is hugely influential in determining significance levels [82] (i.e., when sample sizes are large enough, any effect can produce a small *p*-value if the sample size or measurement precision is high enough [81]; with large sample sizes, often happening at a pixel-level analysis of SITS).

We then computed four indicators of the effects of wildfire disturbances on ecosystem functioning (i.e., one for each dimension) for four selected burned areas (see Figure 2) at a local scale to illustrate the proposed approach for RS-based assessment, mapping, and monitoring of wildfire disturbance severity on multiple dimensions of ecosystem functioning. The selection of EFAs for these final indicators took into consideration the results from both the Random Forest model-based ranking and the pairwise correlations between EFAs. The translation of the selected EFAs into indicators was done by taking the absolute values of the anomalies of the selected EFAs for the four burned areas. Finally, spatial patterns between all indicators, across years, were further analyzed, local-scale, through local Spearman correlations, with a neighborhood window size of 5 pixels.

#### **3. Results**

#### *3.1. EFA Ranking*

In terms of the performance of the Random Forest (RF) models for our test case, the results obtained for leave-group-out cross-validation can be considered very good across all sample sets, with around 0.9815 ± 0.0004 for AUC, 0.9080 ± 0.0015 for sensitivity, and 0.9647 ± 0.0010 for specificity.

As for variable importance scores, EFAs extracted from TCTG had overall higher values than the other dimensions (Figure 3; see Supplementary Information S2 for the full plot of variable importance scores for all non-aggregated variables). Considering each of the three different EFA components, the ones measuring aspects of "quantity"—particularly the mean and either one of the extreme values (minimum or the maximum, depending on the specific case)—obtained overall higher variable importance scores than "seasonality" or "timing" EFAs, with the notable exception of seasonality EFAs extracted from TCTG at year 0. Conversely, the same exception was observed within all seasonality EFAs, where otherwise, LST was the dimension from which EFA anomalies obtained higher variable importance scores overall, especially for standard deviation. As for timing EFAs, the scores obtained were generally lower than for the two other groups, except for the ones extracted from either TCTG or LST at year 0, particularly the springness of the time of the year of the minimum TCTG value and the winterness of the time of the year of the maximum LST value, closely following the corresponding quantity EFAs.

Across the four different dimensions of ecosystem functioning, those extracted for the year of fire (i.e., year 0) generally obtained better scores than those of the two subsequent years. However, for quantity EFAs extracted from TCTW, the first year after the fire (i.e., year +1) obtained better scores. Furthermore, quantity EFAs extracted from TCTB also achieved higher scores for the second year after the fire (i.e., year +2), although not as high as for year 0. Finally, in terms of individual EFAs, the top-ranked ones were extracted from TCTG, with the minimum value (TCTG-min) and the standard deviation of TCTG (TCTG-std) at year 0 holding the highest scores among all EFAs tested (see Supplementary Information S1 for additional information on the ranking results).

**Figure 3.** Distributions of variable importance scores obtained from Random Forest models to assess the effects of wildfire disturbances on four dimensions of ecosystem functioning: primary productivity ("productivity"), vegetation water content ("water"), albedo, and sensible heat ("heat"), in the same year as the respective fire ("year 0") and the two years after ("year +1" and "year +2"), for each of the three components considered: quantity (**left**); seasonality (**center**), and timing (**right**) (Note that the scale of the *y*-axis is logarithmic).

#### *3.2. Analysis of Effects*

Considering the RF-based importance ranking, 12 EFAs (Table 4) were selected, one for each dimension (i.e., "primary productivity", "vegetation water content", "albedo", and "sensible heat") and each component of intra-annual dynamics (i.e., "quantity", "seasonality", and "timing"). As for "timing" EFAs, it should be noted that the selected metrics do not correspond to the highest-ranked ones for each case. In this instance, parsimony and interpretability criteria were also considered, since linearized versions of timing EFAs (i.e., "winterness" and "springness") often require the use of both at the same time to better understand in which month/season the minima/maxima occur (see Supplementary Information S3 for partial dependence plots of the predictive variables of the models, extracted from the 12 selected EFAs).

Table 4 summarizes the results for the effect sizes (both magnitudes and directions) between the response variable and each of the selected EFA anomalies (see Supplementary Information S4 for additional information on the effect size results). These showed overall strong negative effects (i.e., lower values associated with wildfire disturbances) for TCTG-min and TCTW-avg (especially at years 0 and +1 and years +1 and +2, respectively), while positive (i.e., raised values associated with wildfire disturbances) for LST-max (especially at years 0 and +1). For TCTB-avg, however, the effects changed from strongly negative at year 0 to strongly positive at year +2. As for seasonality EFAs, the effects were, overall, positive, with the strongest magnitude associated with TCTG-std and LST-std. Finally, timing EFAs had overall weaker effects, except for TCTG-tmn (with large negative effects) and LST-tmx (with large positive effects), both only at year 0.

**Table 4.** List of selected EFAs with the corresponding effects categories, based on Cohen's d effect size statistic. The number of arrows in the "Effect category" columns symbolizes the magnitude of the effect size, with |d| < 0.2 "negligible", |d| < 0.5 "small", |d| < 0.8 "medium", and otherwise "large". The direction and color of the arrows symbolize the sign of the effect size, with ascending arrows depicting positive (non-"negligible") effect sizes while descending arrows depict negative (non-"negligible") effect sizes ("negligible" effects are symbolized by a horizontal dash).


#### *3.3. Main Patterns in EFA Anomalies*

The distributions of the deviations from the normal variation (i.e., anomalies) for each of the twelve selected EFAs, grouped by year after fire (i.e., years 0, +1, and +2) and by burned vs. unburned pixels, are shown in Figure 4 (see also Supplementary Information S5 for the distributions of the reference values of all pixels in the NW-IP).

Overall, there were notable differences in the distributions obtained for burned vs. unburned pixels, with the most notable contrast having been observed for TCTG-tmn at year 0, with median anomalies corresponding to occurrences of 80–88 days earlier than the reference. Within the same dimension, differences were also notable for TCTGmin (especially at years 0 and +1), and also for TCTG-std (especially at year 0), with values of burned pixels considerably below/above (respectively) those of unburned pixels. Considerable differences between burned vs. unburned were also obtained for TCTW-avg, but only at years +1 and +2, while for TCTB-avg, they were observable at years 0 and +2, with opposite directions (i.e., negative vs. positive). As for heat-related EFAs, differences in the distributions of anomalies between burned and unburned pixels were notable for LST-max and LST-std across all years after fire, but especially at years 0 and +1. Finally, LST-tmx showed a contrast between burned and unburned pixels, with median anomalies corresponding to delays of 8–16 days.

#### *3.4. Multi-Dimensional Assessment of Wildfire Disturbance Severity*

For the RS-based assessment of wildfire disturbance severity on the four dimensions of ecosystem functioning considered, at the local level, one indicator for each dimension was used, based on the previously selected anomalies in quantity EFAs (i.e., TCTG-min, TCTWavg, TCTB-avg, and LST-max). This selection took into consideration the results from both the model-based ranking (i.e., quantity EFAs ranked better overall than those of seasonality or timing) and pairwise correlations (up to ±0.83; see Supplementary Information S6).

**Figure 4.** Distributions of the anomalies of the 12 selected EFAs, in the same year as the respective fire ("year 0") and the two years after ("year +1" and "year +2"), for each of the three components considered: (**a**) quantity; (**b**) seasonality, and (**c**) timing. The specific intra-annual metrics selected for each of the four dimensions of ecosystem functioning considered were the following: Primary productivity ("Productivity")—Tasseled Cap Transformation "Greenness" feature minimum value (TCTG-min), standard deviation (TCTG-std), and time of the year of the minimum value (TCTG-tmn); Vegetation water content ("Water")—Tasseled Cap Transformation "Wetness" feature average value (TCTW-avg), standard deviation (TCTW-std), and time of the year of the minimum value (TCTW-tmn); Albedo—Tasseled Cap Transformation "Brightness" feature average value (TCTB-avg), standard deviation (TCTB-std), and time of the year of the maximum value (TCTB-tmx); Sensible heat ("heat")—land surface temperature maximum value (LST-max), standard deviation (LST-std), and time of the year of the maximum value (LST-tmx).

Figure 5 shows the median (±0.5 median absolute deviations) profiles of each of the four selected indicators, across the three years after fire, for each of the four burned areas of interest (see Supplementary Information S7 for additional information on the spatial and temporal patterns of the indicators of wildfire disturbance severity in the four burned areas A–D). This figure illustrates different post-fire trajectories, with varying relations between the four dimensions of ecosystem functioning across individual burned areas. Specifically, different areas had different ratios between "Productivity" and "Heat" (e.g., A vs. D). Furthermore, "Water" seemed to be less affected in area D than in the other areas, while "Albedo" in area D exhibited a different post-fire trajectory relative to the other areas, with the value for year +1 being slightly higher than for year 0.

Finally, Figure 6 shows the variation in the spatial and temporal patterns of each of the four selected indicators of primary productivity, vegetation water content, albedo, and sensible heat, for each of the four individual burned areas of interest (A–D). Overall, some wildfire disturbance severity "hotspots" are observable, while not necessarily coinciding

across dimensions (see also Supplementary Information S7). Moreover, different burned areas exhibited different patterns across time between dimensions of ecosystem functioning, with temporal trends being divergent in some cases while convergent in others (e.g., "Water" and "Albedo" dimensions in the areas of interest C vs. D), or with different relative preponderances (e.g., "Productivity" vs. "Heat" in the areas of interest A vs. D), further illustrating the patterns observed in Figure 5.

**Figure 5.** Profiles of wildfire disturbance severity on the four dimensions of ecosystem functioning: Primary productivity ("Productivity"), vegetation water content ("Water"), albedo, and sensible heat ("Heat"), at short to medium term (i.e., between years 0 and +2 after the fire event), for four individual burned areas (letters A–D; see Figure 2 for their location within the study area). Dots connected by thick lines represent median values across all pixels of each individual burned area, while shaded areas of the corresponding color represent ±0.5 × median absolute deviation. Values for the "heat" dimension are scaled by a factor of 0.005 for visual comparability purposes.

**Figure 6.** *Cont.*

**Figure 6.** Maps of wildfire severity on primary productivity ("Productivity"), vegetation water content ("Water"), albedo ("Albedo"), and sensible heat ("Heat"), at short to medium term (i.e., between years 0 and +2 after the fire event), for four individual burned areas (letters **a**–**d**; see Figure 2d, labels A–D), based on indicators derived from inter-annual anomalies in quantity EFAs extracted from satellite image time-series.

#### **4. Discussion**

#### *4.1. Fire Severity Patterns in the NW Iberian Peninsula*

#### 4.1.1. Effects of Wildfires across Dimensions and Components

In our test case, we found that the associations between burned areas and inter-annual EFA anomalies were particularly strong for primary productivity and sensible heat, although albedo and vegetation water content also revealed important effects caused by wildfire disturbances. This is in line with previously published research, as the sudden removal of green vegetation due to fire usually translates into abrupt breaks that are commonly observable in time-series of spectral vegetation indices (e.g., TCTG). This removal can result in substantial carbon losses, both in terms of aboveground biomass [8] as well as soil organic matter [9,83]. Moreover, fire-mediated changes in nutrient concentrations can ultimately limit productivity over long periods [10]. Because vegetation is a regulator of land surface energy fluxes [27], the removal of vegetation causes an increase in observed LST induced by a reduction in evapotranspiration [25], thus increasing the sensible-tolatent heat ratio [17]. Furthermore, the observed effects in albedo are also consistent with studies reporting either a darkening or a brightening (or both) effect after a fire (e.g., [22,84]) due to the vegetation removal as well as the presence (or not) of black carbon in soot, which absorbs visible solar radiation [84]. Finally, effects on vegetation water content could be related to loss of moisture in canopy foliage due to damage caused by fire, eventually leading to vegetation mortality [15]. Indeed, our results show that using information from multiple dimensions of ecosystem functioning can yield improved results over the use of fewer—or only one—indices for enhanced detection, mapping, and evaluation of ecological effects of disturbances such as wildfires.

Across the three different components of ecosystem functioning considered, anomalies in quantity EFAs (e.g., mean, minimum, maximum) had overall stronger associations with wildfire disturbances than seasonality (e.g., standard deviation) or timing (e.g., time of minimum or maximum) ones. For severity assessment purposes, quantity and seasonality EFAs seemed to be somewhat redundant, as there can be, in some cases, high pairwise correlations between EFA anomalies of those two types. Conversely, for other types of disturbances such as land clearing for agriculture and ranging, the impact on seasonality (i.e., seasonal variability) can be much stronger than on quantities (see [85]). While some studies show that wildfire disturbances can have profound impacts on aspects of timing of key dimensions of ecosystem functioning (e.g., by inducing phenological changes [11]), our approach may not be best-suited to assess such modifications. This is because the intra-annual descriptors used (i.e., EFAs) seemed to more likely capture the location of extreme values (i.e., minimum/maximum) within the year corresponding to direct effects of the wildfire disturbance rather than translating changes (e.g., delays or lags) in the annual cycles of the underlying ecological processes.

As for the specific statistical measures used to extract EFAs, stronger effects can result from the use of parametric measures over that of their non-parametric equivalents. This may be because parametric measures are generally more sensitive (and susceptible) to extreme or abnormal values than non-parametric ones, which may or may not be desirable, depending on the specific application purpose. On the other hand, non-parametric measures may be more adequate for extracting long-term inter-annual trends than parametric ones (e.g., [5]).

The results obtained point to a strong added value of the proposed approach to discriminate the effects of wildfire disturbances on different aspects of ecosystem functioning.

#### 4.1.2. Temporal Effects of Wildfires

Short-term effects of wildfire disturbances on ecosystem functioning seemed to attain better overall performance than medium-term effects, since inter-annual anomalies in intraannual EFAs for the year of the fire occurrence obtained much better scores, overall, than those for the first and second years following fire. This is partially in line with other studies, as the majority of studies analyzing RS-based effects of fires focus on either short-term (i.e., abrupt) wildfire disturbance severity (e.g., [32,56]) or on longer-term effects more related to post-fire recovery (e.g., [3,34]), while the medium-term effects (i.e., from one or two years to five years after the fire) are often overlooked.

EFAs derived from TCTW (i.e., vegetation water content) seemed to be more affected by wildfire disturbances in the first year after the fire than the year in which the fire occurred, suggesting a temporal lag associated with this dimension relative to both the fire event and other dimensions. This observation (which, to our knowledge, has not been previously reported) could be linked to changes in hydrological dynamics, such as increased soil water repellency, precipitation run-off, and increased quantities of impervious materials such as ashes and soot [84], which can clog soil pores [86]. Alternatively, this effect could be due to post-fire changes in foliar moisture leading to mortality occurring over an extended period after fire [15], since TCTW is quite sensitive to the moisture content of vegetation [59].

As for the observed effects of wildfire disturbances on albedo (i.e., TCTB), this dimension of ecosystem functioning seemed to exhibit a shift in directionality, from negative in year 0 to positive in years +1 and +2, denoting a change in the relationships with the wildfire disturbance at short vs. medium term. This is also in line with the findings reported in other studies [22,87], which observed both a decrease in albedo (i.e., a darkening effect) immediately after the fire, and, usually, for a brief period, as well as a small increase (i.e., brightening) one year after the fire, which can be persistent. On the other hand, this darkening, which can be due to changes in the relative abundance of surfaces with distinct reflective properties [88] (e.g., ash, char, soot, bare soil), can, in some cases, persist for multiple years after a fire [23].

Finally, the effects of wildfire disturbances on both primary productivity and sensible heat (through TCTG and LST, respectively) seemed to be reflected across the three time periods analyzed. However, they were especially noticeable for the year of the fire occurrence. This also goes in the same direction as other studies addressing the impacts of fire on ecosystems, which report usually short-term effects for satellite-derived proxies of primary productivity (e.g., [41,89]) after fire. An increase in observed LST immediately following fires is also frequently reported in the literature [25]. This can sometimes still be observed one year following the fire [19], after which LST anomalies tend to become smaller, and seasonality starts governing the LST time-series [27].

Together, these results highlight the importance of taking into consideration both shortand medium-term effects of wildfire disturbances on multiple dimensions of ecosystem functioning.

#### 4.1.3. General Patterns

Overall, EFA anomalies were able to efficiently discriminate burned and unburned areas, with very high values obtained for both effect sizes and performance measures of the Random Forest (RF) models between EFA anomalies (predictive variables) and burned areas (response variable). This allowed supporting the selection of EFAs from which a compact set of indicators were employed to assess the main patterns of wildfire disturbance severity in our test area.

Diverse spatio-temporal patterns of wildfire disturbance severity on the four dimensions of ecosystem functioning—primary productivity, vegetation water content, albedo, and sensible heat—could be observed both between as well as within individual burned areas across years (i.e., from the year of the fire to two years after). These results showcase the added value of (i) analyzing the effects of wildfire disturbances for multiple dimensions of ecosystem functioning, as opposed to addressing only overall burn severity in commonly used approaches (e.g., based on the Differenced Normalized Burn Ratio, ΔNBR); and (ii) taking into consideration both the short and medium term, thus enabling to obtain insights on potentially lagged effects of fires on ecosystem functioning, as different dimensions can exhibit different severities, and with different timings.

#### *4.2. General Considerations about the Proposed Framework*

#### 4.2.1. Satellite Image Time-Series

While in the presented test case, we used data extracted from MODIS products exclusively, other sources of satellite image time-series (SITS) are available and can equally be applied using the proposed approach. Such data sources include platforms with diverse characteristics in terms of spatial, temporal, and spectral resolutions and length of the historical archive, such as SPOT-VEGETATION, PROBA-V, Landsat 5/7/8, or Sentinel-2/3 (e.g., [90]).

As for the base SITS to be extracted from those sources, we recommend, in the proposed approach, the use of LST and the TCT features of "Brightness", "Greenness", and "Wetness" as RS-based proxies of albedo, primary productivity, and vegetation water content, respectively. This option allows to derive information on the four above-mentioned dimensions of ecosystem functioning, from the smallest number of sources (the three TCT features are derived from the same source), compactly and coherently, and is often available for the major satellite data sources (e.g., [91]). However, in the case of primary productivity and vegetation water content, if data are not available to derive all four RS-based variables, alternative spectral indices could be used instead, such as the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI) for "productivity", or the Normalized Difference Water Index (NDWI) or Land Surface Water Index (LSWI) for "water".

#### 4.2.2. Additional Data Sources

Well-designed field-based measurements could be used, if available, to validate the information provided by SITS, further enhancing the proposed approach, although it is not a mandatory requirement. To that end, data can be collected in the field, following robust sampling designs, such as spectral/radiometric readings, aerial surveys using cameras and sensors mounted on unmanned aerial vehicles (UAVs), or data collected under already existing and commonly used protocols (or adapted versions from these), such as the Composite Burn Index (CBI; e.g., [92]) or the improved geometrically structured Composite Burn Index (GeoCBI; e.g., [93]), to complement the information derived from SITS. Such data can, however, be difficult and/or costly to obtain, even more so if multiple observations are required to produce a coherent time-series. These types of field-level measurements are usually compared to commonly used indicators of burn severity (e.g., [40]) such as the differenced normalized burn ratio (ΔNBR)—or any of its enhanced forms, such as the relative normalized burn ratio (RdNBR) and the relativized burn ratio (RBR).

Unlike the more commonly used RS-derived methods to estimate burn severity, the approach proposed in this study can provide information on the effects of wildfire disturbances on multiple dimensions of ecosystem functioning—namely primary productivity, vegetation water content, albedo, and sensible heat—simultaneously, towards a comprehensive evaluation of fire severity.

#### 4.2.3. Applicability and Future Directions

Multi-dimensional approaches to wildfire disturbances using EFA-like metrics are, to our knowledge, scarce (e.g., [94]). Very few studies have adopted a multi-dimensional approach using EFAs, regardless of the specific application (e.g., [48,50,95]). More studies should be conducted in the future, with diverse contexts in terms of baseline environmental conditions and fire disturbance regimes, to further test and improve the proposed approach. It is important to highlight that although the results obtained for the test case presented in this study may be somewhat specific to the study area and the analyzed period, the general approach could be applied, with relative ease, to other contexts and timeframes to conduct informed selections of the most efficient and informative EFA-based indicators for wildfire disturbance severity assessment.

The proposed approach makes use of annual descriptors of the dynamics of multiple dimensions of ecosystem functioning—so-called Ecosystem Functioning Attributes (EFAs) extracted from SITS, to derive indicators of wildfire disturbance severity. However, another possible approach could be to derive indicators directly from the time-series, without the need to compute EFAs first. This would open the possibilities for other kinds of indicators, with finer temporal (i.e., <1 year) resolutions, on a (semi-)continuous basis, to be derived from SITS (e.g., [3]). On the other hand, although additional types of indicators could be obtained in this way, others would not—e.g., metrics of timing based on annual cycles, such as the date of the maximum or minimum, would only make sense on an annual basis. Besides, the annual nature of EFAs—and the derived indicators—facilitates both their computation as well as their interpretation, as it offers a better translation of spectral indices into informative ecosystem variables. This can be advantageous for some application purposes, such as for reporting from official authorities. Notwithstanding, approaches adopting either an annual or continuous basis could complement each other to characterize the effects of wildfire disturbances on ecosystem functioning.

Finally, in this study, we focused on short- to medium-term impacts of wildfire disturbances. However, long(er)-term effects of wildfire disturbance severity have been observed (e.g., [96]). Indeed, long-term impacts of wildfire disturbances on ecosystem functioning could also be monitored using the approach proposed in this study, provided that longenough SITS are available. However, this was not tested in this study, as that would considerably increase both computational requirements as well as methodological complexity.

#### **5. Conclusions**

Here, we described a framework to assess the effects of wildfire disturbances on four key dimensions of ecosystem functioning—primary productivity, vegetation water content, albedo, and sensible heat. This approach is based on indicators derived from several descriptors of the intra-annual dynamics of ecosystem functioning—called Ecosystem Functioning Attributes (EFAs)—extracted from remotely sensed satellite image time-series

(SITS). We found that all four dimensions of ecosystem functioning suffered important effects of wildfire disturbances—especially primary productivity and sensible heat. We also found that quantity EFAs held the highest potential for indicating wildfire severity. Finally, we found important differences in short- and medium-term effects between the four different dimensions of ecosystem functioning, suggesting temporally lagged effects in vegetation water content as well as directionality shifts in albedo. Together, these results highlight the added value of the proposed framework to enhance RS-based assessment, mapping, diagnostics, and monitoring of wildfire disturbances. Using this approach, a more comprehensive evaluation of fire severity can be attained, which is not captured by indicators derived from only one spectral index, measuring only overall burn severity. As such, this multi-dimensional framework can contribute to a deeper and more detailed understanding of the impact of fires in ecosystems, with the potential to support assessments of severity, recovery, and resilience and their interactions with biodiversity and ecosystem services.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2072-429 2/13/4/780/s1: Supplementary Information S1—Summary table of results of the Random Forest models, for all 180 variables obtained from EFA anomalies, tested against burned vs. unburned areas; Supplementary Information S2—Distributions of variable importance scores from Random Forest models, for individual predictive variables; Supplementary Information S3—Partial dependence plots of the 36 predictive variables extracted from the 12 selected EFAs; Supplementary Information S4—Summary table of results, for each of the twelve EFAs selected for further analyses, tested against burned vs. unburned areas, of both the Random Forest models, and the effect sizes analysis; Supplementary Information S5—Distributions of the reference values (i.e., inter-annual median of values from unburned years) of the selected EFAs; Supplementary Information S6—Summary of maximum Spearman rank pairwise correlations between the 12 selected EFA metrics; Supplementary Information S7—Main spatial and temporal patterns obtained for each of the four burned areas of interest.

**Author Contributions:** Conceptualization, B.M., J.G., M.C., D.A.-S. and J.P.H.; methodology, B.M., J.G.; formal analysis, B.M.; data curation, B.M.; writing, B.M., with the assistance of J.G., M.C., D.A.-S. and J.P.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Portuguese national funds through FCT—Foundation for Science and Technology, I.P., under the GreenRehab project (PCIF/RPG/0077/2017). This work was partially supported by the collaboration between D.A.S. and J.H. under the projects RESISTE (P18-RT-1927), funded by Consejería de Economía, Conocimiento, Empresas y Universidad from the Junta de Andalucía, and DETECTOR (A-RNM-256-UGR18), with the contribution of the European Union Funds for Regional Development. B.M. was supported by FCT, the Ministry of Education and Science, and the European Social Fund, within the 2014–2020 EU Strategic Framework, through FCT (PhD scholarship SFRH/BD/99469/2014). J.G. was funded by the Individual Scientific Employment Stimulus Program (2017), through FCT (contract nr. CEECIND/02331/2017).

**Acknowledgments:** We acknowledge the use of MODIS imagery obtained from NASA's Land Processes Distributed Active Archive Center (LP DAAC), available free of charge. The authors would like to thank the anonymous reviewers for providing comments and suggestions that helped to improve the quality of the original manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


### *Article* **Wildland Fire Tree Mortality Mapping from Hyperspatial Imagery Using Machine Learning**

**Dale A. Hamilton \*, Kamden L. Brothers, Samuel D. Jones, Jason Colwell and Jacob Winters**

Department of Mathematics and Computer Science, Northwest Nazarene University, 623 S University Blvd, Nampa, ID 83686, USA; kbrothers@nnu.edu (K.L.B.); samueljones@nnu.edu (S.D.J.); jcolwell@nnu.edu (J.C.); jwinters@nnu.edu (J.W.)

**\*** Correspondence: dhamilton@nnu.edu; Tel.: +1-(208)-467-8769

**Abstract:** The use of imagery from small unmanned aircraft systems (sUAS) has enabled the production of more accurate data about the effects of wildland fire, enabling land managers to make more informed decisions. The ability to detect trees in hyperspatial imagery enables the calculation of canopy cover. A comparison of hyperspatial post-fire canopy cover and pre-fire canopy cover from sources such as the LANDFIRE project enables the calculation of tree mortality, which is a major indicator of burn severity. A mask region-based convolutional neural network was trained to classify trees as groups of pixels from a hyperspatial orthomosaic acquired with a small unmanned aircraft system. The tree classification is summarized at 30 m, resulting in a canopy cover raster. A post-fire canopy cover is then compared to LANDFIRE canopy cover preceding the fire, calculating how much the canopy was reduced due to the fire. Canopy reduction allows the mapping of burn severity while also identifying where surface, passive crown, and active crown fire occurred within the burn perimeter. Canopy cover mapped through this effort was lower than the LANDFIRE Canopy Cover product, which literature indicated is typically over reported. Assessment of canopy reduction mapping on a wildland fire reflects observations made both from ground truthing efforts as well as observations made of the associated hyperspatial sUAS orthomosaic.

**Keywords:** mask region-based convolutional neural network; small unmanned aircraft system; canopy cover; tree mortality

#### **1. Introduction**

This study investigates the use of machine learning to identify burn severity as a measure of canopy reduction. Tree crowns are classified from hyperspatial (sub-decimeter resolution) imagery using a mask region-based convolutional neural network (MR-CNN). The tree crown classification is then summarized at 30-m resolution, resulting in a canopy cover map. A post-fire canopy cover mapping product was compared to a pre-fire canopy cover mapping layer from the LANDFIRE project. The difference between the pre-fire and post-fire canopy cover enables a calculation of canopy reduction resulting as an effect of the fire, which gives an indication of tree mortality which is a primary indicator of post-fire effects resulting from a wildland fire.

Forest management practices, including fire suppression efforts, have led to forests with higher density than were found under pre-European settlement conditions. As a result, wildlands are experiencing a higher frequency of catastrophic fires, with millions of acres burned across the western US every year, often burning with much higher severity than experienced under historic conditions. These higher intensity fires result in increased incidence of soil erosion, hydrologic degradation, and deteriorating ecosystem resilience. The ability to detect trees in hyperspatial imagery enables the calculation of canopy reduction, which is a measure of wildland fire severity. The availability of canopy reduction data enables better decision making for ecosystem management and fire effect mitigation.

**Citation:** Hamilton, D.A.; Brothers, K.L.; Jones, S.D.; Colwell, J.; Winters, J. Wildland Fire Tree Mortality Mapping from Hyperspatial Imagery Using Machine Learning. *Remote Sens.* **2021**, *13*, 290. https://doi.org/ 10.3390/rs13020290

Received: 2 December 2020 Accepted: 8 January 2021 Published: 15 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### *Background*

The term "wildland fire severity" can refer to many different effects observed through a fire cycle, from how intensely an active fire is burning to the response of the ecosystem to the fire over the subsequent years. As opposed to the intensity of an active fire, which is often referred to as fire severity, this study investigates direct or immediate effects of a fire as observed in the days and weeks after the fire is contained, which we refer to as burn severity [1,2]. Previous studies in this effort examined the effect of the fire on biomass consumption which indicates how much surface vegetative material has been consumed by the fire [1,3–5]. By contrast, this study investigates the effect that fire has had on canopy vegetation in forested ecosystems. In situations where the amount of current post-fire canopy has not been significantly altered from what can be observed from the pre-fire conditions, it can be assumed that a fire was only consuming surface vegetation (or possibly the lower canopy) but did not consume the upper portion of the tree crowns in the canopy [6]. This low intensity fire low intensity surface fire results in no mortality of the trees in the upper canopy would result in no noticeable change of canopy cover between pre-fire and post-fire conditions. Another situation results from where all the upper tree crowns in the forest canopy have been either been consumed by a fire or the crown has been scorched by a more intense fire on the surface which extends into the tree crowns [6]. In these active crown fires, the canopy has been consumed or scorched by the fire, resulting in the death of all the trees provided that the trees are not resprouting deciduous species. This high severity active crown fire will result in a complete reduction of canopy cover as compared from pre-fire conditions. In between surface fires and crown fires, lie a situation where fire burning through surface vegetative material having enough intensity to burn or scorch the crowns of individual trees or clusters of trees, resulting in a moderate severity passive crown fire where some of the tree crowns are completely consumed or scorched [7]. While the post-fire canopy cover will be reduced from pre-fire canopy conditions as a result of these passive crown fires, the post-fire canopy cover will not be reduced to zero because some of the trees that made up the pre-fire canopy will survive, having been neither burned nor scorched.

Most land managers that use canopy cover mapping layers currently rely on the canopy cover layers that are available through the LANDFIRE project, a collaborative effort between the United States Forest Service and the US Department of the Interior. LANDFIRE produces geospatial products that describe wildland fire fuel conditions across the United States. From the LANDFIRE portal, land managers have access seamless nation-wide geospatial layers with a spatial resolution of 30 m (30 m) [8]. These canopy cover layers were generated from imagery acquired by the Enhanced Thematic Mapper Plus sensor on board the Landsat 7 satellite using a decision tree algorithm that was trained using plot data stored in the LANDFIRE Reference Database which was collected from across the US [9].

The utilization of small unmanned aircraft systems (sUAS) as a remote sensing tool provides land managers an affordable way to acquire imagery with much higher spatial resolution (sub-decimeter) than was previously available. This hyperspatial imagery produces homogeneous pixels which are comprised of a single class, resulting in higher mapping accuracy. By contrast, heterogeneous pixels which are large enough to spatially contain more than one class are common with lower resolution imagery [1,4].

Burned vegetation, unburned canopy vegetation (conifer and deciduous trees), and unburned surface vegetation (brush and grass) have been shown to have distinct spectral responses in both the visible and near infra-red spectra, though the spectral difference between canopy vegetation and surface vegetation is not as large as that between burned and general unburned vegetation [10]. This spectral difference is apparent when observing the two unburned vegetation classes in hyperspatial imagery. Persons who had participated in the acquisition of the sUAS imagery found that they could observe the visual difference between surface and canopy vegetation, with this difference in vegetation type observed in the imagery matching what they had observed in the field while conducting image acquisition flights with the sUAS.

Hamilton [1] found that tree crowns obstruct surface features in imagery acquired with a sUAS, preventing classification of surface features. Obstructed features include surface biomass consumption, a primary indicator of surface fire severity. When a tree crown is surrounded by burned surface fuel in the image, it is reasonable to assume that the ground under the tree has also burned but is obscured by the unburned crown [1,11]. Unburned tree crowns within a burned area cause an under-reporting of fire extent due to the misclassification of unburned tree crowns within the fire perimeter as being unburned areas. The first step to identifying and mapping obscured burned surface fuels is to identify trees within the image. When a tree surrounded by burned surface fuels is identified, it can be inferred that the surface fuels under the tree are in fact burned and should be included within the fire extent. Wildland fire can also affect an ecosystem through tree mortality. Mortality may result from either the ignition and consumption of a tree crown or from scorching of the tree crown caused by heat from surface fire beneath the tree. The extent of tree mortality can be estimated by calculating the change in canopy cover before and after the fire. Post-fire canopy cover can be calculated as the percentage of an image that consists of tree crowns. Pre-fire canopy cover can be retrieved from the Canopy Cover rasters that were published by the LANDFIRE Fire Fuels project [8].

#### **2. Materials and Methods: Mapping Canopy Cover and Tree Mortality**

Canopy cover describes the percentage area occupied by the upper tree crowns in forested ecosystems. In this effort, canopy cover is calculated as the percent of upper layer canopy cover within a stand, which is represented as a 30-m pixel. Once the tree crown classification has been completed, canopy cover is determined by summarizing the percentage of hyperspatial (5 cm) pixels that were classified as canopy vegetation within each 30-m Landsat pixel. Finally, the post-fire canopy cover is compared against pre-fire canopy cover from the LANDFIRE project to determine how much upper layer canopy cover reduction occurred due to the fire, giving an indication of tree mortality, which is a primary indicator of post-fire effects in forested ecosystems.

#### *2.1. Classification of Tree Crowns from Hyperspatial Imagery*

The initial step in mapping canopy cover and tree mortality enables the mapping of individual tree crowns within a hyperspatial orthomosaic acquired with the use of an sUAS. By classifying the tree crowns within an orthomosaic, we will be able to distinguish between landcover lifeforms that are indicative of canopy vegetation as opposed to surface vegetation.

The mask region-based convolutional neural network (MR-CNN) is an algorithm for instance segmentation. It detects objects in an image and determines which pixels comprise each object. When classifying an image, MR-CNN begins by applying a convolutional neural network (CNN) to extract features from the image. The features extracted from the image are used as input for a region proposal network, which slides over the feature map and detects regions of interest (RoIs) which are likely to contain objects. Each RoI is evaluated further to determine what type of object, if any, is inside (in this case, tree and non-tree); the object's bounding box; and what pixels inside the bounding box comprise the object [12].

Training data for the MR-CNN model consists of images in which the boundaries of all trees are marked with individual polygon. Unmarked areas are assumed to be non-tree. The model is trained using mini-batch gradient descent to minimize a loss function that measures the difference between the model's output in response to a training image and the training data.

A MR-CNN model with a ResNet-101 backbone that had been pre-trained on the Microsoft COCO dataset [13] was retrained to detect trees. Training images were cropped from orthomosaics of study areas flown with sUAS over the Boise National Forest. These

forested sites primarily contained Pinus Ponderosa (Ponderosa Pine), Populus Trichocarpa (Black Cottonwood), Pseudotsuga Menzieii (Douglas Fir), and Pinus Contorta (Lodgepole Pine). A polygon was traced around each tree occurring in each training image using the VGG Image Annotator [14] as shown in Figure 1. A total of 1148 trees were labeled over 14 training images. Questions of whether visually ambiguous objects were trees or surface vegetation were answered by referring to a 3D model created from the same imagery as the orthomosaics. To discourage overfitting, training images were rotated and mirrored randomly during the training process. The training data was created using many different forest types which also included images that contained scorched trees. These scorched trees were put in the "nontree" class to help the MR-CNN distinguish between burnt and non-burnt trees.

**Figure 1.** Mask Region-based Convolutional Network training data.

As an instance segmentation algorithm, MR-CNN can answer questions that a pixelbased approach cannot, such as, "How many trees are in the image?" and "What is the distribution of tree sizes?" To facilitate calculation of canopy cover, the MR-CNN's output was used to perform semantic segmentation. All pixels that were determined to be part of a tree were marked as tree crown and all other pixels were marked as surface, an example of which is shown in Figure 2. Only photosynthesizing foliage was annotated in the training data for the MR-CNN, allowing the trained MR-CNN to distinguish between live trees and scorched trees, whose foliage was no longer photosynthesizing, during classification. Scorched trees would be classified as "no tree" by the MR-CNN.

**Figure 2.** (**a**) is pre-classified image. (**b**) is an image that has been classified by the MR-CNN.

#### *2.2. Calculating Canopy Cover*

Canopy cover describes the vertical projection of the tree canopy onto a horizontal plane in forested ecosystems. In this effort, canopy cover is calculated as the percent canopy cover within a stand. Once the tree crown classification has been completed as described earlier, canopy cover is determined by summarizing the percentage of 5 cm pixels that were classified as canopy vegetation within each 30-m (30 m) Landsat pixel. Because the MR-CNN is trained to only classify live trees, scorched canopy vegetation is not included in the canopy cover calculations.

#### 2.2.1. Calculate Five Centimeter Class Densities for Each 30 Meter Pixel

The density of the tree class is calculated for each 30-m Landsat pixel from the 5 cm tree classification pixels obtained from the MR-CNN. A tool was written which calculates the density of 5 cm pixels classified as tree crown within each 30 m pixel, outputting a 30 m resolution grayscale TIFF image. This Density Tool counts the total number of hyperspatial tree class pixels occurring in the 30 m pixel. The total number of tree pixels that occurred is then divided by the total number of 5 cm pixels within the 30 m pixel and multiplied by 100. The class densities for each 30 m pixel are recorded in the class density output raster. This class density raster contains the percentage of each 30-m pixel that is comprised of 5 cm pixels that are classified as being part of a tree crown, which effectively records Canopy Cover for each 30-m pixel.

#### 2.2.2. Canopy Cover Analysis

Previous efforts including Hamilton [1,4] and Matese [15] have shown that mapping with higher resolution imagery results in increased accuracy than using lower resolution imagery. This section summarizes the analysis performed to assess the effectiveness of the methods specified previously for calculating canopy cover.

Setting and notation: Suppose that there is a 30-m pixel depicting a region of ground. The same 30 m × 30 m region of ground is viewed as a 600 × 600 region of 5 cm pixels in the hyperspatial orthomosaic. Then the canopy cover is calculated from the hyperspatial orthomosaic, and this calculation is subject to error. The probability of correctly identifying a "tree" pixel (sensitivity) is observed from training data to be at least 0.70, and the probability of correctly identifying a "non-tree" pixel (specificity) is observed from training data to be at least 0.95. Both values are usually higher (examples are shown in Section 3.2), but we will use the values of 0.70 and 0.95 as conservative estimates.

Question: What level of accuracy can one expect for the value of canopy cover calculated from the 600 × 600 pixel region of the hyperspatial orthomosaic corresponding to a Landsat pixel? Precisely, what is the standard error of this value? Is this significantly smaller than what one might hope to obtain from LANDFIRE data?

Answer: Using The value of canopy cover can be estimated with bias close to 0 and standard error less than 0.001080. In most cases, the standard error is less than 0.000824. These values of standard error are significantly smaller than what one might hope to obtain from LANDFIRE data. They are obtained using mathematics beyond just routine statistical inference. The details are omitted here so as not to disturb the flow of the article. However, they are included in the Appendix A for the interested reader.

#### *2.3. Canopy Cover: sUAS Derived vs. LANDFIRE*

The LANDFIRE project is a collaborative effort between the United States Forest Service and the US Department of the Interior, producing geospatial products that describe wildland fire fuel conditions supporting wildland fire modeling across the United States. LANDFIRE's Existing Vegetation Cover and Forest Canopy Cover products record the vertically projected percentage cover of the upper vegetation canopy layer [8]. For this effort in evaluating tree mortality, the Existing Vegetation Cover raster was updated to only contain canopy cover in forested ecosystems. Additionally, the Existing Vegetation Cover was updated (via the ArcGIS Spatial Analyst Con tool) with the Forest Canopy Cover for forested pixels where the predominant lifeform appeared to be either shrub or herbaceous. This conversion was accomplished by reclassing the vegetation cover for all non-forested lifeforms to zero, then substituting the Forest Canopy Cover for the Existing Vegetation Cover on pixels that had a non-forested life form, but also had a Forest Canopy Cover value. This resulted in a geospatial layer with spatial resolution of 30 m representing vertically projected percentage forest canopy cover.

The comparison of the hyperspatial derived canopy cover (HDCC) as accomplished using the ArcGIS combine raster tool. This resulted in an attribute table that had each combination of canopy cover that occurred within pixels in the two rasters as well as how many pixels the combination occurred in. From the combine table, the difference between the HDCC and LANDFIRE forest Canopy Cover (LFCC) can be calculated. This difference can be used to analyze the statistical difference between the sUAS derived and LANDFIRE products as well as the determination of canopy reduction as a measure of tree mortality in areas where the sUAS canopy cover was calculated from imagery acquired after a wildland fire.

#### 2.3.1. Canopy Cover Comparison Methods

When classifying sUAS imagery through the Mask Regional-Convolutional Neural Network (MR-CNN), a new image is created as an output image in which pixels are based on the original image is part of a tree crown (tree) or does not comprise part of a tree crown (no tree). During this process, however, spatial referencing and null data (where a pixel represents no data) values are lost and must be added back.

Spatial referencing is added back using an ArcPy script which takes the world file spatial referencing from the original images that were aligned with the 30 m imagery, inputted into the MR-CNN and adds the spatial referencing to the MR-CNN output.

Band 4 of the classified output from the MR-CNN contains the tree/no tree classification. This band is extracted using the Copy Raster Tool on Band 4 of the image. This newly extracted classification raster is given as input to the density tool, sums the number of 5 cm tree pixels within the 30 m pixel and divides it by the total number of 5 cm pixels within a 30 m pixel. Multiplying this number by one hundred gives the percent of tree pixels or the canopy cover in the 30 m pixel area as expressed in Equation (1). This process is repeated until each pixel in the image has been analyzed giving a 30 m resolution canopy cover (CC) output raster.

$$^{90}\text{\% of CC} = \left(\frac{\text{tree pixels}}{\text{total pixels}}\right) \ast 100\tag{1}$$

Before adding the null data back into the orthomosaic, the 5-cm resolution imagery is resampled to 30 m resolution. Resampling before adding the null data speeds up the process of adding the null data back in. Null data is recorded with a pixel value of 0 (not a tree) by the MR-CNN. To restore its value to NULL, ArcPy's conditional toolset is used to compare pixels in the source image (the MR-CNN output) and with the original data to check if the value was NULL. The density tools output is changed to 255 for every pixel where the source image's value is NULL. To represent NULL, 255 is used because it is outside the range of possible values (0% to 100%).

The various methods of processing the data outlined in the previous sections creates a HDCC raster that is aligned with the LFCC raster. This allows the rasters to be compared pixel by pixel, using the ArcGIS Combine Raster tool. This creates a new raster with a raster attribute table column for each of the combined rasters. Next, the LFCC column was subtracted from the HDCC column into a newly added Difference column. With sUAS data of unburned areas, this shows how close HDCC and LFCC data are. When comparing pre-fire LFCC against post-fire HDCC, this difference can be used in the calculation of canopy reduction, an indicator of tree mortality which is a primary effect of wildland fire.

The workflow which enable the calculation of the HDCC and its camparison to the LFCC is summarized in Figure 3.

**Figure 3.** An overview of the processes used to create the HDCC.

#### 2.3.2. Canopy Cover Comparison Analysis

Working assumption: In Section 2.2.2, it was observed that (assuming 5-cm pixel identification better than 0.80 correct for pixels representing ground with canopy cover, better than 0.95 for pixels representing ground without canopy cover, and no more than 0.80 canopy cover) using adjusted hyperspatial derived data, canopy cover can be calculated with a standard error of less than 0.000824. This is a much smaller standard error than we could hope to obtain from 30 m LANDFIRE pixels. Accordingly, it is reasonable, when

comparing LANDFIRE data to hyperspatial derive data, to make the working assumption that the adjusted HDCC value is correct.

Objective: The method of maximum-likelihood estimation is used to construct a model of the error of the value of canopy cover derived from LANDFIRE data.

Setting and notation: There are *n* LANDFIRE 30 m pixels, and for each, the estimated value *B* ∈ [0,1] of canopy cover. For each there is a corresponding region of 5 cm pixels in the hyperspatial orthomosaic, and the adjusted value of canopy cover *C* ∈ [0,1] calculated from the 600<sup>2</sup> pixels in corresponding region of the hyperspatial orthomosaic. As stated above, the latter value (*C*) is presumed correct, and the objective is to specify a model of the error in the former value (*B*). Write

$$X = B - \mathbb{C} \tag{2}$$

for the error. The data consists of n pairs (*b*1, *c*1), (*b*2, *c*2), ... , (*bn*, *cn*), from which are obtained n error values

*x*<sup>1</sup> = *b*<sup>1</sup> − *c*1,... , *xn* = *bn* − *cn*. (3)

Model for error: A normal model for *X* is used, whose density is

$$fX(\mathbf{x}) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \tag{4}$$

with parameters *μ* (mean) and *σ* (standard deviation). (Note that *μ* and *σ* are the respective mean and standard deviation of the error *B* − *C* only, and not of the canopy cover itself or any estimate thereof). The values of the parameters are unknown, and the next task is to choose the best values for them (so that the model best accounts for the present data).

Maximum-likelihood estimators for parameters: The maximum-likelihood estimations for *μ* and *σ* are

$$
\hat{\mu} = \overline{X} = \frac{\sum X}{n} \tag{5}
$$

and

$$\mathcal{P} = \mathcal{S} = \sqrt{\frac{\sum \left( X - \overline{X} \right)^2}{n}} \tag{6}$$

respectively, the latter being the "population standard deviation" *S*ˆ with denominator *n* instead of the "sample standard deviation" *S* with denominator *n* − 1.

Conclusion: The maximum-likelihood estimation for a normal error *X* of the LAND-FIRE value for canopy cover has mean *X* and standard deviation *S*ˆ, calculated from the *n* values *x*<sup>1</sup> = *b*<sup>1</sup> – *c*1, ... , *xn* = *bn* − *cn* of the difference *X* = *B* − *C* between LANDFIRE Canopy Cover and (adjusted) Hyperspatial Derived Canopy Cover.

#### *2.4. Mapping Tree Mortality via Canopy Reduction*

While it is much more accurate to map canopy cover from hyperspatial sUAS imagery, the mapping extent is limited to how large of a study area can be flown with a sUAS. Current regulatory and technical constraints limit flight operations with a single sUAS in forested environments to under 400 hectares in a single day [16]. For mapping post-fire canopy cover, this constraint can be overcome with more flights on all but the largest fires. Due to the unpredictable nature of fires resulting from unplanned ignitions, it is not feasible to obtain hyperspatial sUAS derived canopy cover data reflecting pre-fire conditions. Fortunately, the LANDFIRE products are generated for the entire United States and are updated every two years.

For mapping canopy reduction, which is a proxy for tree mortality, a comparison was made between the LFCC (pre-fire canopy cover) against the HDCC (post-fire) layer, as described in the previous section.

Once the difference between pre-fire LFCC and post-fire HDCC, canopy cover is calculated, canopy cover reduction (CCR) can be calculated by calculating the difference between LFCC and HDCC, minus the estimated error standard deviation (multiplied by the Standard Normal density) and estimated error mean as shown in Equation (7).

#### 2.4.1. Canopy Reduction Methods

Due to LANDFIRE's canopy cover being less accurate than the HDCC, it cannot simply be stated that reduction is the discrepancy between sUAS Derived data and LANDFIRE data. Instead, we must use the data from unburned areas to see how much this discrepancy varies. Then any numbers outside of the normal variance is a change in canopy cover to a 95% degree of certainty. Once it is determined that there has been a change, there are two possibilities for what that change could be, active crown fire or passive crown fire. If the HDCC is 0% it was an active crown fire in this area otherwise it was a passive crown fire.

The main difficulty with calculating the type of fire based on this method is that there is some uncertainty due to LANDFIRE being 30-m data. For this reason, it is difficult to tell what happened to pixels where there was not enough of a difference to be certain if there was a change. The discrepancy between the data could just be due to LANDFIRE's overreporting or it could be due to a passive crown fire.

#### 2.4.2. Canopy Reduction Analysis

Working assumption: As before, a 30 m-by-30 m square of land is being considered, which corresponds to a single LANDFIRE pixel as well as to a 600 × 600 pixel region in a hyperspatial orthomosaic. Again, *B* ∈ [0,1] is the canopy cover according to LANDFIRE, and *C* ∈ [0,1] is the adjusted value of canopy cover obtained from the hyperspatial orthomosaic, which is assumed to be correct. In the previous section, maximum-likelihood estimation was used to construct a model of the error *X* = *B* − *C*, which was a normal density with estimated parameter values *μ* = *μ*ˆ = *X* and *σ* = *σ*ˆ = *S*ˆ.

Objective: A criterion is developed by which it can be concluded, with the desired level of confidence 1 − *α*, that canopy cover, from past LANDFIRE image to present hyperspatial orthomosaic, has decreased by at least some particular value *r*.

The difficulty: The difficulty here is that there is no available pre-fire value of *C* (canopy cover calculated from a hyperspatial orthomosaic), but only a past value of *B* (canopy cover calculated from LANDFIRE) to compare with the present value of *C*. The notation *X* = *B* − *C*, is still used, while now it is recognized that *B* is a past value and *C* is a present value. This means that *X* (difference in LFCC and (adjusted) HDCC values) may include change in canopy cover (where a positive value might result from a crown fire), or inaccuracy of LANDFIRE data, or both.

The null hypothesis: The null hypothesis H0 that we wish to test is that there was a reduction in canopy cover of no more than *r*. This H0 will be rejected if, even allowing for a canopy reduction of *r*, the value of *B* − *C* is too high to be explained by chance; that is, if

$$B - C > r + \mu + z\_n \sigma,\tag{7}$$

or, equivalently,

$$X \succ r + \mu + z\_{\mathfrak{a}} \sigma. \tag{8}$$

Replacing the inequality with an equality and solving for *r*, the following criterion is obtained.

Criterion: It may be inferred that canopy cover has decreased by at least

$$B - C - \mu - z\_a \sigma\_\prime \tag{9}$$

where *B* is past canopy cover according to LANDFIRE, *C* is the adjusted value of canopy cover from a hyperspatial orthomosaic, *μ* and *σ* are the parameter values obtained as described previously, and *zα* is the critical value of the Standard Normal density at significance level *α*. Intuitively, the amount that we can conclude by which canopy cover has been reduced is the difference "past LANDFIRE minus present hyperspatial" *B* − *C*, minus the average discrepancy *μ* between the LANDFIRE *B* and hyperspatial-derived *C*, minus

an additional margin of error *zασ* which is larger, the higher the level of confidence we wish to have.

#### **3. Results**

#### *3.1. Tree Crown Classification Validation*

To assess the performance of the MR-CNN in mapping tree crowns, validation data was prepared using ArcGIS. Sample areas of canopy and surface vegetation were recorded as polygons in shapefiles digitized from existing orthomosaics. Validation data were made for all orthomosaics used in the research.

The MR-CNN was run on the orthomosaics used for validation. Agreement between a classifier's output and the validation data was calculated using ArcMap's Tabulate Area tool as shown in Table 1. A classifier's accuracy on an orthomosaic was computed as the area where the classifier agreed with the class marked in the validation data divided by the total area marked in the validation data.

**Table 1.** Confusion matrices of MR-CNN tree crown classifications. Classified surface and canopy are shown in the rows while human labeled surface and canopy are in columns.


Orthomosaics where choose with a variety of different conditions to find any anomalies where the classifier underperforms. Trees in the chosen orthomosaics were comprised primarily of Pinus Ponderosa (Ponderosa Pine), Populus Trichocarpa (Black Cottonwood), Pseudotsuga Menzieii (Douglas Fir) and Pinus Contorta (Lodgepole Pine). All orthomosaics were from flights conducted with the Boise National Forest in southwestern Idaho [16,17]. The distinguishing characteristics of the ten chosen orthomosaics were as follows:


The MR-CNN model's overall accuracy was 90%. The specificity (correctly identified tree pixels/total number of tree pixels) was 84% and the sensitivity (correctly identified non-tree pixels/total number of non-tree pixels) was 99%. The MR-CNN was found to correctly identify a high number of tree pixels while only missing a low number of tree pixels as shown in Table 2.


**Table 2.** MR-CNN tree crown classification statistical summary.

#### *3.2. Canopy Cover Analysis*

Running an orthomosaic through the Snapper Tool, MR-CNN, and Density Tool gave a 30 m orthomosaic with canopy cover that also is coregistered with LANDFIRE data. Then the HDCC numbers were adjusted based on the bias of the MR-CNN. The values used to calculate the bias were the accuracy of correctly identifying tree pixels (sensitivity) and the odds of a false positive (specificity). These numbers ended up being 84% and 99%, respectively. Figure 4 shows the HDCC layer for South Placerville.

**Figure 4.** Values are in percent canopy cover.

#### *3.3. Canopy Cover Comparison: Hyperspatial Derived vs. LANDFIRE*

Ten orthomosaics where there was no fire were compared using the methods stated above. LANDFIRE was over reporting by an average of 2.5% over the HDCC layer. This agrees with what Scott [18] stated, that the "canopy cover values [for LANDFIRE] are too high." The standard deviation for the difference between each pixel was 14. This number is a little high and will hopefully be driven down by future improvements. Figure 5 shows a piece of the results for Edna Creek. White represents a higher canopy cover for the HDCC layer while white represents a higher canopy cover reported by LANDFIRE. The grey composing most of the figure is for when the difference is equal or close to zero.

**Figure 5.** A raster representing the comparison between the LANDFIRE and HDCC data. White means LANDFIRE reported more canopy cover while black indicates the inverse.

Table 3 shows that the LFCC over reported for every orthomosaic.


**Table 3.** Mean difference and the standard deviation of that difference for each area.

#### *3.4. Tree Mortality Mapping Assessment*

The standard deviation and mean collected from the ten non-burnt orthomosaics show that there needs to be at least a change of 27% to know there has been a change with 95% confidence. If the difference is high enough to say there has been a change next the HDCC layer can be looked at to see whether it was a passive crown fire or an active crown fire. Pixels that did not contain adequate canopy reduction (more than 26%) from comparing pre-fire LFCC to post-fire HDCC are considered to be inconclusive. This would include

unburned areas, surface fire, passive crown fire, and active crown fire. If an area had less than 27% canopy cover then we cannot say whether any canopy reduction occurred. This method was used on orthomosaics acquired over portions of the Mesa fire (a 16,000-hectare mixed severity fire located on the Payette National Forest in southern Idaho) to calculate burn intensity. This fire worked well as there were areas of active crown fire, passive crown fire, and surface fire. Figure 6 shows the Southwest area of the Mesa fire which had a lot of Active Crown fire, while Figure 7 the east section which had more passive crown fire and possible surface fire.

**Figure 6.** (**a**) Hyperspatial data after a burn. (**b**) Types of fire: Red = Active, Yellow = Passive, Black = inconclusive crown fire activity.

**Figure 7.** (**a**) Hyperspatial data after a burn. (**b**) Types of fire: Red = Active, Yellow = Passive, Black = inconclusive crown fire activity.

#### **4. Discussion**

Post-fire mapping products are commonly used by land managers as they are developing post-fire management plans which specify strategies for mitigating the effects of the fire. In order for this data to be most useful, it needs to be readily available to managers. The advantage of using imagery acquired by sUAS is that it can be available as quickly as

managers can get out to the site and conduct image acquisition flights. The timing of the flights can be tailored to optimize meteorological conditions that are conducive to acquiring image acquisition by enabling managers to acquire imagery as soon as the scene is not obscured by smoke and there is low enough cloud cover to ensure adequate illumination of the scene.

The MR-CNN was found to be very effective in identifying tree crowns. The classifier was found to have virtually no false positives with a specificity of 99%. The false negatives were still minimal, with a sensitivity of 84%.

Scott [18] showed that LANDFIRE Canopy Cover typically is over reported. As expected, when comparisons were made of study areas that had not experienced wildland fire, the HDCC was on average 2.5 percentage points lower than the LANDFIRE Canopy Cover. While the mean difference between the HDCC and LFCC was only a few percentage points, the higher standard deviation which typically ranged a bit more than 10 percentage points indicates that the canopy cover between the HDCC and LFCC products differed more significantly than indicated solely by the means.

The canopy reduction mapping over a burned area was very promising. When ran on the Mesa fire most pixels correctly identified surface fire as well as both active and passive crown fire. These methods work best on areas that had higher pre-fire canopy cover and canopy reduction. There needs to be at least a 27% difference for there to be a known change. Various theories to bring this number down are brought up in the future work Section 5.1.

In addition to higher accuracy mapping of canopy cover as well as tree mortality as measured by canopy reduction, the use of sUAS allows for much more current data acquisition than is possible with the LANDFIRE products. Using imagery acquired with a sUAS, it would be possible to have the HDCC product available within a day of flying a study area. By contrast, the current LANDFIRE 2.0.0 Remap effort is using Landsat imagery that was acquired in 2016. The initial product releases which covered the majority of the Contiguous US did not occur until 2019, meaning even the initial product releases were derived from imagery that was already three years out of date [8]. While canopy cover can be mapped with much higher accuracy using hyperspatial imagery acquired with a sUAS, this research team found that current technical and regulatory constraints on drone usage realistically will only allow for the acquisition of up to 600 hectares per day [16]. To effectively map large class F (>400 ha) and G (>5000 ha) fires [19], larger and more efficiently acquirable imagery satellite system still needs to be utilized to obtain a large enough analysis extent to map the whole fire. As the spatial resolution of satellite imagery increases with some of the new high-resolution sensors that are being put into orbit, it may be possible to leverage this methodology, applying it to new generations of imagery as they approach hyperspatial resolution.

#### **5. Conclusions**

Detection and mapping of trees from hyperspatial drone imagery enables more accurate mapping of canopy cover, a measure of forest density which is used to determine ecosystem departure from historic fire regimes, a primary measure of ecosystem resilience [20]. A comparison of a pre-fire LANDFIRE Canopy Cover product against post-fire HDCC resulted in a very effective measure of canopy fire reduction, which is a measure of tree mortality resulting as an effect of the wildland fire.

#### *5.1. Future Work*

While this effort was highly effective in establishing a timely, effective and costefficient means of mapping canopy cover and tree mortality, there are many areas of future research which could result in an improvement in the effectiveness of this methodology.

An improvement in wildland fire extent mapping can be achieved by identifying unburned trees that are surrounded by surface vegetation. Current classification methods developed through this effort [1,3] which map burn extent are not capable of detecting that

a surface fire burned underneath an unburned tree crown due to the view of the surface vegetation being obscured from the sensor on the sUAS by the tree crown. While it may not be possible to determine whether the biomass consumption under the tree crown is high or low due to the presence of white ash, it would be possible to identify that the surface vegetation was burned and update the burn extent and aerial mapping accordingly based on assumptions made based on whether surface vegetation visible next to the crown was burned.

While the specificity of the MR-CNN was quite high at 99%, the accuracy of the classifier could be improved by refining the training data for the MR-CNN, based on observations of false negatives during validation. This refinement would likely improve tree crown classification, reducing the number of false negatives, thereby increasing the classifier sensitivity, which was 84% in the validation dataset used for this effort.

An improvement in tree crown classification accuracy might be able to be achieved by using pixels' heights as additional input to the classifiers. Per-pixel height values were available from a digital surface model produced during the orthomosaic creation process but were not used as an input to the MR-CNN in this effort.

Due to the need to split each orthomosaic, the MR-CNN classifier can miss trees that are on the border of these splits since they look like half trees. This problem could be eliminated by modifying the snapping and density tool to run larger files than what they are currently limited to due to OpenCV's file size limits. This would allow entire orthomosaics to be processed at once.

A more major problem that arises with splitting up the image is when 30 m pixels overlap. This happens because the LANDFIRE and hyperspatial image cannot line up perfectly. Because of this the 5 cm data has an extra pixel on the bottom and right when it is clipped. The density tool creates a whole 30-m pixel from these few extra 5 cm pixels on the edges. When the 30-m rasters are recombined the data at the overlap is corrupted. This problem can be fixed either creating a more direct method to handle the extra pixel or by getting rid of the split. Getting rid of the split would require OpenCV to be able to handle larger orthomosaics. This problem is most likely a big factor in the high standard deviation. If fixed, the difference needed to be confident there has been a change should drop.

**Author Contributions:** Conceptualization, D.A.H.; Methodology, D.A.H., K. Brothers, S.D.J., J.W.; software, D.A.H., K.L.B., S.D.J., J.W.; validation, K.L.B.; formal analysis, J.C.; investigation, D.A.H., K.L.B., S.D.J.; data curation, D.A.H.; writing, D.A.H., J.C., K.L.B., S.D.J., J.W.; visualization, K.L.B.; supervision, D.A.H.; project administration, D.A.H.; funding acquisition, D.A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the USDA Forest Service Boise National Forest, Forest Service Agreement No, 18-CS-11040200-0025.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data and source code from this project is available at the https: //firemap.nnu.edu/tree-mortality.

**Acknowledgments:** We would like to acknowledge the students in the Northwest Nazarene University Department of Math and Computer Science who have helped with different aspects of this effort, including current and former students Nicholas Hamilton, Jonathan Hamilton, Braelyn Boerner, Gabriel Murphy, Enoch Levandovsky, Gabriel Johnson, Austin White, Aleesha Chavez and Andrew Welk. Additionally, we would like to acknowledge Northwest Nazarene University who funded the research efforts of many of the students mentioned.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Accuracy of Calculated Value for Canopy Cover**

In this appendix, the accuracy of the value for canopy cover calculated from the hyperspatial orthomosaic is shown to be better than that calculated from LANDFIRE. This conclusion was stated in Section 2.2.2, but depends on mathematics beyond just routine statistical inference. The details are given here to justify the conclusion.

Setting and notation: Suppose that there is a 30 m pixel depicting a region of ground with canopy cover value *c* ∈ [0,1]. The actual value *c* would be unknown in practice. The same 30 m × 30 m region of ground is viewed as an *n* × *n* region of smaller pixels in the hyperspatial orthomosaic. In practice, *n* would be equal to 600, and the 600 × 600 orthomosaic would consist of 5 cm pixels, but the symbol *n* will be used for now for greater clarity.) The number of true "tree" pixels is *n*2*c*, and the number of "non-tree" pixels is *<sup>n</sup>*2(<sup>1</sup> − *<sup>c</sup>*). Each of the *<sup>n</sup>*<sup>2</sup> pixels is assigned by the MR-CNN a binary value *<sup>A</sup>* = 1 (tree) or *A* = 0 (non-tree). Then the canopy cover is calculated from the hyperspatial orthomosaic as

$$
\overline{A} = \frac{\sum A}{n^2},
\tag{A1}
$$

where the sum is taken over all *n*<sup>2</sup> pixels.

The assignments are subject to error. Suppose the probability of correctly assigning *A* = 1 to a "tree" pixel (sensitivity) is known to be *p* ∈ [0,1], and the probability of correctly assigning *A* = 0 to a "non-tree" pixel (specificity) is known to be *q* ∈ [0,1].

Question: What level of accuracy can one expect for the value of canopy cover calculated from the *n* × *n* pixel region of the hyperspatial orthomosaic corresponding to a Landsat pixel? Precisely, what is the standard error of this value? Is this significantly smaller than what one might hope to obtain from LANDFIRE data?

Mean and standard deviation of calculated value: The expected value (mean) of the calculated value *A* is

$$E[\overleftarrow{A}] = cp + (1 - c)(1 - q) \ = 1 - q + c(p + q - 1). \tag{A2}$$

The standard deviation is

$$\sqrt{\frac{n^2cp(1-p)+n^2(1-c)q(1-q)}{(n^2)^2}} = \frac{\sqrt{cp(1-p)+(1-c)q(1-q)}}{n} \tag{A3}$$

Adjustment of value for canopy cover: Recall that *c* is the proportion of small pixels representing ground that has canopy cover. It is the number wanted but not known. One solves for *c* in the expression for *E A* above to obtain

$$\varepsilon = \frac{E\left[\overline{A}\right] + q - 1}{p + q - 1},\tag{A4}$$

and thus

$$\mathcal{C} = \frac{\overline{A} + q - 1}{p + q - 1} \tag{A5}$$

is an unbiased estimator of *c*. This will be called the "adjusted canopy cover". The standard error of *C* is just the standard deviation of *A* divided by *p* + *q* − 1, which is

$$\frac{\sqrt{cp(1-p) + (1-c)q(1-q)}}{n(p+q-1)}\tag{A6}$$

Next, this expression will be examined more closely to obtain an upper bound. Derivation of upper bound: The quantity *x*(1 − *x*) is decreasing for on the interval [0.5, 1].

So that since *p* > 0.70 and *q* > 0.95,

$$p(1-p) < (0.70)(1-0.70) \text{ and } q(1-q) < (0.95)(1-0.95). \tag{A7}$$

It also follows that

$$(0.70)(1 - 0.70) \text{ > (0.95)} (1 - 0.95),\tag{A8}$$

whereupon

$$c(0.70)(1 - 0.70) + (1 - c)(0.95)(1 - 0.95) \tag{A9}$$

(which is a weighted mean of (0.70)(1 − 0.70) and (0.95)(1 − 0.95)) is an increasing function of *c*. Combining this fact with the first two inequalities, one obtains

$$\begin{array}{c} c p(1-p) + (1-c)q(1-q) \\ < c(0.70)(1 - 0.70) + (1-c)(0.95)(1 - 0.95) \\ < (0.80)(0.70)(1 - 0.70) + (1 - 0.80)(0.95)(1 - 0.95) \end{array} \tag{A10}$$

for

$$p > 0.70, \quad q > 0.95,\qquad\text{and } c < 0.80. \tag{A11}$$

The actual canopy cover (*c*) in the training data used was less than 0.80, the sensitivity (*p*) was better than 0.70, and the specificity (*q*) better than 0.95. Using also *p* + *q* − 1 > 0.70 + 0.95 − 1*,* one obtains

$$\frac{\sqrt{cp(1-p) + (1-c)q(1-q)}}{n(p+q-1)} < \frac{\sqrt{(0.80)(0.70)(1-0.70) + (1-0.80)0.95(1-0.95)}}{600(0.70 + 0.95 - 1)} \approx 0.001080,\tag{A12}$$

which compares favorably with results obtained with LANDFIRE 30 m pixels. If one assumes that *p* > 0.80 (which is true in most cases) rather than merely *p* > 0.70, the upper bound on the standard error of *C* becomes

$$\frac{\sqrt{(0.80)(0.80)(1 - 0.80) + (1 - 0.80)0.95(1 - 0.95)}}{600(0.80 + 0.95 - 1)} \approx 0.000824,\tag{A13}$$

which would be impossibly low to achieve with LANDFIRE 30 m pixels.

Conclusion: The adjusted value of canopy cover calculated from the hyperspatial orthomosaic has lower standard error than we could hope to obtain from LANDFIRE 30-m pixels.

#### **References**

