Using Landsat Spectral Indices in Time-Series to Assess Wildfire Disturbance and Recovery

Satellite earth observation is being increasingly used to monitor forests across the world. Freely available Landsat data stretching back four decades, coupled with advances in computer processing capabilities, has enabled new time-series techniques for analyzing forest change. Typically, these methods track individual pixel values over time, through the use of various spectral indices. This study examines the utility of eight spectral indices for characterizing fire disturbance and recovery in sclerophyll forests, in order to determine their relative merits in the context of Landsat time-series. Although existing research into Landsat indices is comprehensive, this study presents a new approach, by comparing the distributions of pre and post-fire pixels using Glass’s delta, for evaluating indices without the need of detailed field information. Our results show that in the sclerophyll forests of southeast Australia, common indices, such as the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR), both accurately capture wildfire disturbance in a pixel-based time-series approach, especially if images from soon after the disturbance are available. However, for tracking forest regrowth and recovery, indices, such as NDVI, which typically capture chlorophyll concentration or canopy ‘greenness’, are not as reliable, with values returning to pre-fire levels in 3–5 years. In comparison, indices that are more sensitive to forest moisture and structure, such as NBR, indicate much longer (8–10 years) recovery timeframes. This finding is consistent with studies that were conducted in other forest types. We also demonstrate that additional information regarding forest condition, particularly in relation to recovery, can be extracted from less well known indices, such as NBR2, as well as textural indices incorporating spatial variance. With Landsat time-series gaining in popularity in recent years, it is critical to understand the advantages and limitations of the various indices that these methods rely on.


Introduction
As scientists shift toward viewing the earth as a single interconnected system [1], understanding forest dynamics and the complex relationships with human societies becomes more and more pertinent.This type of analysis requires a multidisciplinary approach integrating local, regional, and global knowledge [2].Remote sensing via satellite is ideally suited to meet these needs, especially when considering the large and often remote areas that forests occupy.Among available satellites, the Landsat program offers an unparalleled historical record stretching back over four decades.The opening of the image archive in 2008, along with the advances in computer processing, has led to a plethora of new and novel applications exploiting Landsat time-series [3].Commonly, in the forest domain, these studies look to establish disturbance and recovery histories, following events such as wildfire, logging, and insect damage [4][5][6].Using a time-series, rather than image pairs, allows for change to be differentiated from background noise, whilst also capturing longer-term ecological trends [4].
Methods for characterizing forest dynamics (abrupt changes and longer term trends) using time-series differ, but a point of similarity is the use of spectral indices.Spectral indices convert multi-spectral satellite data into a single component, so individual pixels can be tracked through time.Spectral indices also have an advantage over single bands by amplifying desired effects (e.g., changes in vegetation condition) and reducing unwanted features, such as atmospheric and topographic noise [7].There are numerous spectral indices in the literature.However, when considering those commonly used in Landsat derived pixel-based time-series, the field narrows significantly.Frequently used is the Normalized Difference Vegetation Index (NDVI) [8].NDVI is a measure of photosynthetic biomass, and has been shown to correlate well with ecological parameters, such as the fraction of green vegetation cover [9] and leaf area index [10].NDVI is sensitive to changes in vegetation condition, and has been shown to accurately detect forest disturbances.However, it is generally considered to be less adept in representing forest recovery, due to grasses and other non-woody vegetation colonizing a site after a disturbance, and consequently returning the NDVI signal to its pre-disturbance state [11].In areas of sparse vegetation, NDVI can also be adversely affected by soil reflectance.To correct for soil brightness, Huete developed the Soil Adjusted Vegetation Index (SAVI), which incorporates a soil correction factor into the NDVI formula [12].
Indices using short-wave infrared (SWIR) bands are commonly used in Landsat time-series, as these wavelengths are sensitive to forest structure, moisture, shadowing, and vegetation density [5].The Normalized Burn Ratio (NBR) is a ratio of the near-infrared and second SWIR band (2.08-2.35µm), and was developed by Key and Benson [13] to identify burned areas following fire and provide a quantitative measure of burn severity.Several authors have found NBR to correlate highly with field-based measurements in forest ecosystems [14][15][16], however Roy et al. [17] suggest caution when using NBR for burn severity mapping as their investigations indicated sub-optimal results.In Landsat time-series NBR is used extensively, and has proven adept at characterizing forest dynamics in the United States of America (USA) [18] and Canada [19].Similar to NBR is the Normalized Difference Moisture Index (NDMI), which uses the near-infrared with the first SWIR band (1.55-1.75µm).NDMI is sometimes favored for tracking disturbances other than fire, and was used by Goodwin et al. [20] for classifying areas that were disturbed by the Mountain Pine Beetle in western Canada.NBR2 is another variation of a ratio/difference index, contrasting the two Landsat SWIR bands.It is provided as a standard product by the United States Geological Survey (USGS), but is rarely used in the literature.Storey et al. [21] found it useful for post-fire recovery assessment in chamise chaparral vegetation in southern California, while Stroppiana et al. [22] used it as part of an ensemble to map burned areas.
The Tasseled Cap (TC) transformation of Landsat Multispectral Scanner (MSS) data was first presented by Kauth and Thomas in 1976, and was later adapted by Crist and Cicone for Landsat TM data [23].The various components of TC are created via linear transformations using defined coefficients.In simplified terms, Brightness (TCB) represents the overall brightness of all bands, Greenness (TCG) is a contrast between the visible and near-infrared bands, and Wetness (TCW) is a contrast of the visible and near-infrared with the SWIR bands, making it sensitive to soil and plant moisture [23].TC Angle (TCA) [24] is calculated as the arctan of TCG/TCB and describes the vegetation cover within the TCB-TCG spectral plane [25].Various time-series studies have shown success with TC components.For instance, Senf et al. [6] used TC components to track insect disturbance in British Columbia, Canada, and found that TCG was useful for detecting Western Spruce Budworm disturbance, whereas TCW and TCB were better indicators of Mountain Pine Beetle disturbance.
The time-series method, and the index (or indices) used, can significantly alter the outcomes of a study, as highlighted recently by Cohen et al. [26].Often, studies evaluating spectral indices look to establish the strength of the relationship between the index and field data [16].An alternative approach, especially when field data are not available, is to use human interpreted reference data as a validation method.In a recent study, Schultz et al. [27] assessed eight spectral indices in their ability to detect deforestation in the tropics, using manually interpreted reference pixels for training and validation.The challenge with using field data and human interpreted imagery to train or validate models is that the data needs to be both spatially representative of the study area and temporally relevant (i.e., collected at appropriate time intervals).Many Landsat time-series studies are retrospective investigations covering large areas, and field data does not exist.Where ancillary data is available, it is more likely to indicate forest disturbance (e.g., maps of fire extent and severity) than forest recovery, which requires multiple data collections over many years.One of the strengths of satellites like Landsat are the consistent re-visit cycles, which provides the data necessary for time-series analysis.
Although research into Landsat spectral indices is comprehensive, this study adds several novel insights to the existing body of literature.Firstly, it presents a simple and robust method for assessing and comparing indices using Glass's delta, which is suitable where limited or no field data are available.Secondly, it looks at how various indices respond to fire disturbance and recovery in sclerophyll forests, which are the dominant forest type in Australia, but are also common elsewhere in the world.Thirdly, it assesses indices in the context of Landsat time-series, but independently of a specific algorithm.

Study Area
The study area contains over three million hectares of public forest in the eastern half of Victoria, Australia (Figure 1).This area was chosen because it has high ecological and economic importance, and it recently experienced three major wildfire events in the space of six years (see Figure 1 for extent of burned area).The area consists primarily of sclerophyll forests, tending to be wet in some areas and dry in others.At the wetter end, trees can attain heights over 75 m, while at the dryer end, trees are typically shorter than 40 m [28].Forest classes can be further subdivided according to Australia's National Vegetation Information System, as outlined in Table 1 (refer to Mellor and Haywood [29] for further details).The burned area is primarily contained within three major bioregions-the Victorian Alps and the Northern and Southern Highlands [30].The Alps have mild summers and cool winters, reach elevations up to 2000 m, and typically experience over 1400 mm of annual precipitation.The Highlands are located on both the northern and southern sides of the Alps, at elevations between 200 m and 1300 m, and typically experience annual rainfall between 500 and 1200 mm [28].

Tree Height (m) Canopy Cover (%)
Low (<10) Woodland (<50) Medium  Open (50-80) Tall (+30) Closed (>80) In 2003, wildfires in the northeast of Victoria burned over 1.3 million hectares of forest.Three years later, over the summer of 2006-2007, major wildfires again burned a further 1 million hectares of forest, mostly southwest of the 2003 fires.In February 2009, the devastating 'Black Saturday' fires burned 400,000 hectares across the state of Victoria [31], much of it in the Highlands region.Figure 1 indicates the extent of the burned area, which forms the study area for this research.

Landsat Data and Pre-Processing
We obtained all available Landsat TM and ETM+ surface reflectance products with less than 70% cloud-cover from 1 January to 31 March (representing southern hemisphere summer) for years 1992-2016 (paths 91/92 and rows 85/86), from the USGS archive.Surface reflectance products were processed using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [32], and include a cloud mask that was calculated with the FMask algorithm [33].We created annual summer composites using a Best Available Pixel (BAP) method of image compositing, which has been used by other studies for preparing Landsat data for use in long time-series [4,34,35].Typically, it involves choosing the first clear pixel from an image stack that is closest to a preferred day of the year, in order to minimize the effects of phenology and variations in sun angle.We chose an anniversary date of February 15 and seasonal window of plus/minus 45 days.A late summer date was used to capture fires in the year they occurred.A slight penalty (five days) was applied to ETM+ images with Scan Line Corrector errors (SLC-off), so that preference was given to TM data if available.This resulted in a time-series stack of 25 years with over 98% coverage.

Establishment of Candidate Reference Pixels
Fire maps that were maintained by the state of Victoria's land management agency [36] were used to indicate the general extent of the three large fires.Candidate reference pixels were chosen via a systematic sampling process based on the Victorian Forest Monitoring Program (VFMP) plot network [37].The VFMP plot network consists of 786 2km by 2km plots that are distributed throughout public land in Victoria, stratified by bioregion and land tenure.In each plot, 10 random pixels were selected (resulting in 7860 pixels), and a team of six worked to manually interpret each pixel to establish its disturbance history (Figure 1 shows an example of the reference pixel sampling method).This was achieved by interrogating multiple lines of evidence, such as state fire records and high resolution imagery from Google Earth.Quality assurance was performed by an independent operator, who assessed 10% of all the pixels to evaluate the accuracy of the dataset (for details see

Landsat Data and Pre-Processing
We obtained all available Landsat TM and ETM+ surface reflectance products with less than 70% cloud-cover from 1 January to 31 March (representing southern hemisphere summer) for years 1992-2016 (paths 91/92 and rows 85/86), from the USGS archive.Surface reflectance products were processed using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [32], and include a cloud mask that was calculated with the FMask algorithm [33].We created annual summer composites using a Best Available Pixel (BAP) method of image compositing, which has been used by other studies for preparing Landsat data for use in long time-series [4,34,35].Typically, it involves choosing the first clear pixel from an image stack that is closest to a preferred day of the year, in order to minimize the effects of phenology and variations in sun angle.We chose an anniversary date of February 15 and seasonal window of plus/minus 45 days.A late summer date was used to capture fires in the year they occurred.A slight penalty (five days) was applied to ETM+ images with Scan Line Corrector errors (SLC-off), so that preference was given to TM data if available.This resulted in a time-series stack of 25 years with over 98% coverage.

Establishment of Candidate Reference Pixels
Fire maps that were maintained by the state of Victoria's land management agency [36] were used to indicate the general extent of the three large fires.Candidate reference pixels were chosen via a systematic sampling process based on the Victorian Forest Monitoring Program (VFMP) plot network [37].The VFMP plot network consists of 786 2km by 2km plots that are distributed throughout public land in Victoria, stratified by bioregion and land tenure.In each plot, 10 random pixels were selected (resulting in 7860 pixels), and a team of six worked to manually interpret each pixel to establish its disturbance history (Figure 1 shows an example of the reference pixel sampling method).This was achieved by interrogating multiple lines of evidence, such as state fire records and high resolution imagery from Google Earth.Quality assurance was performed by an independent operator, who assessed 10% of all the pixels to evaluate the accuracy of the dataset (for details see Soto-Berelov et al. [38]).A total of 1391 pixels fell within the fire boundaries that were considered in this study.Of these, 1056 were classified as being disturbed by one or more of the three wildfires, and were subsequently used for the bulk of the analysis presented in this paper.In the section investigating different forest classes, an additional 5000 random pixels (with a minimum distance of 100 m) were selected inside the VFMP plots that fell within the fire history polygons, to ensure an adequate number of samples in each class.Visual inspection of the imagery indicated that, on balance, the majority of these were fire affected.

Landsat Spectral Indices
From the composite Landsat images, we generated the spectral indices that are shown in Table 2.
Table 2. Landsat spectral indices used in this paper, and a selection of pixel-based time-series studies using these indices (band numbers refer to Landsat TM and ETM+ bands).

Data Distributions of Pixels Pre and Post-Fire
To assess the sensitivity of each index in its response to fire, we conducted a number of tests on the 1056 disturbed reference pixels.For each index, we created image stacks covering 25 years and extracted the underlying values for each pixel of interest using the raster package [49] in R [50].The data was then grouped by relative years (e.g., year before fire, year of fire, year after fire, etc.).The aim was to compare the distributions of the pre-fire values and the post-fire values, and how they differ across indices (see Figure 2 for a conceptual diagram).To quantify the magnitude of the change between pre and post-fire values, we used the concept of effect size.Effect size refers to a family of statistical measures that are commonly used to measure the difference between two distributions in a standardized way, independent of sample size.For this study, we used Glass's delta, which is simply the difference in means between two groups, divided by the standard deviation of the control group [51].
where µ 1 is the mean of group 1 and µ 2 is the mean of group 2, and σ 1 is the standard deviation of group 1.In this exercise, the mean of group 1 (the control group) is the average value in a given index for all of the reference pixels in the 10 years prior to the fire (e.g., NDVI of 0.7).The mean of group 2 is the average value of all the pixels post-fire (e.g., NDVI of 0.3).The difference of −0.4 is then divided by the standard deviation of the control group (e.g., 0.1), which gives an effect size of −4.We chose to use the standard deviation of only the pre-fire values, rather than all of the values (as with Cohen's d), as this reflects the natural range of values for undisturbed forest in our study area.The effect size that is significant (practically speaking) will differ study to study.Cohen loosely defined effect sizes equaling 0.2 as small, 0.5 as medium, and 0.8 (or greater) as large [51].In our case, the variation in the mean values in the 10 pre-fire years is an indication of what effect size has practical significance, as this captures the natural fluctuations inherent in each index.
aim was to compare the distributions of the pre-fire values and the post-fire values, and how they differ across indices (see Figure 2 for a conceptual diagram).To quantify the magnitude of the change between pre and post-fire values, we used the concept of effect size.Effect size refers to a family of statistical measures that are commonly used to measure the difference between two distributions in a standardized way, independent of sample size.For this study, we used Glass's delta, which is simply the difference in means between two groups, divided by the standard deviation of the control group [51].

∆ = −
where μ 1 is the mean of group 1 and μ 2 is the mean of group 2, and σ 1 is the standard deviation of group 1.In this exercise, the mean of group 1 (the control group) is the average value in a given index for all of the reference pixels in the 10 years prior to the fire (e.g., NDVI of 0.7).The mean of group 2 is the average value of all the pixels post-fire (e.g., NDVI of 0.3).The difference of −0.4 is then divided by the standard deviation of the control group (e.g., 0.1), which gives an effect size of −4.We chose to use the standard deviation of only the pre-fire values, rather than all of the values (as with Cohen's d), as this reflects the natural range of values for undisturbed forest in our study area.The effect size that is significant (practically speaking) will differ study to study.Cohen loosely defined effect sizes equaling 0.2 as small, 0.5 as medium, and 0.8 (or greater) as large [51].In our case, the variation in the mean values in the 10 pre-fire years is an indication of what effect size has practical significance, as this captures the natural fluctuations inherent in each index.As well as the mean, we looked at how the standard deviation (SD) changed post-fire, the hypothesis being that a larger dispersion of post-fire values could be an indicator of which index may more accurately map fire severity (i.e., more classes or a greater range of values).The change in SD was calculated by dividing the SD post-fire by the SD pre-fire.We also calculated the percentage of the post-fire values that overlapped with the pre-fire values, with a lower percentage indicating better separation.Although this is somewhat captured in the effect size already, it is nevertheless interesting to consider the percentage of overlapping pixels, especially in terms of the change immediately after the fire when compared with that of one year later.

Spectral Response in Different Forest Classes
To determine how indices performed across different forest classes, we split the data based on tree height and canopy cover (as outlined in Table 2).The original classification was performed as part of the VFMP [52] and is used in State of the Forest reporting [53].As outlined earlier, to ensure an adequate number of samples in each forest class, in addition to the 1056 reference pixels, we generated a further 5000 random pixels in the VFMP plots, with the fire year being determined by As well as the mean, we looked at how the standard deviation (SD) changed post-fire, the hypothesis being that a larger dispersion of post-fire values could be an indicator of which index may more accurately map fire severity (i.e., more classes or a greater range of values).The change in SD was calculated by dividing the SD post-fire by the SD pre-fire.We also calculated the percentage of the post-fire values that overlapped with the pre-fire values, with a lower percentage indicating better separation.Although this is somewhat captured in the effect size already, it is nevertheless interesting to consider the percentage of overlapping pixels, especially in terms of the change immediately after the fire when compared with that of one year later.

Spectral Response in Different Forest Classes
To determine how indices performed across different forest classes, we split the data based on tree height and canopy cover (as outlined in Table 2).The original classification was performed as part of the VFMP [52] and is used in State of the Forest reporting [53].As outlined earlier, to ensure an adequate number of samples in each forest class, in addition to the 1056 reference pixels, we generated a further 5000 random pixels in the VFMP plots, with the fire year being determined by the fire history polygons [36].After removing those that did not fall within a class, we were left with 5759 pixels (Table 3, note that low tree height was uncommon in this study area and therefore not used).We then calculated the standardized means for each index in each forest class to establish the sensitivity of each index in the different forest systems.In addition, Analysis of Variance (ANOVA) tests between all of the pairs of forest classes (e.g., Medium Wood versus High Open, etc.) were conducted on a subsample of 250 pixels per class (to maintain class balance), to test for statistical significance between forest class distributions.Understanding forest regrowth and recovery following fire is essential for land managers to make informed management decisions.Free and open access to the long archive of Landsat data has created new opportunities for assessing the post-disturbance recovery of vegetation in terms of spectral response [19].Researchers have approached spectral recovery in different ways.Kennedy et al. [18] use a measure of recovery based on the difference between the each pixel's disturbance value and that of five years after disturbance, while Pickell et al. [11] look at recovery in terms of the number of years for the spectral index value to reach 80% of its pre-disturbance value.Although Landsat time-series cannot capture the full complexity of forest recovery, it enables large area assessments that are beyond the practicality of field based methods.In this context, we looked to evaluate the merits of each index in tracking post-fire spectral recovery.We compared indices by grouping reference pixels by year and by considering how each year's distribution post-fire relates to the overall pre-fire (undisturbed) distribution.The length of recovery, according to each index, was determined by calculating when the distribution mean of a post-fire year first reaches the lowest mean from the 10 years pre-fire.

Changes in Texture Pre and Post-Fire
Texture is not widely used in time-series studies [54], however, it has been a recognized image processing technique for many years [55,56].With advances in computer power, there has recently been interest in considering the spatio-temporal variables [57] in time-series studies.A change in texture pre-fire to post-fire is interesting in that it may assist in image classification; and, it may indicate ecological changes in the underlying forest.For example, following fire the forest may become more diverse (hence have greater textural variation), or it may become more homogenous (less variation).To capture some of the spatial variations ('texture') and how this manifests in different indices, we looked at each candidate pixel in relation to its neighbours.This was achieved by creating a 60 m buffer around each pixel and calculating the standard deviation of all the pixels in the buffer area pre and post-fire.A 90 m buffer was also trialed, and produced similar results; therefore, 60 m was considered adequate.Again, these values were standardized to delta using the distribution means and standard deviations, as outlined in Section 2.5.

Data Distributions of Pixels Pre and Post-Fire
Density histograms indicate the relative distribution of values pre-fire, directly after fire, and one year post-fire (Figure 3).Three different methods were used to quantify the information that is shown in the histograms.These results (Table 4) include the change in mean (standardized to delta ∆), the percentage overlap, and the change in SD.The lowest mean value from the 10 years prior to the fire is also presented as an indication of the natural undisturbed variation.For example, the standardized mean for NDVI for the fire year is −4.3, which is significantly lower than the lowest mean from the 10 years prior to the fire, which is −0.4.Our results show that the mean of the NDVI values changes the most directly after a fire, by −4.30, followed by TCA with −3.90, and NBR with −3.58.One year after fire, NBR has the greatest mean change, with −2.04, followed by NBR2 with −1.86, and NDMI with −1.63.NDVI has the smallest percentage of overlapping pixels directly following fire, with 14%, while one year later NBR2 shows greatest separation, with 39%.NBR shows the greatest change in SD both directly following fire and one year later, changing by a factor of 1.96 and 1.46, respectively.

Data Distributions of Pixels Pre and Post-Fire
Density histograms indicate the relative distribution of values pre-fire, directly after fire, and one year post-fire (Figure 3).Three different methods were used to quantify the information that is shown in the histograms.These results (Table 4) include the change in mean (standardized to delta Δ), the percentage overlap, and the change in SD.The lowest mean value from the 10 years prior to the fire is also presented as an indication of the natural undisturbed variation.For example, the standardized mean for NDVI for the fire year is −4.3, which is significantly lower than the lowest mean from the 10 years prior to the fire, which is −0.4.Our results show that the mean of the NDVI values changes the most directly after a fire, by −4.30, followed by TCA with −3.90, and NBR with −3.58.One year after fire, NBR has the greatest mean change, with −2.04, followed by NBR2 with −1.86, and NDMI with −1.63.NDVI has the smallest percentage of overlapping pixels directly following fire, with 14%, while one year later NBR2 shows greatest separation, with 39%.NBR shows the greatest change in SD both directly following fire and one year later, changing by a factor of 1.96 and 1.46, respectively.

Spectral Reponse in Different Forest Classes
For each of the eight indices, the change in mean directly following a fire was calculated for each forest class, again being standardized using the mean and standard deviation of the pre-fire values (Figure 4).In general, these results show that all of the indices are most responsive in woodland systems (low canopy cover).The wetness indices, particularly TCW, showed much less distinction in closed forest systems.As in the previous section, NDVI and TCA displayed the greatest changes, with NDVI shifting by as much as −4.8 in high woodland systems.In contrast, TCW only shifted by −0.5 in medium closed forests.NBR and NBR2 consistently occupied positions 3 and 4 in all of the forest classes.ANOVA tests on a subsample of 250 pixels per class showed that all the indices displayed significant differences between most forest classes, with p < 0.001 for all combinations, except for the following: Medium Closed and High Closed, and Medium Open and High Woodland, which none of the indices were able to clearly distinguish between; and, Medium Woodland and High Woodland, which SAVI was unable to distinguish between, with all of the other indices having p-values < 0.02.

Spectral Reponse in Different Forest Classes
For each of the eight indices, the change in mean directly following a fire was calculated for each forest class, again being standardized using the mean and standard deviation of the pre-fire values (Figure 4).In general, these results show that all of the indices are most responsive in woodland systems (low canopy cover).The wetness indices, particularly TCW, showed much less distinction in closed forest systems.As in the previous section, NDVI and TCA displayed the greatest changes, with NDVI shifting by as much as −4.8 in high woodland systems.In contrast, TCW only shifted by −0.5 in medium closed forests.NBR and NBR2 consistently occupied positions 3 and 4 in all of the forest classes.ANOVA tests on a subsample of 250 pixels per class showed that all the indices displayed significant differences between most forest classes, with p < 0.001 for all combinations, except for the following: Medium Closed and High Closed, and Medium Open and High Woodland, which none of the indices were able to clearly distinguish between; and, Medium Woodland and High Woodland, which SAVI was unable to distinguish between, with all of the other indices having p-values < 0.02.

Spectral Recovery Post-Fire
Figures 5 and 6 show the mean values for five years prior to a fire (for context) and nine years after (indicating recovery), for greenness and wetness based indices, respectively (note that pixels burned in 2009 did not contribute to the distributions of years 8 and 9). Figure 5 shows the greenness

Spectral Recovery Post-Fire
Figures 5 and 6 show the mean values for five years prior to a fire (for context) and nine years after (indicating recovery), for greenness and wetness based indices, respectively (note that pixels burned in 2009 did not contribute to the distributions of years 8 and 9). Figure 5 shows the greenness indices (particularly SAVI and TCG) almost returning to pre-fire levels three years after fire, although they do not technically pass the lowest pre-fire mean until year five.For wetness indices (Figure 6), the time to recover is longer, with NDMI reaching the lowest pre-fire mean at year seven, NBR and TCW at year eight.NBR2 does not reach pre-fire levels even after nine years, and interestingly, this index has a more consistent (smooth) recovery.Table 5 outlines the average number of years that each index takes to recover, defined by the year when the mean reaches the lowest mean from the ten years prior to fire disturbance.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 17 indices (particularly SAVI and TCG) almost returning to pre-fire levels three years after fire, although they do not technically pass the lowest pre-fire mean until year five.For wetness indices (Figure 6), the time to recover is longer, with NDMI reaching the lowest pre-fire mean at year seven, NBR and TCW at year eight.NBR2 does not reach pre-fire levels even after nine years, and interestingly, this index has a more consistent (smooth) recovery.Table 5 outlines the average number of years that each index takes to recover, defined by the year when the mean reaches the lowest mean from the ten years prior to fire disturbance.indices (particularly SAVI and TCG) almost returning to pre-fire levels three years after fire, although they do not technically pass the lowest pre-fire mean until year five.For wetness indices (Figure 6), the time to recover is longer, with NDMI reaching the lowest pre-fire mean at year seven, NBR and TCW at year eight.NBR2 does not reach pre-fire levels even after nine years, and interestingly, this index has a more consistent (smooth) recovery.Table 5 outlines the average number of years that each index takes to recover, defined by the year when the mean reaches the lowest mean from the ten years prior to fire disturbance.Results of the texture analysis are shown in Figures 7 and 8 for greenness and wetness based indices, respectively.Of the greenness indices, NDVI and TCA follow a similar trend, showing an increase in textural variation following a fire.In contrast, SAVI and TCG show less textural variation directly after a fire, but an increase after one year.NDVI and TCA, in particular, show that in the years following a fire, there is an increase in pixel variation, which gradually returns to pre-fire levels around eight or nine years after the fire.In the wetness indices, NBR, NDMI, and TCW all show an increase in textural variation directly after a fire, and return to pre-fire levels at around year four.These results are quite different to the recovery metrics that are presented earlier (based on individual pixels), where greenness indices returned to pre-fire levels before the wetness indices.NBR2 appears insensitive to textural variation, maintaining a similar level for the entire time series.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 17 increase in textural variation following a fire.In contrast, SAVI and TCG show less textural variation directly after a fire, but an increase after one year.NDVI and TCA, in particular, show that in the years following a fire, there is an increase in pixel variation, which gradually returns to pre-fire levels around eight or nine years after the fire.In the wetness indices, NBR, NDMI, and TCW all show an increase in textural variation directly after a fire, and return to pre-fire levels at around year four.These results are quite different to the recovery metrics that are presented earlier (based on individual pixels), where greenness indices returned to pre-fire levels before the wetness indices.NBR2 appears insensitive to textural variation, maintaining a similar level for the entire time series.

Discussion
For all of the indices, values that were measured directly following a fire were markedly different from the pre-fire (undisturbed) averages.With the exception of TCW, all of the indices showed a high degree of separation and a low percentage of overlap between pre and post-fire distributions.Greenness indices showed high sensitivity directly after a fire; however, one year later, increase in textural variation following a fire.In contrast, SAVI and TCG show less textural variation directly after a fire, but an increase after one year.NDVI and TCA, in particular, show that in the years following a fire, there is an increase in pixel variation, which gradually returns to pre-fire levels around eight or nine years after the fire.In the wetness indices, NBR, NDMI, and TCW all show an increase in textural variation directly after a fire, and return to pre-fire levels at around year four.These results are quite different to the recovery metrics that are presented earlier (based on individual pixels), where greenness indices returned to pre-fire levels before the wetness indices.NBR2 appears insensitive to textural variation, maintaining a similar level for the entire time series.

Discussion
For all of the indices, values that were measured directly following a fire were markedly different from the pre-fire (undisturbed) averages.With the exception of TCW, all of the indices showed a high degree of separation and a low percentage of overlap between pre and post-fire distributions.Greenness indices showed high sensitivity directly after a fire; however, one year later,

Discussion
For all of the indices, values that were measured directly following a fire were markedly different from the pre-fire (undisturbed) averages.With the exception of TCW, all of the indices showed a high degree of separation and a low percentage of overlap between pre and post-fire distributions.Greenness indices showed high sensitivity directly after a fire; however, one year later, they displayed much less distinction.In contrast, wetness indices experienced smaller differences directly following a fire event, but one year later had greater separation.These results were somewhat expected, and they align with findings in other studies [5,11].Furthermore, they suggest that in sclerophyll forests, vegetation quickly regains photosynthetic activity at the canopy level following a fire, with a large proportion of pixels returning to pre-fire levels within one year.This is most likely attributable to a combination of epicormic growth, as well as understory vegetation, such as grasses and non-woody plant matter.It is worth noting that TCA, which we have classed as a greenness index, appears to be more capable than the other greenness indices in capturing fire disturbances one year after the event.
Our results also clearly indicate greater dispersion of values in most indices following fire, with the exception of SAVI and TCG, which hardly change at all, and NDMI, which changes little.The standard deviations for both NBR and NDVI, for example, almost double following fire, and still maintain relatively high levels of dispersion one year later.Given that fires impact forests across a range of severity levels, a greater dispersion of values post-fire may indicate that the index is more suitable for mapping burn severity.Indeed, NBR (as noted earlier), has been used extensively for this purpose [14][15][16], although these authors concur that best results are usually found only in forested ecosystems.
When evaluating the performance of the indices across different forest systems, relatively speaking they performed similarly.That is, post-fire changes immediately after a fire were greatest in NDVI and TCA in all forest classes, followed by NBR and NBR2.We were not able to detect major differences between some forest classes (High Closed-Medium Closed and Medium Open-High Wood) in any index.However, there is clearly a distinction with regards to pre-and post-fire values between closed canopy forests versus woodland or open systems.This may indicate that Victoria's closed sclerophyll forests are more resilient to fire than their open counterparts; however, it could equally be a function of Landsat only capturing spectral changes of the canopy.More research into forest types (in terms of tree species) could provide further information in this domain.
In agreement with other studies [11], we found that wetness indices take longer than greenness indices to return to pre-disturbance levels (eight years vs. five years).Depending on the ecological variable of interest, there may be a preference to adopt the longer timeframes as more accurately representing forest recovery.While an index such as NDVI captures the initial return of vegetation, and correlates with biophysical parameters, such as the fraction of green vegetation cover and green leaf biomass [9], it is limited in its ability to represent structural attributes, which are often more important indicators of factors, such as biodiversity and carbon [11].In contrast, NBR and the other wetness indices are more closely aligned with forest moisture and structure through the utilization of SWIR bands.Other studies suggest that TCW is well suited to observe forest recovery due to its ability to track overall moisture content [43], however, in our study, we found it less reliable because of its low level of separation directly following a fire.Like Storey et al. [21], we found that NBR2 has extended recovery timeframes, and may be worth considering for future post-fire recovery studies.In southeast Australia, many eucalypts have the ability to survive low and moderate fire through epicormic resprouting (Figure 9), whereas after high intensity stand replacement fires, forest regrowth is dependent on new seedlings (Figure 10) [58], which naturally thin out as the forest matures.However, these recovery patterns are also species and location dependent.In this study, we looked at relatively few pixels across a very large area (3 million hectares), thus grouping all of the pixels together may not be the ideal approach for considering forest recovery dynamics.The texture analysis produced unexpected results.Whereas, in the pixel-based analysis it was the greenness indices that quickly returned to pre-fire levels, in the texture analysis, it was the wetness based indices.NBR2's lack of textural variation makes it unsuitable for this type of analysis, perhaps due to the high correlation between the SWIR bands.NDVI and TCA both indicated a relatively long recovery time in terms of the textural variation, returning to pre-fire levels in eight to nine years.This time period agrees with the wetness indices in the individual pixel analysis.This finding has some potentially useful ramifications.One is that there may be some additional information in terms of forest recovery that can be unlocked through the consideration of spatial variation, and two, given that variation appears in the red and near-infrared bands, this facilitates the use of a greater range of available data (e.g., Landsat MSS data going back to 1972, before the SWIR bands were introduced, or other satellite systems).Studies demonstrating improved classification accuracies with texture typically include a range of variables [59].In this study, the only texture variable investigated was that of standard deviation, which is one of many variables that are found in the literature; additional information may be available in other metrics.Including texture in pixel-based time-series is an unexplored area and there are further research opportunities in this domain, however there are considerable technical challenges to overcome in order to incorporate spatial components into a temporal analysis.

Conclusions
This paper presents a straight-forward method for comparing the merits of various spectral indices by considering all of the pixels as a single distribution.In this research, we made use of existing reference data to select our candidate pixels, but the method itself does not rely on detailed  The texture analysis produced unexpected results.Whereas, in the pixel-based analysis it was the greenness indices that quickly returned to pre-fire levels, in the texture analysis, it was the wetness based indices.NBR2's lack of textural variation makes it unsuitable for this type of analysis, perhaps due to the high correlation between the SWIR bands.NDVI and TCA both indicated a relatively long recovery time in terms of the textural variation, returning to pre-fire levels in eight to nine years.This time period agrees with the wetness indices in the individual pixel analysis.This finding has some potentially useful ramifications.One is that there may be some additional information in terms of forest recovery that can be unlocked through the consideration of spatial variation, and two, given that variation appears in the red and near-infrared bands, this facilitates the use of a greater range of available data (e.g., Landsat MSS data going back to 1972, before the SWIR bands were introduced, or other satellite systems).Studies demonstrating improved classification accuracies with texture typically include a range of variables [59].In this study, the only texture variable investigated was that of standard deviation, which is one of many variables that are found in the literature; additional information may be available in other metrics.Including texture in pixel-based time-series is an unexplored area and there are further research opportunities in this domain, however there are considerable technical challenges to overcome in order to incorporate spatial components into a temporal analysis.

Conclusions
This paper presents a straight-forward method for comparing the merits of various spectral indices by considering all of the pixels as a single distribution.In this research, we made use of existing reference data to select our candidate pixels, but the method itself does not rely on detailed The texture analysis produced unexpected results.Whereas, in the pixel-based analysis it was the greenness indices that quickly returned to pre-fire levels, in the texture analysis, it was the wetness based indices.NBR2's lack of textural variation makes it unsuitable for this type of analysis, perhaps due to the high correlation between the SWIR bands.NDVI and TCA both indicated a relatively long recovery time in terms of the textural variation, returning to pre-fire levels in eight to nine years.This time period agrees with the wetness indices in the individual pixel analysis.This finding has some potentially useful ramifications.One is that there may be some additional information in terms of forest recovery that can be unlocked through the consideration of spatial variation, and two, given that variation appears in the red and near-infrared bands, this facilitates the use of a greater range of available data (e.g., Landsat MSS data going back to 1972, before the SWIR bands were introduced, or other satellite systems).Studies demonstrating improved classification accuracies with texture typically include a range of variables [59].In this study, the only texture variable investigated was that of standard deviation, which is one of many variables that are found in the literature; additional information may be available in other metrics.Including texture in pixel-based time-series is an unexplored area and there are further research opportunities in this domain, however there are considerable technical challenges to overcome in order to incorporate spatial components into a temporal analysis.

Conclusions
This paper presents a straight-forward method for comparing the merits of various spectral indices by considering all of the pixels as a single distribution.In this research, we made use of existing reference data to select our candidate pixels, but the method itself does not rely on detailed reference data.The main advantage in using these particular pixels was that they had been systematically sampled, based on plots stratified by bioregion and forest tenure [37].Thus, they are an accurate reflection of the entire forest estate in the study area.However, by considering all of the pixels as equal participants to a single distribution, detailed information in individual pixels may be lost.Nonetheless, the purpose of the exercise was not to derive detailed information about forest dynamics, but to determine which indices may be best suited for this task.Of the indices that were tested, we consider NBR as the most reliable index for tracking fire disturbance and recovery in sclerophyll forests, due to its consistently high performance across the range of tests performed.Although NDVI and TCA showed greater discrimination between pre and post-fire pixels directly after a fire, NBR was better one year after a fire event.In addition, it presented longer recovery time-frames (an average of eight years), which is a better indicator of forest structure and the return of biomass.As computing power increases, it conceivably becomes less important to choose only one or a few indices, with ensembles of indices offering improved results [26].However, for large area mapping the literature suggests that we are not yet at that stage.In addition, as more and more data from various satellite sensors becomes available, the selection of appropriate indices will remain important.The methods that are presented here offer a simple solution for an evidence based selection process.Although we have primarily chosen simple indices that are used commonly in Landsat time-series, our methods are equally applicable to other indices and satellites, as well as being transferable to other ecosystems.

Figure 1 .
Figure 1.Study area (as indicated by the cross-hatched fire area), showing the location of the Victorian Forest Monitoring Program (VFMP) plots and example reference pixels.

Figure 1 .
Figure 1.Study area (as indicated by the cross-hatched fire area), showing the location of the Victorian Forest Monitoring Program (VFMP) plots and example reference pixels.

Figure 2 .
Figure 2. Conceptual diagram showing distributions of pre-fire and post-fire values.

Figure 2 .
Figure 2. Conceptual diagram showing distributions of pre-fire and post-fire values.

Figure 3 .
Figure 3. Density histograms showing the distributions of pre-fire values (blue), directly after fire (red), and one year post-fire (green).

Figure 3 .
Figure 3. Density histograms showing the distributions of pre-fire values (blue), directly after fire (red), and one year post-fire (green).

Figure 4 .
Figure 4. Change in mean directly following a fire, according to forest class (note values converted to positive).

Figure 4 .
Figure 4. Change in mean directly following a fire, according to forest class (note values converted to positive).

Figure 5 .
Figure 5. Mean values for greenness indices from five years prior to fire to nine years post-fire.

Figure 6 .Table 5 . 3 . 4 .
Figure 6.Mean values for wetness indices from five years prior to fire to nine years post-fire.

Figure 5 .
Figure 5. Mean values for greenness indices from five years prior to fire to nine years post-fire.

Figure 5 .
Figure 5. Mean values for greenness indices from five years prior to fire to nine years post-fire.

Figure 6 .Table 5 . 3 . 4 .
Figure 6.Mean values for wetness indices from five years prior to fire to nine years post-fire.

Figure 6 . 5 . 3 . 4 .
Figure 6.Mean values for wetness indices from five years prior to fire to nine years post-fire.

Figure 7 .
Figure 7. Mean values of textural variation for greenness indices from five years prior to fire to nine years post-fire.

Figure 8 .
Figure 8. Mean values of textural variation for wetness indices from five years prior to fire to nine years post-fire.

Figure 7 .
Figure 7. Mean values of textural variation for greenness indices from five years prior to fire to nine years post-fire.

Figure 7 .
Figure 7. Mean values of textural variation for greenness indices from five years prior to fire to nine years post-fire.

Figure 8 .
Figure 8. Mean values of textural variation for wetness indices from five years prior to fire to nine years post-fire.

Figure 8 .
Figure 8. Mean values of textural variation for wetness indices from five years prior to fire to nine years post-fire.

Figure 9 .
Figure 9. Open forest two years after moderate severity fire.

Figure 10 .
Figure 10.Closed forest nine years after high severity (stand replacement) fire.

Figure 9 . 17 Figure 9 .
Figure 9. Open forest two years after moderate severity fire.

Figure 10 .
Figure 10.Closed forest nine years after high severity (stand replacement) fire.

Figure 10 .
Figure 10.Closed forest nine years after high severity (stand replacement) fire.

Table 1 .
Native forest structural classes in Australia.

Table 3 .
Number of reference pixels in each forest class used in this study.

Table 4 .
Post-fire response of each index, shown as a standardized change in mean, percentage overlap, and relative change in standard deviation, with best results indicated in bold.

Table 4 .
Post-fire response of each index, shown as a standardized change in mean, percentage overlap, and relative change in standard deviation, with best results indicated in bold.