MAPPING DISTURBANCE DYNAMICS IN WET SCLEROPHYLL FORESTS USING TIME SERIES LANDSAT

In this study, we characterised the temporal-spectral patterns associated with identifying acute-severity disturbances and low-severity disturbances between 1985 and 2011 with the objective to test whether different disturbance agents within these categories can be identified with annual Landsat time series data. We analysed a representative State forest within the Central Highlands which has been exposed to a range of disturbances over the last 30 years, including timber harvesting (clearfell, selective and thinning) and fire (wildfire and prescribed burning). We fitted spectral time series models to annual normal burn ratio (NBR) and Tasseled Cap Indices (TCI), from which we extracted a range of disturbance and recovery metrics. With these metrics, three hierarchical random forest models were trained to 1) distinguish acute-severity disturbances from low-severity disturbances; 2a) attribute the disturbance agents most likely within the acute-severity class; 2b) and attribute the disturbance agents most likely within the low-severity class. Disturbance types (acute severity and low-severity) were successfully mapped with an overall accuracy of 72.9%, and the individual disturbance types were successfully attributed with overall accuracies ranging from 53.2% to 64.3%. Low-severity disturbance agents were successfully mapped with an overall accuracy of 80.2%, and individual agents were successfully attributed with overall accuracies ranging from 25.5% to 95.1. Acute-severity disturbance agents were successfully mapped with an overall accuracy of 95.4%, and individual agents were successfully attributed with overall accuracies ranging from 94.2% to 95.2%. Spectral metrics describing the disturbance magnitude were more important for distinguishing the disturbance agents than the post-disturbance response slope. Spectral changes associated with planned burning disturbances had generally lower magnitudes than selective harvesting. This study demonstrates the potential of landsat time series mapping for fire and timber harvesting disturbances at the agent level and highlights the need for distinguishing between agents to fully capture their impacts on ecosystem processes. * Corresponding author


INTRODUCTION
Disturbance regimes play an important role in wet sclerophyll forests in South-East Australia by renewing old and susceptible forests, recycling nutrients and supporting habitat structures (Attiwill, 1994).The two key forms of disturbance that occur within these forests are natural disturbance (primarily wildfire) and human disturbance (primarily timber harvesting).Other types of natural disturbance in these forests include windthrow and mechanical damage -particularly in the understoryresulting from snowstorms.Human disturbance also includes planned burning and a range of silvicultural management regimes.
There has been increasing debate in the literature on whether human disturbance through timber harvesting has altered interactions between wildfire and wet sclerophyll forests, resulting in more widespread fire outbreaks.Some studies support the hypothesis that timber harvesting increases fire risk and severity (Lindenmayer, 2010;Lindenmayer et al., 2011Lindenmayer et al., , 2009) ) and others oppose (Attiwill and Adams, 2013;Attiwill et al., 2014;Ferguson and Cheney, 2011).
In Victoria, the most comprehensive disturbance information at the landscape level is found in the State Fire History Database (SFHD) (Department of Environment Land Water and Planning, 2015a) and the State Logging History Database (SLHD) (Department of Environment Land Water and Planning, 2015b).These databases have employed a range of methods to document and map fire and logging disturbances.Unfortunately both of these databases have significant documented positional and attributional limitations (Department of Sustainability and Environment, 2009a;GHD, 2012;Phan and Kilinc, 2015).These limitations may be overcome by mapping disturbances using Landsat time-series data.It is hoped that this will increase the knowledge and understanding of landscape-causes and consequences of both natural and anthropogenic disturbances within these forests and better inform the debate.
Previous studies have shown that Landsat's spectral bands can be used to discriminate fire severity and logging intensity in wet sclerophyll forests in South-East Australia.Victorian studies utilising Landsat for fire severity mapping (Department of Sustainability and Environment, 2009b, 2007, 2003;Haywood and Sparkes, 2009) or timber harvesting (Lehmann et al., 2013;Miller et al., 1994;Woodgate and Black, 1988) have used spectral information from one or two images.However approaches based on single years or binary maps are often restricted in their ability to characterise the complex dynamics between wildfire, climate change and timber harvesting.Thus a more comprehensive mapping approach utilising longer time series and characterising the disturbance magnitude and duration would be beneficial.
Following the opening of the United States Geological Survey (USGS) Landsat archive and the related increase in capacity to produce time series (Wulder et al., 2012), Landsat time series have increasingly been used at the regional scale to map a range of disturbances (timber harvesting, wildfires and insect outbreaks) using pixel-based time-series methods (White et al., 2014).The adoption of these techniques by Australian forest agencies has been limited.This has been partly due to the computational complexity of some of the procedures, use of proprietary software and the empirical nature of the customised requirements such as the trial and error basis for determining the optimal parameterization for the segmentation of the pixel time series (e.g.(Kennedy et al., 2007)).Nevertheless, with increasing availability of Landsat imagery and cloud computing (Wulder and Coops, 2014), coupled with diminishing availability of skilled photo interpreters (Haywood and Stone, 2011), there is increasing interest in South-East Australia for analysing pixel time series to better understand the ecological dynamics of fire, logging and their interactions within wet sclerophyll forests.
As significant research remains to be done before fully automated landscape level forest disturbance mapping can be achieved, the general approach adopted here has been to develop a semi-automated pixel-time-series-based method which is as practical as possible.Thus the interim goal -rather than trying to replace existing databases and associated methods -should be to support them in generating more timely, consistent (temporal and spatial) and accurate products.New and/or better tools are required to produce incremental improvements in these areas.It is not necessary for these tools to provide final solutions or 100% correct results, they simply need to be tools that are useful and that can be easily corrected when things go awry.They should be simple to apply, not require expensive equipment, not substantially alter the existing mapping workflow, nor involve inordinate fine-tuning by the interpreter.
Although there are a number of existing pixel-level disturbance mapping tools available in the literature for identifying forest dynamics (Karfs et al., 2004;Kennedy et al., 2010;Udelhoven, 2011;Verbesselt et al., 2010), they have all been developed overseas for non-eucalypt forests and significant effort is required to become familiar with these algorithms and associated proprietary software modules.As a consequence, an alternative approach that develops an integrated workflow process utilising standard (maintained) open-source software and packages was applied to this study.To ease the computational burden and storage requirements it was decided to limit the approach to utilize annual Landsat time series.The overall goal was to determine the capacity of readily available open-source tools to model spectral-temporal pixel time-series from annual Landsat time series to map fire and timber disturbance dynamics within wet sclerophyll forests in South-East Australia.Specific objectives were to: 1.
test how well fire and timber harvesting disturbances can be distinguished with annual Landsat time series and open-source software; 2. characterise the spectral-temporal pixel-time-series of fire and timber harvesting disturbances with respect to severity magnitude and spectral recovery; and 3. map the spatial and temporal pattern of fire and timber harvesting disturbances using open-source software.

OPEN-SOURCE SOFTWARE
By adopting an open-source approach for spatial data management, processing and analysis, users such as forest management agencies can benefit from freely available software products and access to source code through which new algorithms can be integrated and manipulated.The key opensource software utilised within this study are outlined below.

GRASS
The Geographical Resources Analysis Support System (GRASS) platform (GRASS Development Team, 2012) was chosen due to its popularity within the open-source community and because it fully integrates with the open-source statistical software package, R (R Development Core Team, 2012), along with the python scripting language (van Rossen, 1995).
It is an open-source geographical information system (GIS) capable of handling raster, topological vector, image processing and graphic data.Released under the GNU General Public License (GPL), GRASS is developed by a multi-national group of developers and is one of the eight initial software projects of the Open Source Geospatial Foundation.GRASS has a modular structure into which may be plugged new routines programmed in a variety of languages (e.g., Python, C, shell), and there are over 300 modules and more than 100 addon modules for the creation, manipulation and visualisation of both raster and vector data.The GRASS modules are designed under the UNIX philosophy (i.e., that programs work together and handle text streams) and can be combined using scripting to create more complex or specialized modules by a user.GRASS supports an extensive range of raster and vector formats through GDAL/OGR libraries, including OGC-conformal (Open Geospatial Consortium) Simple Features for interoperability with other GIS.

R
R is an open-source language and software environment commonly used in research fields for statistical computing and graphics.One of the main advantages of R is its objectorientated approach, which allows results of statistical procedures to be stored as objects and used as input in further computations.R is a simple and effective formal complete programming language, and the R environment is, therefore, highly extensible.GRASS and R software can be integrated through the R package, spgrass (Bivand, 2007), an interface allowing GRASS functions to be implemented within R code and data to be easily exchanged between the two software packages.In addition, R package, raster (Hijmans and van Etten, 2012), has functions for creating, reading, manipulating, and writing raster data.The package also implements raster algebra and most functions for raster data manipulations that are common in GIS.

Python
Python is an object-orientated high-level programming language that is widely used as a scripting language in the spatial analysis environment.Python's popularity has led to the creation of many useful libraries, increasing its flexibility and interoperability, and it has well developed modules for linking with GRASS and R.

Geographic and biophysical characteristics
Our study area is the  1988).

Natural and anthropogenic disturbances
Wildfire is the major natural disturbance associated with the study area.Several fires have occurred within the study area over the past 150 years, the most extensive being in 1926; 1939 and the recent extreme fire event of 2009 (Price and Bradstock, 2012).
The study area is also subject to intensive hardwood timber harvesting.Large-scale timber cutting, generally selective harvesting and sawmilling occurred in these forests in the latter part of the nineteenth and early twentieth centuries.Large-scale salvage operations followed major wildfires, particularly the extensive 1939 fires.Since the 1960s, clearfellling has been the major silvicultural system practised (Squire et al., 1991).
Figure1: Study area located in Central Highlands of Victoria, Australia.

DATA AND METHODS
A general overview of the methods used in this study is shown in Figure 2.Each of the steps taken in the study are described in detail below.

Landsat data and pre-processing
We downloaded all level-1 terrain-corrected (L1T) Landsat data acquired between 1 January 1984 to 28 February 2011 with cloud cover < 90% from the USGS archive for path 92 and row 086.Each image was first screened for cloud and cloud shadow using Fmask (Zhu and Woodcock, 2012) and converted to surface reflectance using LEDAPS for the 23 year time period (Masek et al., 2006).
To minimise the effect of phenology and data gaps caused by atmospheric interference, we constructed annual anniversarydate, best observation composites using all cloud free observations within a pre-defined seasonal window, following the method of Kennedy et al. (2010).For building the best observation composites, we defined the seasonal window as ± 60 days around February 15.
Using the outlined selection criteria we had 100% coverage for the annual composite stack.All pre-processing steps were conducted in GRASS using a range of standard and custom modules.

Landsat vegetation indices
In this study we utilised four indices which are responsive to different vegetation cover/disturbance properties including vegetation greenness, moisture content, canopy structure, and exposed soil signal.We generated Landsat time series stacks using the Normalised Burn Ratio (NBR, Key and Benson, 2005), Tasseled Cap Wetness (TCW, Crist and Cicone, 1984)), Tasseled Cap Brightness (TCB, Crist and Cicone, 1984) and Tasseled Cap Angle (TCA, Powell et al., 2010).The creation of the vegetation indices was conducted in GRASS using the mapcalc function.

Landsat time series analysis
We carried out the time series analysis using a number of standard R packages.The time series analysis was conducted to: 1) extract spectral time series for each pixel 2) statistically identify and fit structural equations and 3) extract summary information from trends.

Extraction of time series for each pixel
Once the vegetation indices were calculated, the image stack was loaded into R using the raster package (Hijmans and van Etten, 2012) using the RasterStack function.Spectral time series were then extracted as a vector for processing using the calc function within the raster package.Spectral values for each year can be taken from any arbitrary window kernel centred on the pixel of interest; in this study, we chose to use the mean value in a 3 x 3 window as a compromise between spatial detail and robustness to pixel misregistration across images in the stack.

Statistically identify and fitting trends
Once a consistent spectral time series is extracted from the image stack we used the bfast package (Verbesselt et al., 2010) to fit a structural breakpoints from a linear regression model.Following previous studies using the bfast package (DeVries et al., 2015;Verbesselt et al., 2012), we assigned a value of 0.25n to h.A structural breakpoint is declared when the null hypothesis of structural stability (i.e.stability of the seasonality pattern) is rejected (Verbesselt et al., 2012;Zeileis et al., 2005).
The decision to reject this null hypothesis is based on a boundary condition which is set according to a 5% probability level following the Functional Central Limit Theorem (see (Leisch et al., 2000) for more information on how this boundary function is computed).

Extract summary information from trends
Once a trend was fitted to the vegetation indices time series we derived the following set of metrics: 1.For pixels without a breakpoint detected, the slope and intercept of the linear trend of the time series was extracted 2. For pixels with a breakpoint detected, the magnitude of the breakpoint was calculated, the date of the breakpoint, the slope and intercept for the line segments before and after the breakpoint were extracted.

Forest disturbance mapping
We followed a two-phase classification approach based on Senf et al. (2015) to map spatial and temporal patterns of natural and anthropogenic disturbance: First, we classified the Landsat time series disturbance and recovery metrics into three classes: 1) acute high severity disturbances; 2) low severity disturbances; and 3) undisturbed areas.We refer to this classification phase as disturbance type classification.Second, we assigned all pixels identified within the acute high severity in the first classification phase a likelihood of being disturbed by either wildfire or clearfell timber harvesting; we assigned all pixels identified within the low severity disturbance a likelihood of being disturbed by planned burning, selective timber harvesting, insects or drought.We refer to this classification phase as disturbance agent attribution.

Phase one: disturbance type classification:
In the first classification phase, we used the Landsat time series disturbance and recovery metrics to classify forest changes into 1) acute high severity disturbances, 2) Low severity disturbances, and 3) undisturbed forest.High severity disturbances (such as clear-fell timber harvest and wildfires) behave differently in spectral and temporal space than low severity disturbances (such as selective harvesting and planburning disturbances) which makes them distinguishable with Landsat time series.While some of the low severity disturbances can eventually lead to complete stand mortality, spectral change magnitudes associated with clear-fell timber harvesting and wildfire disturbances are usually significantly higher (Schroeder et al., 2011).As a reference data, we randomly selected and labelled 500 pixels closely following the approach by Cohen et al. ( 2010).
For identifying and labelling disturbances in the reference pixels, we used Landsat imagery, Landsat spectral time series plots, high resolution imagery (Rapideye or GoogleEarth Imagery), the SFHD and SLHD databases.
The SFHD (1903SFHD ( -2015) ) contains polygon-level data on fire perimeter, type (wildfire or planned burning), and for a limited number fires, mapping methodology and severity information.
The SLHD (1879-2015) collects polygon-level data on extent, silvicultural operation, forest type, start/end dates of logging event and mapping methodology.Phan and Kilnic (2015) found that almost 40% of the state fire history database contained missing or incorrect information regarding the date stamps on the fires.However, they did find that recent records (2006-2014) have a higher quality, with only 7% containing missing or incorrect data.To reduce uncertainties in our analysis we only utilised data in the state fire history database that could also be linked with entries in the state bushfire ignitions point database (Department of Environment Land Water and Planning, 2015c), the state planned burning ignitions database (Department of Environment Land Water and Planning, 2015d) or the Country Fire Incident Reporting System (Country Fire Authority, 2015).Some of the variability in quality within the database can be attributed to the wide diversity in base data and mapping methodologies utilised (on-screen digitising using aerial photography, field GPS data capture, ground observations, , thermal line scanner mapping, automated image interpretation using Rapideye, Spot and Landsat imagery, transfer from hard copy maps).As the extreme wildfire event of 2009 covered in excess 60% of the study area, this single disturbance event has been removed from the analysis.
The state logging history database also has a range of documented positional and attributional limitations (Department of Sustainability and Environment, 2009a;GHD, 2012).The data is subject to a certain observer bias, variations in base data, and interpreter/analysis experience.It has not been uncommon to find 5-10% of the records omitted or duplicated.For the last 15-20 years the state logging history database base data has been sourced from GPS ground survey of the logging boundary.The accuracy of the resultant mapped polygon has improved but is still reliant on several factors such as GPU unit specifications, satellite positions, atmospheric conditions and natural barriers to the signal.Prior to the use of GPS, the logging boundaries were estimated from sketch mapping on 1:10000 hard copy maps.
In total, 47 pixels were identified as acute high severity disturbance, 70 pixels were identified as low severity disturbances and 363 were identified as undisturbed.A small proportion (20 pixels) could not clearly be assigned to one of these categories.Using the reference pixels, we trained a random forest classification model (Breiman, 2001) provided by the randomForest package (Liaw and Wiener, 2002) within R. The random forest model was validated using the out-of-bag confusion matrix (Breiman, 2001), from which we estimated overall, user's and producer's accuracies, as well as errors of omission and commission.

Phase two: disturbance agent attribution:
Following the disturbance type mapping in phase 1 (Section 4.3.1)we estimated for: 1.Each acute high severity disturbance pixel the probability of being disturbed by wildfire or clear-fell timber harvesting, respectively; and 2.
Each low severity disturbance pixel the probability of being disturbed by selective timber harvesting or low severity wildfire/planned burning (for the purposes of this paper low severity wildfire and planned burning were collapsed into a single class).Creating a continuous probability of class presence can offer greater flexibility from a forest management perspective than discrete classes (Wulder et al., 2006).For this purpose we calibrated two additional random forest models with additional reference datasets from the state fire and logging history databases.

Attribution of acute high severity disturbance
For the acute high severity disturbance attribution, we selected all pixels covered by either wildfire or clear-fell harvesting polygons from the SFHD and SLHD databases.

Attribution of low severity disturbance
For the low severity disturbance attribution, we selected all pixels covered by either planned burning or selective harvesting polygons from the SFHD and SLHD databases.

Classification of disturbance types
The disturbance classification yielded an overall accuracy of 72.9% (Table 1), 1), with the highest user's and producer's accuracies in the undisturbed class (92.7% and 77.1%, respectively), slightly lower user's and producer's accuracies for the acute severity class (56.3% and 64.3%, respectively), and moderate accuracies for the low severity class (25.5% and 53.2%, respectively).Class confusion was highest between low severity disturbance areas and undisturbed areas.In total, 14.6% of the forested area contained acute severity and 9.8% contained low severity disturbances.Most of the forested area in the study area (75.6%) was stable over the study period.The classification map (Figure 4) was used to identify acute severity and low severity areas for the following results.
The confusion matrix is derived from the out-of-bag sample of the random forest model.

Disturbance agent attribution
The binary classification of low severity wildfire/planned burning and selective timber/thinning harvesting disturbances (using a probability threshold of p = 0.5) achieved an overall accuracy of 80% (Table 2), indicating that the attribution of these two agents is much more difficult that the acute severity agents (Table 3).The user's accuracy for the selective logging was quite low, which means this disturbance agent was overestimated in the final mapped product.The binary classification of clearfell timber harvesting and wildfire disturbances (using a probability threshold of p=0.5) achieved an overall accuracy of 95.4% (Table 3), indicating that these two agents can be reliable distinguished using disturbance and recovery metrics derived from Landsat time series.

DISCUSSION
This study demonstrates the feasibility of using an open-source framework for constructing and evaluating a spectral pixel time series model and its implementation to produce an accurate operational land management agency forest disturbance map.The framework established successfully integrates freely available spatial data-pre-processed and collated in GRASS-into the R statistical analysis environment.After construction and validation of an spectral time series segmentation, the resulting model was implemented in GRASS using an R-GRASS interface package, spgrass (Bivand, 2007), before finally using GRASS to filter the forest prediction map and apply the minimum mapping unit of the adopted forest definition to the final forest extent spatial product.

CONCLUSION
In this study we characterised acute high severity and low severity disturbance in South-East Australia, using a wellestablished Landsat-based time series technique.From our results, we conclude that Landsat can be utilised to reliably distinguish between acute severity disturbance agents (clearfellng and wildfire) in our study region, using specific spectra time-series features.However, more research is needed in distinguishing between the low severity disturbance agents (low severity wildfire/planned burning and selective logging).
The resulting maps and estimates offer a combined and detailed picture of disturbance dynamics in our study region through quantifying both the temporal and spatial dynamics.These otherwise unavailable spatially explicit and quality assured maps can help inform science and management needs.

Figure 2 :
Figure 2: Flowchart outlining the main steps implemented in this study 4.1 Forest population mask Similar to Hansen et al. (2013) we used a forest/non-forest mask to avoid confusion between forest disturbances and other land cover dynamics.The forest/non-forest mask used was that created by Mellor et al. ( 2013).

Figure 3 :
Figure 3: Example of pixel trajectories with related Landsat imagery for (A) no disturbance (B) acute high severity (clearfell) (C) low severity (planned burning)

Figure 4 :
Figure 4: Map derived in the disturbance classification phase showing undisturbed areas, acute severity disturbances, and low severity disturbances.

Figure 5 :
Figure 5: Mapped probability of (a) selective harvesting and (b) low severity wildfire/planned burning.
Toolangi State Forest and surrounding area.This forest is located approximated 80 km north east of Melbourne, in the Victorian Central Highlands, South-East Australia.The total area of the Toolangi State forest is approximately 40,000 hectares.The total area of the study area is 180,000 hectares.A high proportion of this mountainous area supports wet sclerophyll forests, dominated by Eucalyptus regnans (Mountain Ash).The area was selected to represent a variety of ash-forest types, forest conditions and disturbances.As mentioned previously, these forests are currently at the centre of a debate on whether timber harvesting in the region increases fire risk and severity.
The area experiences a cool temperate climate, with mild summers and cool winters.Average annual rainfall exceeds 1200 mm over most of the area.Soils tend to be free draining, friable, brown gradational, have high water holding capacities, and have developed on a variety of volcanic parent rock materials (Department of Natural Resources and Environment,

Table 1 :
Validation of the first classification phase (disturbance type classification), which distinguishes acute high severity disturbances, low severity disturbances and undisturbed areas

Table 2 :
Confusion matrix for predicting disturbance agents in low severity disturbance classes

Table 3 :
Confusion matrix for predicting disturbance agents in acute severity disturbance