Spatial Knowledge and Information Canada, 2019, 7(2), 4 A functional data analysis approach for characterizing spatial-temporal patterns of landscape disturbance and recovery from remotely sensed data MATHIEU L. BOURBONNAIS MICHAEL A. WULDER NICHOLAS C. COOPS Earth, Environmental and Canadian Forest Service Integrated Remote Sensing Geographic Sciences Natural Resources Canada Studio University of British University of British Columbia Okanagan JOANNE C. WHITE Columbia Canadian Forest Service TRISALYN A. NELSON Natural Resources Canada FAROUK NATHOO School of Geographical Mathematics & Statistics Sciences and Urban Planning GEORDIE W. HOBART University of Victoria Arizona State University Canadian Forest Service Natural Resources Canada CHRIS T. DARIMONT GORDON B. STENHOUSE Applied Conservation Grizzly Bear Program TXOMIN HERMOSILLA Science Lab, Department of Foothills Research Institute Integrated Remote Sensing Geography Studio University of Victoria University of British Columbia ABSTRACT mixture model. The resulting eight watershed clusters were mapped with mean Contemporary landscape regionalization functions representing unique temporal approaches, frequently used to summarize trajectories of disturbance and recovery. and visualize complex spatial patterns and There was considerable variability in disturbance regimes, often do not account disturbance amplitude among the clusters for the temporal component which may which increased markedly in the mid-1990’s provide important insight on disturbance, while remaining low in parks and protected recovery, and change in ecological areas. The regionalization highlights unique processes. The objective of this research was temporal trajectories of disturbance and to employ novel statistical approaches in recovery driven by anthropogenic and functional data analysis to quantify and natural disturbances and enables insight cluster spatial-temporal patterns of regarding how cumulative spatial landscape disturbance and recovery in 223 disturbance patterns evolve through time. watersheds using a Landsat disturbance time series from 1985 – 2011 in western 1. Introduction Alberta, Canada. Cumulative spatial patterns of disturbance, representing the Terrestrial ecosystems are subject to a range proportion, arrangement, size, and number of natural and anthropogenic disturbances of disturbances per watershed, were that influence landscape dynamics and modelled as functions and scores from a heterogeneity. In North America, the functional principal component analysis frequency, extent, and severity of natural were clustered using a Gaussian finite disturbances, including forest fires and 2 Functional data analysis regionalization insect infestations, has been increasing due Alberta, Canada from 1985 to 2011 using to anthropogenic influences and climate Landsat disturbance time series data change (Turner, 2010). Similarly, (Hermosilla et al., 2015). Methods in anthropogenic activities and anthropogenic functional data analysis (FDA) are pressures on many terrestrial ecosystems specifically designed to characterize are growing as resource extraction activities, multivariate high-dimensional time series including forest harvest, road network data (Ramsay & Silverman, 2005). Using development, and energy development and the FDA framework, our regionalization mining contribute to land use change and identifies unique temporal trajectories of landscape fragmentation (Pickell et al., cumulative disturbance patterns 2016). Cumulatively, landscape disturbance representing underlying distributions and is temporally dynamic given post- spatial-temporal dynamics of specific disturbance recovery, regeneration, and natural and anthropogenic disturbance succession. As such, monitoring and types, including forest fires, harvest, and quantifying how spatial patterns of natural roads (Bourbonnais et al., 2017). and anthropogenic landscape disturbance change over time is critical for 2. Methods and Data understanding how ecological processes are influenced by disturbance and recovery. 2.1 Landsat data and disturbance Change detection and attribution of pattern metrics disturbance from remotely sensed time series data provide opportunities to develop The study used a novel Canada-wide new hypotheses on disturbance recovery landscape disturbance time series derived and land cover change. The spatial from a best-available pixel Landsat data resolution and longevity of the Landsat product where disturbances, including mission, in particular, allows detection of forest harvest, oil and gas well-sites, roads, landscape alterations that are the result of a forest fires, and non-stand replacing given management or land use decision over disturbances (e.g., insects and drought) large areas in a systematic fashion (Wulder were detected and attributed annually from et al., 2012). While regionalization 1985 – 2011 (Hermosilla et al., 2015). Using approaches, where geographic entities are the Landsat disturbance time series, spatial grouped based on common factors to patterns of landscape disturbance were summarize complex landscape and quantified annually using the proportion environmental factors (Hargrove and area disturbed, the probability of Hoffman 2004), have been developed to disturbance adjacency, the mean characterize spatial patterns of landscape disturbance patch area, and the number of disturbance (e.g., Long et al., 2010), the disturbance patches in 223 watersheds in temporal dynamics of disturbance and western Alberta. Watersheds were selected recovery are often left unaccounted which as the landscape unit of analysis for the can influence interpretation of resulting regionalization as they are commonly used patterns (Pickell et al., 2016). as an environmentally relevant scale for monitoring forest and land cover changes The goal of this study is to characterize (Wulder et al., 2009). The disturbance disturbance as a temporally dynamic, pattern metrics were adjusted annually to allowing us to quantify and map cumulative account for recovery by comparing the patterns of landscape disturbance while normalized burn ratio (NBR = (B4- simultaneously accounting for recovery. To B7)/(B4+B7) where B4 and B7 correspond this end, we develop a novel functional data to Landsat bands 4 – near-infrared – and 7 analysis regionalization of landscape – short-wave infrared, respectively) from disturbance in 223 watersheds in western the pre- and post-disturbance periods (Key & Benson, 2006). A disturbance pixel was Functional data analysis regionalization 3 considered recovered, and subsequently our regionalization. We regionalized masked from the annual disturbance watersheds with common disturbance pattern metrics, when the post-disturbance patterns by clustering the FPC scores using NBR values reached 80% of the mean pixel Gaussian finite mixture models with the NBR values from the two years preceding optimal number of groups selected using the disturbance (Pickell et al., 2016). negative of the Bayesian Information Criterion (Fraley & Raftery, 2002). The 2.2 Functional data analysis clustered watersheds were then mapped and regionalization compared using the mean disturbance pattern metric curves by region. We further In the FDA framework, discrete time series explored variability in pattern metrics of observations (i.e., disturbance pattern attributed disturbances (fire, harvest, roads, metrics) are considered to arise through the well-sites, and non-stand replacing) for each regular sampling of a smooth function (i.e., watershed cluster using a functional curve) rather than thought of as a analysis of variance (FANOVA) by realization from a multivariate distribution comparing the mean curves based on shape (Ramsay & Silverman, 2005). Following the and temporal variability (Ramsay & FDA approach, the time series of discrete Silverman, 2005). disturbance pattern metrics in each watershed were converted to curves using B- 3. Results splines as the basis function. We used a functional principal component analysis Three FPCA scores were required to explain (FPCA), which estimates a set of eigenvalue- 90% of the variance in the proportion eigenfunction pairs, to quantify the primary disturbance, probability of disturbance modes of temporal variation among the adjacency, and mean disturbance patch curves for each of four disturbance pattern area, and two scores for the number of metrics (Ramsay & Silverman, 2005). For disturbance patches (Figure 1). Amplitude each of the four disturbance pattern metrics, in the first FPCA score, representing the we computed the minimum number of greatest deviation of the curve from the FPCA scores, representing the difference mean, generally increased markedly from the mean disturbance pattern curve, beginning in the mid-1990’s characterizing required to explain 90% of the functional differing disturbance pattern trajectories variance in the curves. The FPCA scores (n = resulting from a rapid increase in resource 11), which represent the primary modes of extraction and industrial activity in the temporal variation in the disturbance region. pattern metric curves, formed the basis for 4 Functional data analysis regionalization Figure 1. Functional principal components analysis (FPCA) for disturbance pattern metric curves of proportion disturbance (A), probability of disturbance adjacency (B), mean disturbance patch size (C) and number of disturbance patches (D). The upper plot maps the score with the highest absolute value for each watershed. The lower panel shows the mean of the fitted disturbance pattern metric curve (solid black line) and how the amplitude of the mean curve varies if the FPCA curve is added (+) or subtracted (-). The Gaussian finite mixture model throughout the study period with notable incorporating the eleven FPCA scores with recovery beginning circa 2005. Conversely, the greatest support (BIC = -2095.35) regions occurring primarily in parks and resulted in eight disturbance pattern regions protected areas (clusters 4, 7, and 8) had the (Figure 2). Watersheds in clusters 1 lowest overall disturbance amplitude. (35.92%), 5 (16.33%), and 4 (13.46%) Interestingly, while the mean proportion represented the greatest proportion of the disturbance, probability of disturbance study area. The amplitude (i.e., vertical) and adjacency, and mean disturbance patch size phase (i.e., horizonal) variability of the four- curves demonstrated periods of recovery disturbance pattern metric mean curves (i.e., periods of decreasing amplitude in the characterized periods of increasing curve), the number of disturbance patches disturbance (i.e., increasing curve generally increased over time suggesting amplitude) and spectral recovery (i.e., spatial variability in recovery may result in a decreasing curve amplitude) among the complex spatial mosaic of patches in watershed regions. Watersheds in clusters 1, different successional states (Gómez et al., 2, 3, 5, and 6 had increasing amplitude 2011). Functional data analysis regionalization 5 Figure 2. Mean curves by cluster for the proportion disturbance (A), the probability of disturbance adjacency (B), the mean disturbance patch area (C), and the number of disturbance patch (D) pattern metrics. Each mean curve is associated with the watersheds mapped by cluster membership (E). Parks and protected areas are shown in green. and recovery as a single continuous function The watershed clusters also revealed can reveal properties of the underlying variability in the amplitude and phase of the ecological processes and how patterns of mean disturbance pattern metrics for the landscape disturbance evolve over a attributed disturbance types compared continuum (Pickell et al., 2016). However, it using FANOVA (Figure 3). Mean forest fire is difficult to quantify landscape disturbance curves were significantly different (p < 0.05) cumulatively and to account for the among the watershed clusters, with large temporal dynamics of disturbance and forest fires prevalent in clusters 2 and 5. recovery, as well as the interaction of Trajectories representing the mean multiple sources of natural and proportion area disturbed and number of anthropogenic disturbance. Using the novel disturbed patches for forest harvest, roads, FDA approach described here, regional and well-sites were also significantly spatial-temporal disturbance patterns can different among the clusters, and had the be interpreted through representative greatest amplitude in clusters 6, 1, 5, and 2 disturbance trajectories which illuminate representing watersheds primarily outside the different disturbance processes and of parks and protected areas. indicate where and when anthropogenic disturbance is the dominant driver of 4. Conclusion observed patterns in the study area. As new time series of disturbance and land cover data become increasingly available, FDA- While piece-wise properties of curves have based approaches can be useful for been employed to detect and quantify quantifying and summarizing complex disturbance patterns (Gómez et al., 2011), spatial-temporal landscape patterns. modelling patterns of landscape disturbance 6 Functional data analysis regionalization Figure 3. Results of the functional analysis of variance showing mean curves of the proportion area disturbed (Pd), probability of disturbance adjacency (Pdd), mean disturbance patch area (Mdpa), and number of disturbance patches (Ndp), for the attributed disturbances (fire, harvest, non-stand replacing, roads, and well-sites) by cluster. Fraley, C., & Raftery, A. E. (2002). Model- based clustering, discriminant analysis, Acknowledgements and density estimation. Journal of the American Statistical Association, 97(458), This work was supported by the Natural 611–631. Sciences and Engineering Research Council of Canada through a Collaborative Research Gómez, C., White, J. C., & Wulder, M. A. and Development Grant (CRDPJ 486174- (2011). Characterizing the state and 15), a Canadian Graduate Scholarship, and processes of change in a dynamic forest by the Foothills Research Institute Grizzly environment using hierarchical spatio- Bear Project and its many funding partners. temporal segmentation. Remote Sensing of Environment, 115(7), 1665–1679. References Hermosilla, T., Wulder, M. A., White, J. C., Bourbonnais, M. L., Nelson, T. A., Coops, N. C., & Hobart, G. W. (2015). Stenhouse, G. B., Wulder, M. A., White, J. Regional detection, characterization, and C., Hobart, G. W., Hermosilla, T., Coops, attribution of annual forest change from N. C., Nathoo, F., & Darimont, C. T. 1984 to 2012 using Landsat-derived time- (2017). Characterizing spatial-temporal series metrics. Remote Sensing of patterns of landscape disturbance and Environment, 170, 121–132. recovery in western Alberta, Canada using a functional data analysis approach and Key, C.H., Benson, N.C., 2006. Landscape remotely sensed data. Ecological Assessment: Sampling and Analysis Informatics, 39, 140–150. Methods. USDA for. Serv. Gen. Tech. Rep. RMRS-GTR-164-CD, 1–55. Functional data analysis regionalization 7 Long, J. A., Nelson, T. A., & Wulder, M. A. Wulder, M.A., White, J.C., Andrew, M.E., (2010). Regionalization of landscape Seitz, N.E., Coops, N.C., 2009. Forest pattern indices using multivariate cluster fragmentation, structure, and age analysis. Environmental Management, characteristics as a legacy of forest 46(1), 134–142. management. Forest Ecology and Management. 258, 1938–1949. Pickell, P. D., Coops, N. C., Gergel, S. E., Andison, D. W., & Marshall, L. (2016). Wulder, M. A., Masek, J. G., Cohen, W. B., Evolution of Canada’s boreal forest spatial Loveland, T. R., & Woodcock, C. E. (2012). patterns as seen from space. PLoS ONE, Opening the archive: How free data has 11(7), e0157736. enabled the science and monitoring Ramsay, J. O., & Silverman, B. W. (2005). promise of Landsat. Remote Sensing of Functional data analysis. New York, NY: Environment, 122, 2–10. Springer. Turner, M. G. (2010). Disturbance and landscape dynamics in a changing world. Ecology, 91(10), 2833–2849.