Forest Ecology and Management 259 (2010) 2355–2365 Contents lists available at ScienceDirect Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco Assessing changes in forest fragmentation following infestation using time series Landsat imagery Nicholas C. Coops a,∗ , Steve N. Gillanders a,1 , Michael A. Wulder b,2 , Sarah E. Gergel c,3 , Trisalyn Nelson d,4 , Nicholas R. Goodwin e,5 a Department of Forest Resources Management, 2424 Main Mall, University of British Columbia, Vancouver V6T 1Z4, Canada Canadian Forest Service (Pacific Forestry Center), Natural Resources Canada, 506 West Burnside Road, Victoria V8Z 1M5, Canada c Department of Forest Sciences, 2424 Main Mall, University of British Columbia, Vancouver V6T 1Z4, Canada d Department of Geography, PO Box 3060 STN CSC, University of Victoria, Victoria, BC V8W 3R4, Canada e Queensland Government, Department of Natural Resources and Water, Climate Building, 80 Meiers Road, Indooroopilly, Qld 4068, Australia b a r t i c l e i n f o Article history: Received 7 January 2010 Received in revised form 9 March 2010 Accepted 12 March 2010 Keywords: Multi-temporal Time series Landsat Spatial pattern Landscape pattern index Mountain pine beetle Insect Fragmentation Connectivity Heterogeneity Disturbance a b s t r a c t The current epidemic of mountain pine beetle (Dendroctonus ponderosae Hopkins) in British Columbia, Canada, has impacted an area of over 13 million hectares presenting a considerable challenge to provincial forest resource managers. Remote sensing technologies offer a highly effective tool to monitor this impact due to very large areas involved and its ability to detect dead and dying tree crowns. Conventionally, change detection procedures based upon spectral values have been applied; however, analysis of landscape pattern changes associated with long-time series change detection approaches present opportunities for the generation of unique and ecologically important information. This study is focussed on the detection and monitoring of the shape and area characteristics of lodgepole pine stands during mountain pine beetle infestation to quantify the progression of forest fragmentation and related loss of landscape connectivity. A set of landscape pattern indices were applied to a set of images consisting of six Landsat satellite images spanning the period from 1993 to 2006. Our results indicate that the impacts of the mountain pine beetle infestation on forest spatial pattern consist of an increase in the number of patches, an increase in forest patch shape complexity, a reduction in forest patch size, an increase in forest patch isolation, and a decrease in interspersion. These findings demonstrate the unique information available from long-time series satellite imagery combined with pattern analysis to better understand the combined effects of insect infestation and forest salvage and harvesting. © 2010 Elsevier B.V. All rights reserved. 1. Introduction The role of disturbance in altering landscapes and modifying ecosystems over a range of spatial and temporal scales is well recognized in ecology (Perry, 2002; Drever et al., 2006; Noss and Lindenmayer, 2006; Jentsch, 2007). Pickett and White (1985) define disturbance as “any relatively discrete event in time that dis- ∗ Corresponding author. Tel.: +1 604 822 6452; fax: +1 604 822 9106. E-mail addresses: nicholas.coops@ubc.ca (N.C. Coops), sgilland@interchange.ubc.ca (S.N. Gillanders), mwulder@pfc.cfs.nrcan.gc.ca (M.A. Wulder), sgergel@interchange.ubc.ca (S.E. Gergel), tnelson@office.geog.uvic.ca (T. Nelson), Nicholas.Goodwin@nrw.qld.gov.au (N.R. Goodwin). 1 Tel.: +1 604 822 4148; fax: +1 604 822 9106. 2 Tel.: +1 250 363 6090; fax: +1 250 363 0775. 3 Tel.: +1 604 827 5163; fax: +1 604 822 9102. 4 Tel.: +1 250 472 5620. 5 Tel.: +61 07 3896 9650. 0378-1127/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.foreco.2010.03.008 rupts ecosystem, community, or population structure and changes resources, substrate availability, or the physical environment”. Disturbance events can contribute to the maintenance of biodiversity (Connell, 1978) and heterogeneity (Turner et al., 2003), as well as be the primary drivers of declines in biodiversity and species endangerment (Hansen et al., 2001). As a result, a greater understanding of the role of disturbance regimes promotes betterinformed management decisions, by gaining insight related to landscape dynamics and the historic range of variability in ecosystems. One of the key impacts of disturbance on forested landscapes is fragmentation where forested habitat is reduced into an increasing number of smaller, more isolated, patches (Wilcove et al., 1986). This can then result in a modification of the microclimate within and surrounding the remnant, intact forest patches (Saunders et al., 1991) and a change in forest ecosystem function and condition (Wickham et al., 2008). Natural disturbances tend to alter forest landscape pattern differently from anthropogenic impacts such as timber harvesting (Mladenoff et al., 1993). For instance, natural disturbances often 2356 N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 result in patches with less edge effect between patches as compared with timber harvesting (Tinker et al., 1998). Hudak et al. (2007) characterized a range of patch characteristics of stand replacing harvest and fire disturbances in a coniferous forest landscape, and found that clear-cutting practices resulted in smaller patch sizes, smaller patch perimeter lengths, greater inter-patch distances, more edge habitat, and less interior habitat when compared to landscape patterns created by natural disturbance processes. Similarly, Tinker et al. (1998) observed that timber harvesting fragmented the landscape through a distinct suite of structural changes including a decrease in forest patch core area and mean patch size, and an increase in edge density and patch density. Furthermore, timber harvesting tends to remove a larger amount of biomass from a forest than most natural disturbances, and results in the removal of those select stands with high timber volume and quality (Tinker et al., 1998). The role of non-stand replacing disturbances such as insect infestation on forest fragmentation, however, is not as well understood in forested landscapes. The current mountain pine beetle (Dendroctonus ponderosae Hopkins) epidemic in the central interior of British Columbia, Canada, represents a critical management and ecological concern. The epidemic is believed to have resulted from an absence of extreme winter temperatures, an abundance of suitable hosts, and a moderating trend in temperature extremes. As of January 2008, the cumulative area of British Columbia impacted (red or grey-attack) by the beetle was 13.5 million hectares (ha) (British Columbia Ministry of Forests and Range, 2008). This current area of infestation is beyond that of any previously recorded and has grown rapidly from a surveyed area of 164,000 ha in 1999 (Raffa et al., 2008). The extent and severity of this on-going infestation results in an expectation that this disturbance agent will have a significant impact on forest fragmentation and related connectivity (Gillanders et al., 2008). Prior to this outbreak-level mountain pine beetle infestation, large contiguous patches of forest were common over the landscape, characterized by smaller, isolated patches of natural disturbance including wildfire, windthrow, endemic insect infestation, as well as areas of anthropogenic disturbance such as roads and forestry activity. Functional connectivity may be compromised depending on the extent of these disturbances, however forests would be the dominant matrix component. As beetle populations increase, small isolated clusters of infested trees impacted by the beetle appear on the landscape. Over time, as beetle populations continue to expand, trees neighbouring these small patches of beetle-impacted timber were infested, and new isolated patches appear. In the case of the central interior of British Columbia, where lodgepole pine is a predominant coniferous species, it is expected that these patches of beetle-impacted trees will eventually coalesce and replace largely undisturbed forest as the dominant matrix component. Based on the severity of the current epidemic in British Columbia, this trend will continue until the beetle had largely exhausted the availability of its food source or was impacted by sufficiently cold and appropriately timed winter temperatures. However, in concert with the beetle infestation are anthropogenic disturbances including timber harvesting and salvage logging which consists of the removal of dead or dying trees in an effort to recover economic value that may otherwise be lost (Lindenmayer and Franklin, 2008). When the impacts of timber harvesting and salvage activities are combined with those of the current mountain pine beetle infestation, the cumulative impact of these landscape disturbances is expected to result in an increasingly fragmented landscape with a significantly reduced degree of connectivity (Gillanders et al., 2008). Due to the spatial extent of the current infestation, remotely sensed data provide opportunities to effectively monitor and evaluate the impacts of the mountain pine beetle over a range of scales (Wulder et al., 2005a). Typically, the mapping of mountain pine beetle using remote sensing technologies relies on the spectral response of vegetation indices (Price and Jakubauskas, 1998; Leckie et al., 2005; Skakun et al., 2003; Wulder et al., 2005b, 2006a), which detects changes in pixel-level vegetation conditions. In this research, we evaluate changes in select landscape pattern indices over forest stands that have undergone disturbance by both mountain pine beetle and timber harvest. Analysis of a historical archive of Landsat data at the forest stand level allowed for the development of temporal trajectories of a number of Landsat-derived landscape pattern indices, providing an ability to quantify the impact of fragmentation for both infestation and harvesting activity in order to characterize forest stand conditions during infestation. The specific objectives of this research are to (i) assess how the mountain pine beetle changes landscape spatial pattern and (ii) determine the relative impact of mountain pine beetle infestation and timber harvesting on forest fragmentation and connectivity. The results from this analysis will help to inform resource managers of the changes in landscape pattern occurring in the British Columbia interior as a result of mountain pine beetle infestation and timber harvesting, and will provide much needed information to help guide future forest and land base management decisions. 2. Methods 2.1. Study area The study area is situated in the Morice Timber Supply Area (TSA) which is part of the British Columbia Ministry of Forests and Range (MoFR) Nadina Forest District located in the central interior of British Columbia (Fig. 1). The Morice TSA is located on the western edge of the Central Interior Plateau and covers approximately 1.5 million ha (Tesera Systems Inc., 2006). The main forest species in the area include lodgepole pine (Pinus contorta var. latifolia), hybrid spruce (Picea engelmannii x glauca), and subalpine fir (Abies lasiocarpa). Trembling aspen (Populus tremuloides), silver (or amabilis) fir (Abies amabilis), western hemlock (Tsuga heterophylla), and mountain hemlock (Tsuga mertensiana) are also present in lesser amounts. In the north and central areas of the Morice TSA, mountain pine beetle infestation occurred in the mid-1990s while in the southern region the infestation occurred in the late 1990s (Nelson et al., 2006). The Morice TSA has undergone extensive harvesting and Fig. 1. Location of Morice Timber Supply Area (TSA) and Landsat path 51 row 22 in relation to the province of British Columbia, Canada. N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 2357 Table 1 Image dates and Landsat sensor types used for the data stack. Image acquisition date Landsat sensor 3 October, 1993 24 August, 1996 29 July, 1998 14 August, 2001 29 September, 2003 20 August, 2006 TM TM TM ETM+ TM TM management including mountain pine beetle treatment and salvaging efforts. Traditionally, subject to lower winter temperatures, the Morice TSA represented the northern extent for mountain pine beetle. 2.2. Data sources Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) satellite imagery, with 28.5 m pixels, were the primary data source for this study. The Landsat scene chosen was path 51/row 22, which covers the extent of the area of interest. Scenes were chosen which were largely cloud-free and within the seasonal window of July to September when available. The data stack consists of five Landsat-5 TM and one Landsat-7 ETM+ scenes covering the range of 1993–2006 (Table 1). 2.3. Ancillary data for stratification For stratification purposes, a 25-m digital elevation model (DEM) (BCMSRM, 2002) covering the extent of the study area was re-sampled to match the 28.5-m Landsat imagery and used to derive elevation values. In addition, vegetation resource inventory (VRI) data, derived from standard forest type mapping, were used to analyze forest cover attributes and enable stratification (British Columbia Ministry of Forests, 1998). In general, lowelevation stands are more likely to be infested by the beetle due to the warmer temperatures being more favourable for survival. Stands between 60 and 100 years are considered to be highly susceptible to mountain pine beetle (Carroll and Safranyik, 2004) and with crown closure ranges between 66 and 75% identified as highly susceptible to mountain pine beetle attack (Wulder et al., 2005a, 2006b). Areas within the Morice TSA which did not meet these above criteria were masked from further analysis. From the remaining area post-stratification, 23 subsets representative of the range of forest and fragmentation conditions were randomly selected. Analyzing a sample of the study area using landscape subsets allowed us to characterize trends associated with specific forest conditions and processes. Each subset was 200 × 200 pixels (5.7 km × 5.7 km or 3249 ha) (Fig. 2) in order to maximize the homogeneity of the landscape while still enabling a broad range of disturbance responses to be captured. This is double the minimum window size recommended by previous studies (Turner et al., 2001) enabling good spatial coverage of the study and representation of a range of forest conditions from heavily managed in the north to lightly managed in the central and southern portions of the TSA. 2.4. Pre-processing and classification Landsat data pre-processing included image-to-image geometric registration, radiometric calibration, and radiometric normalization. Each are critical processing steps required to ensure that detected land cover changes are not artefacts of atmospheric conditions, imaging and viewing conditions, sensor degradation, or pixel misalignment but actually represent changes to surface con- Fig. 2. Location of the twenty-three 200 × 200 pixel (3249 ha) image subset samples in relation to the Morice TSA boundary. Data were stratified by elevation and forest inventory attributes in order to collect replicates with similar initial conditions and mountain pine beetle susceptibility. ditions (Schott et al., 1988; Furby and Campbell, 2001; Roberts et al., 2002; Coops et al., 2006). Image registration was performed using a nearest-neighbour 2nd order polynomial transformation using the 2001 Landsat TM scene as the reference image, selected due to a lack of cloud cover. The remaining five images were co-registered with a root mean square error (RMSE) < 0.5 pixels (in both x and y directions). Following registration, any cloud and/or shadow found in the image stacks were manually removed via masking. Radiometric normalization was then applied, to account for differences in atmospheric conditions, solar angle, and satellite sensor characteristics. Radiometric normalization reduces the variability between images ensuring that detectable changes between image dates correspond to land cover changes rather than variability due to differences in satellite solar acquisition. To achieve this, the 2001 image date was atmospherically corrected using the COST (cosine of the solar zenith angle) model (Chavez, 1996). A relative normalization procedure was then applied to the remaining images in the stack using the Multivariate Alteration Detection (MAD) algorithm (Canty et al., 2004; Schroeder et al., 2006) which automatically normalizes multiple images via regression analysis based on stable targets in a base image. Following image pre-processing, the normalized difference moisture index (NDMI) was computed which has been shown to be sensitive to the foliage water content and fraction of dead leaf material (Hunt et al., 1987). The NDMI is the normalized ratio of the shortwave (TM band 5) and near infrared (TM band 4) bands. A low-pass filter (3 × 3 pixel window) was then applied to the NDMI images in order to minimize the influence of any pixel mis-registration. Data were then classified based on the spectral trajectory of the NDMI through time as developed and detailed by Goodwin et al. (2008). In brief the approach utilized a decision tree 2358 N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 Table 2 A selected listing and description of the landscape pattern indices (LPI) applied to the classified image subsets. LPI Description Category Reference Number of patches (NP) Total number of patches in a particular class or an entire landscape Total class area (ha) Patch area multiplied by proportional abundance of the patch (or patch type) (ha) Ratio of total edges (number of cells at patch boundary) and total area (total cells) (m/ha) Provides a measure of patch shape complexity by quantifying the mean fractal dimension of patches of the corresponding patch type, weighted by patch area The degree of isolation and fragmentation of the corresponding patch type The degree of aggregation or ‘clumpiness’ of a map based on adjacency of patches Area McGarigal and Marks (1995), Turner et al. (1989) McGarigal and Marks (1995) McGarigal and Marks (1995) Class area (CA) Area–weighted mean patch size (AREA AM) Edge density (ED) Area-weighted mean fractal dimension (FRAC AM) Mean proximity index (PROX MN) Interspersion and juxtaposition index (IJI) analysis which applied a series of thresholds and rules based on the index responses at field locations, through the time sequence of images, to define five classes, forest re-growth, harvest, healthy forest, mountain pine beetle infestation and no disturbance. The accuracy of discriminating beetle attack from healthy forest stands was assessed both spatially and temporally using 3 years of aerial survey data (1996, 2003, and 2004) with results indicating overall classification accuracies varying between 71 and 86% (Goodwin et al., 2008). Area Area Shape Shape Isolation/proximity Contagion/interspersion McGarigal and Marks (1995), Hargis et al. (1998), Li et al. (2005) McGarigal and Marks (1995) McGarigal and Marks (1995), Gustafson and Parker (1994) McGarigal and Marks (1995) 2.5. Landscape pattern indices Landscape pattern indices can be grouped into categories of area, shape, isolation/proximity, contagion/interspersion, and diversity (McGarigal and Marks, 1995). Based on a previous review, Gillanders et al. (2008) recommended a range of landscape pattern metrics be examined to determine the impact of infestation on the landscape. These included the number of patches, class area, and area-weighted mean patch area. In con- Fig. 3. Landsat TM and ETM+ image subset windows (RGB bands 543) displaying mountain pine beetle infestation (reddish-brown) and timber harvest (bright green = regenerating clearcut; magenta = recent clearcut). (a) and (b) land cover changes over a window representing a managed forest between 1996 and 2006, respectively; (c) and (d) land cover changes over a window representing an unmanaged forest between 1996 and 2006, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 2359 Fig. 4. Change in the landscape over the 13 years of Landsat imagery with respect to area of patches within the three classes: harvest (light green), forest (dark green) and mountain pine beetle (brown). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) cert with edge density, these indices provide an indication of the degree of fragmentation for the different land cover classes. Conversely, area-weighted mean fractal dimension, mean proximity index, and interspersion/juxtaposition index provide a means to characterize patch shape, patch isolation/proximity, and contagion/interspersion, respectively (Table 2). From the time series analysis of the Landsat imagery for each of the 23 subsets we have information regarding the area of each land cover type (forest, MPB infestation and harvest) at each of the six time steps. In addition we have 12 patch statistics averaged within each subset at each of the six time steps. Our analysis approach is as follows. First we investigate the overall changes in class area, and then average patch characteristics through time across the TSA. The other patch-level statistics allow analysis at the class level to assess the role that forest, timber harvest, and mountain pine beetle infestation play in contributing to overall fragmentation and loss of connectivity in the TSA. After this initial analysis, we apply factor analysis to reduce the variance in the 12 metrics and assess changes in the forest patches as a result of harvest and infestation. Analysis was undertaken in Fragstats (McGarigal and Marks, 1995) and Statistica (StatSoft Inc., 2005). 3. Results Evidence of changes to the landscape resulting from mountain pine beetle infestation and timber harvesting can be observed by viewing regions of the imagery where these processes are identified to have occurred. For example, Fig. 3(a) shows a managed landscape in 1996, with forest patches intersected by roads linking harvest clear-cuts, most of which are undergoing reforestation as indicated by the bright NIR reflectance. The adjacent image window shows the same landscape in 2006 (Fig. 3(b)) with recent harvest activ- ity as bright pink polygons as well as some new roads. Also visible in the 2006 image subset is the replacement of much of the formerly dark green forest by reddish-brown pixels, which represent the red or grey-attack stage of tree mortality caused by mountain pine beetle. Likewise, Fig. 3(c) shows an unmanaged landscape characterized by mostly contiguous forest in 1996. The 2006 image (Fig. 3(d)) shows evidence of mountain pine beetle infestation, particularly in the northeast of the subset, but also distributed as smaller patches amongst the otherwise green forest. The change in the landscape over the 13 years of Landsat imagery is shown diagrammatically in Fig. 4 with respect to area of patches within the three classes: harvest, forest, and infestation by mountain pine beetle. In 1993 the largest patches in the landscape were forest, with all patches greater than 1401 ha found in that class. Moderate patches between 700 and 1400 ha were harvested polygons; the smallest patches, with less than 700 ha, were either harvest or mountain pine beetle patches. In contrast, by 2001 mountain pine beetle-infested forest stands had increased in size up to 2100 ha with only the largest class having forested pixels. By 2006 no large forest polygons remain with mountain pine beetle the majority in the 1401–2100 ha class and the harvest class represented at a range of patch sizes from less than 700 to 2100 ha. Pattern analysis at the class level allows for an assessment of the role that forest, timber harvest, and mountain pine beetle infestation play in contributing to fragmentation and loss of connectivity. For instance, Fig. 5(a) shows the changes in class area (CA), (b) edge density (ED), (c) area-weighted mean fractal dimension (FRAC AM) and (d) standard deviation in fractal dimension (FRAC SD) for the forest, harvest, and mountain pine beetle infestation classes through the time series. In 1993 the forest class is clearly the dominant class with a value of 2282 ha compared to harvest with 542 ha and mountain pine beetle infestation with a 2360 N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 Fig. 5. Pattern analysis allows for an assessment of the role that forest, timber harvest, and mountain pine beetle infestation play in contributing to fragmentation and loss of connectivity; (a) shows the changes in class area (CA), (b) edge density (ED), (c) area-weighted mean fractal dimension (FRAC AM) and (d) standard deviation in fractal dimension (FRAC SD) for the forest, harvest, and mountain pine beetle infestation classes through the time series. value of 125 ha. While the harvest class tends to show a gradual increase through the time series with a value of 831 ha by 2006, the forest and mountain pine beetle infestation classes are much more dynamic. For instance, class area for the forest class steadily decreases to a value of 774 ha in 2006, while the mountain pine beetle infestation class increases to a maximum value of 1351 ha in 2006. The edge density values for the forest and the mountain pine beetle classes increase over the time period of the study, as expected from increasing fragmentation, while the harvest class edge density remains lower and constant throughout the 12-year period. The area-weighted mean fractal dimension plot shows a relatively stable trajectory for the harvest, an overall increase for the mountain pine beetle-infested class, and decrease in the N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 2361 Fig. 7. Individual patch loadings for each factor for the forested, harvested and mountain pine beetle-infested patches for the 13-year time period. Fig. 6. Factor loadings for the 12 patch metrics. forest class, while the standard deviation in fractal dimension remains constant for harvest and forest, however increases over the time period of mountain pine beetle-infested stands. Increasing FRAC AM, the area-weighted fractal dimension of all patches within that class, implies patches are becoming increasingly complex with larger perimeters, for a given area. As a result, as the beetle infestation increases, forest patches are being broken up into more complex shapes. Alternatively, the area of healthy forest is reducing, and the shapes of these remaining patches are not as complex. The variation in the fractal dimension is relatively constant across all patches in the case of harvested and forest patches, however the variation increases as the mountain pine beetle infestation increases over time and space. Factor analysis on the 12 fragmentation metrics indicates, as expected, the variation in the metrics can be reduced into a smaller number of overall factors. The first two factors explain 60% of the variation (38 and 22%, respectively). The weightings for the two factors is shown in Fig. 6 and indicates that ED, FRAC AM, and FRAC SD, all contribute positive loadings to the first factor, indicating that as the factor increases the patches have increasing numbers of edge, and increasing amount of edge per area and increasing variation in both configurations. The second factor has negative loadings of contributing area, interspersion and juxtaposition index (IJI) and positive loadings with the number of patches (Fig. 6). The factor weights generally correspond with the overall nature of the metrics as described in Table 2, with factor 1 generally driven by variation in shape-based patch metrics whereas factor 2 is broadly more aligned with area-based metrics. The individual patch loadings for each factor are shown in Fig. 7 for the forested, harvested, and mountain pine beetle-infested patches for the 13-year time period. The figure clearly shows the distinction in the patches based on the 2 major factors at the class level. Harvested patches have both negative factor 1 and 2 loadings with patches located exclusively in the lower right of the figure. This indicates harvested patches, regardless of their location within the Morice TSA study area, and the date of harvest, are generally consis- tent in shape and size with a small number of patches, and a simple shape as they tend to consist of linear patterns including rectangles and simple polygons. In contract, forest patches have positive factor 1 loadings indicating the shapes of these patches are more complex, with increase edge density and fractal dimensions. The variation in the second factor indicates that the number and area of the forest patches varies greatly over the region and over the 13 years of analysis. Finally the mountain pine beetle patches have mid-range factor 1, and higher positive factor 2 scores indicating many patches are similar in shape and area metrics to the forest parameters, however some of the patches also have unique dimensions not similar to either the forest or the harvesting patches. To investigate this further the patches are displayed in reference to the image data in Fig. 8. There appears to be no consistent patterns to the harvest patches from 1993 to 2006. The forest patches maintain the same shape parameters although are found to change in area properties (factor 2) with the number of patches increasing and area decreasing over time. Most noticeable however, is the trajectory of mountain pine beetle patches. Initially in the two Landsat scenes acquired in 1993 and 1996, the positive factor 2 scores indicates small patches, of small area, with simple shape metrics. As the infestation increases over the study area there is a clear alteration in the shape and area of these patches into fewer and larger patches with increasingly complex shape parameters. The trajectory of the mountain pine beetle patch characteristics from 1993 to 2006 form a trend that appears over time to be increasingly similar to the healthy forest patches in 1993. Summarising the fragmentation metrics up to the full landscape levels shows the overall impact of harvesting and mountain pine beetle on the landscape. As expected, Fig. 9(a) shows an overall increase in the number of patches, and a corresponding increase in the complexity of the patches from 1993 to 2006. Both trajectories clearly indicate the metrics change most quickly from 1993 to 2001 with the past 6 years being relatively constant with respect to the addition of new patches or changes in the existing patch characteristics. The edge density follows similar trajectories (Fig. 9(b)), which the interspersion and juxtaposition index has reduced in response to the increased number of edges and patches in the landscape. The decrease in interspersion/juxtaposition index values reflects 2362 N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 Fig. 8. Factor 1 and 2 loadings for the subsets displayed by image data from 1993 to 2006 by class. Fig. 9. Changes in patch metrics summarised over entire study area (a) number of patches and edge density follows similar trajectories, while (b) interspersion and juxtaposition index remains relatively constant and edge density increases. N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 2363 a trend in which patch types are becoming more disproportionally distributed or clumped. Contrary to what we expected, this trend represents a reduction in patch complexity. However, the stable period of the trajectories following 1998 suggests that while fragmentation continues to occur on the landscape, the relative distribution of classes remains constant. 4. Discussion Our purpose in analyzing these landscapes through time has been to determine how mountain pine beetle infestation changes landscape spatial pattern, and assess to what degree mountain pine beetle and timber harvesting contribute to forest fragmentation and loss of connectivity. By viewing the image subsets we can gain insight into how tree mortality, resulting from mountain pine beetle infestation, is manifested on the landscape. For example, Fig. 3 shows a landscape fragmented by timber harvest activities. By 2006, this landscape has been further impacted by the mountain pine beetle in addition to recent harvest activity. We can observe the patchy distribution of beetle-impacted stands and can see that what began as forest and harvest being the dominant matrix components in 1996, has changed to beetle-impacted forest and timber harvest; intact contiguous forest is now a minor component of the matrix. Likewise, we can observe a reduction in forest patch size, an increase in forest patch complexity, and an increase in the number of patches resulting from mountain pine beetle infestation. The results of the landscape pattern indices applied to the 23 image subsets suggest a trend toward a more fragmented landscape. This is consistent with what we would expect based on the considerable historic inventory of mature lodgepole pine in the study area, the presence of operational and salvage logging in the study area, and the severity of the current mountain pine beetle epidemic in the central interior of British Columbia. In general, the results of this analysis confirm that the impacts of the mountain pine beetle on forest spatial pattern in vulnerable areas of the Morice TSA study area consist of: 1. 2. 3. 4. 5. an increase in the number of patches; an increase in forest patch shape complexity; a reduction in forest patch size; an increase in forest patch isolation; and a decrease in interspersion. The results provide detailed information related to the relative contribution that mountain pine beetle infestation and forest harvest makes to overall landscape fragmentation and patch connectivity. Based on our results, in regions identified as vulnerable to mountain pine beetle infestation, forest harvest plays a minor role when compared to the impacts of the mountain pine beetle. For instance, the number of patches and area-weighted mean patch area indicate that the landscape has become more fragmented through the time series due to a decrease in forest patch size and an increase in the number of forest patches. This observation is further supported by the trends of the mountain pine beetle-infested class, which shows a similar trend in magnitude but results in an inverse trajectory when compared to the forest class. After 2000, the number of mountain pine beetle-infested patches declines suggesting that patches of beetle-impacted forest are coalescing, resulting in a more contiguous distribution as the mountain pine beetle-infested class becomes the dominant matrix component. In contrast, the timber harvest class remains relatively stable both in terms of the number of patches and patch size. The class level plot for edge density (Fig. 9(b)) shows a relatively low and stable trajectory for timber harvest when compared to the forest and mountain pine beetle-infested classes. The increase in Fig. 10. Number of patches through time confirms a threshold between 2001 and 2003 at which the mountain pine beetle class trend line crosses the forest trend line. This represents the point at which the mountain pine beetle-infested class becomes the dominant matrix component within the image subset samples. edge for both the forest and the mountain pine beetle-infested classes is consistent with the impacts of mountain pine beetle infestation on a pine forest; namely, mostly contiguous forest being fragmented by small isolated patches of beetle-impacted forest, which continue to multiply and spread across the landscape. Although the trajectories for the forest and mountain pine beetle classes follow an almost identical trend, a threshold is crossed at which point the beetle-impacted forest replaces the forest as the major landscape matrix component. This is represented as the levelling off, and is followed by the decline which suggests that these small patches of beetle-impacted forest are merging, resulting in patches with a greater core/edge ratio. The trajectories of these metrics are similar to what we had predicted in a previous study (Gillanders et al., 2008), although we anticipated that timber harvest would display a considerable increase through the time series for both landscape pattern indices. This discrepancy can in part be explained by how the classifier merges adjacent polygons, thereby exhibiting an increase in patch area but not an increase in number of patches. However, even with the use of a classifier that treated harvest events as distinct objects, the trajectories for the timber harvest class would have remained relatively static compared to those of the forest and mountain pine beetle-infested classes. The analysis of the number of patches of the dominant class through time confirms a threshold between 2001 and 2003 at which the mountain pine beetle class trend line crosses the forest trend line (Fig. 10). This represents the point at which the mountain pine beetle-infested class becomes the dominant matrix component within the image subset samples. While this phenomenon is an indicator of the severe impacts of mountain pine beetle infestation, the inverse trajectory could serve as an indicator of forest regeneration and recovery from disturbance. While both timber harvesting and the impacts from mountain pine beetle infestation contribute to forest fragmentation and loss of connectivity, each represent different disturbance agents with different impacts. Although the results of this study clearly indicate that the mountain pine beetle has a greater relative impact on forest fragmentation and loss of connectivity than timber harvest in the vulnerable areas, mountain pine beetle infestation repre- 2364 N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 sents a natural disturbance and a natural fragmentation event. Thus although the impacts of mountain pine beetle infestation to spatial pattern may be greater than timber harvesting in the study area, the implications to biotic and abiotic processes are markedly different. While the results of these analyses provide a general indication of changes to the ecological conditions of the landscape due to changes in spatial pattern, we recognize the inherent limitations of the spatial resolution of a data source such as Landsat satellite imagery. The use of these data do not allow for an investigation of patch-level dynamics, which are expected to be very different for natural and anthropogenic disturbances. For instance, tree mortality caused by mountain pine beetle infestation leaves an abundance of snags and coarse woody debris, which represent valuable habitat for cavity nesters and a range of forest vertebrates (Chan-McLeod, 2006), while timber harvesting does not. The classification accuracy of the thresholding technique on the time series of Landsat data was generally very good, however Goodwin et al. (2008) found that evidence of spectral confusion exists, especially for forest stands with low levels of infestation. For example, the 1996 aerial survey locations had a lower accuracy (49%, relative to the >78% accuracy for the 2003 and 2004 surveys), consistent with previous research which has demonstrated that stands with small numbers of infested trees are more difficult to classify than sites with over 30 attack trees in each Landsat pixel (Skakun et al., 2003). Furthermore, the complexity and differences in edge effects between a natural and anthropogenic disturbance event are generalized, preventing the observation of patch-level differences. 5. Conclusions This research represents a methodological approach to monitoring disturbance by applying pattern indices to classified maps derived from multi-temporal Landsat imagery. Pattern indices quantify landscape structure and provide a means to infer ecological conditions. In this case, the use of multi-temporal data provided a historical perspective and contributed to our understanding of how the mountain pine beetle interacted with the landscape of the study area over a period of 13 years, beginning with low levels of infestation (pre-outbreak) progressing to epidemic levels (outbreak) where mountain pine beetle-impacted forest had become the dominant matrix component. The outcome of this research reveals that the study area has undergone considerable land cover changes. The impacts of the mountain pine beetle on landscape spatial pattern are notable in terms of both the immediate impacts to biota and ecological processes in the region but also the legacy that these changes will have in forming future forest composition and structure. The implications of these structural changes to ecological processes are as yet unclear, but changes in structural diversity would be expected to result in changes to both compositional and functional diversity (Mladenoff et al., 1993). The potential for landscape pattern indices to be applied to a monitoring program to evaluate land cover dynamics for a range of natural and anthropogenic disturbances has been demonstrated. These methods can easily be applied to monitor other landscapescale forest disturbances such as wildfire, other forest insects, and disease to allow for the characterization of the disturbance as well as to provide an indication of its ecological implications. This utilization of multi-temporal Landsat data to monitor disturbance dynamics via changes in spatial pattern further supports the wide-ranging applicability of the Landsat suite of sensors for monitoring the Earth’s surface. With the recent announcement by the USGS (United States Geological Service) of unrestricted global access to the entire Landsat data archive (Woodcock et al., 2008), it is expected that research utilizing multiple satellite images will continue to proliferate. In addition to the utility and availability of Landsat data, this research highlights the data-rich nature of satellite imagery in general; namely, the value of assessing spatial pattern rather than solely relying on spectral information. The assessment of the distribution and arrangement of patches on a given landscape provides opportunities to infer broad-scale ecological conditions, and when used in a multi-temporal setting, allows for the detection of changes that can affect a wide range of both biota and ecological processes. Further work is required to determine the degree to which our results are characteristic of other landscapes undergoing epidemic mountain pine beetle infestation. This would require applying these methods to other forested landscapes in the British Columbia interior undergoing beetle infestation in order to determine the natural variability inherent in both the landscape and the nature of the beetle infestation. Acknowledgements This project is funded by the Government of Canada through the Mountain Pine Beetle Program administered by Natural Resources Canada, Canadian Forest Service, with additional information available at: http://mpb.cfs.nrcan.gc.ca. The Landsat data used in this study was contributed by the U.S. Geological Survey Landsat Data Continuity Mission Project through participation of Wulder on the Landsat Science Team. References British Columbia Ministry of Forests, 1998. Vegetation Resources Inventory: Photo Interpretation Standards, Version 1.0. Resources Inventory Committee, Victoria, British Columbia, Canada, 65 pp. British Columbia Ministry of Forests and Range, 2008. Beetle Facts. URL: http://www.for.gov.bc.ca/hfp/mountain pine beetle/facts.htm (date accessed 9 February 2009). British Columbia Ministry of Sustainable Resource Management, 2002. Gridded Digital Elevation Model Product Specifications, Edition 2.0. December 2002, Base Mapping and Geomatics Services Branch, Victoria, British Columbia, Canada, 23 pp. Canty, M.J., Nielsen, A.A., Schmidt, M., 2004. Automatic radiometric normalization of multitemporal satellite imagery. Remote Sensing of Environment 91, 441–451. Carroll, A.L., Safranyik, L., 2004. The bionomics of the mountain pine beetle in lodgepole pine forests: establishing a context. In: Shore, T.L., Brooks, J.E., Stone, J.E. (Eds.), Mountain Pine Beetle Symposium: Challenges and Solutions. October 30–31, 2003, Kelowna, British Columbia. Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Information Report BC-X-399, Victoria, BC, pp. 21–32, 298. Chan-McLeod, A.C.A., 2006. A review and synthesis of the effects of unsalvaged mountain-pine-beetle-attacked stands on wildlife and implications for forest management. BC Journal of Ecosystems and Management 7, 119–132. Chavez, P.S., 1996. Image-based atmospheric corrections—revisited and improved. Photogrammetric Engineering and Remote Sensing 62, 1025–1036. Connell, J.H., 1978. Diversity in tropical rain forests and coral reefs. Science 199, 1302–1310. Coops, N.C., Wulder, M.A., White, J.C., 2006. Identifying and describing forest disturbance and spatial pattern: data selection issues and methodological implications. In: Franklin, S.E., Wulder, M.A. (Eds.), Understanding Forest Disturbance and Spatial Pattern: Remote Sensing and GIS Approaches. Taylor and Francis, Boca Raton, FL, USA, pp. 31–62, 246. Drever, C.R., Peterson, G., Messier, C., Bergeron, Y., Flannigan, M., 2006. Can forest management based on natural disturbances maintain ecological resilience? Canadian Journal of Forest Research 36, 2285–2299. Furby, S.L., Campbell, N.A., 2001. Calibrating images from different dates to ‘like value’ digital counts. Remote Sensing of Environment 77, 186–196. Gillanders, S.N., Coops, N.C., Wulder, M.A., Gergel, S.E., Nelson, T., 2008. Multitemporal remote sensing of landscape dynamics and pattern change: describing natural and anthropogenic trends. Progress in Physical Geography 32, 503–528. Goodwin, N.R., Coops, N.C., Wulder, M.A., Gillanders, S., Schroeder, T.A., Nelson, T., 2008. Estimation of insect infestation dynamics using a temporal sequence of Landsat data. Remote Sensing of Environment 112, 3680–3689. Gustafson, E.J., Parker, G.R., 1994. Using an index of habitat patch proximity for landscape design. Landscape and Urban Planning 29, 117–130. Hansen, A.J., Neilson, R.P., Dale, V.H., Flather, C.H., Iverson, L.R., Currie, D.J., Shafer, S., Cook, R., Bartlein, P., 2001. Global change in forests: responses of species, communities, and biomes. BioScience 51, 765–779. Hargis, C.D., Bissonette, J.A., David, J.L., 1998. The behavior of landscape metrics commonly used in the study of habitat fragmentation. Landscape Ecology 13, 167–186. N.C. Coops et al. / Forest Ecology and Management 259 (2010) 2355–2365 Hudak, A.T., Morgan, P., Bobbitt, M., Lentile, L., 2007. Characterizing stand-replacing harvest and fire disturbance patches in a forested landscape: a case study from Cooney Ridge Montana. In: Wulder, M.A., Franklin, S.E. (Eds.), Understanding Forest Disturbance and Spatial Pattern: Remote Sensing and GIS Approaches. CRC Press (Taylor and Francis), Boca Raton, FL, pp. 209–231, 246. Hunt, E.R., Rock, B.N., Nobel, P.S., 1987. Measurement of leaf relative water content by infrared reflectance. Remote Sensing of Environment 22, 429–435. Jentsch, A., 2007. The challenge to restore processes in face of nonlinear dynamics—on the crucial role of disturbance regimes. Restoration Ecology 15, 334–339. Leckie, D.G., Cloney, E., Joyce, S.P., 2005. Automated detection and mapping of crown discolouration caused by jack pine budworm with 2.5 m resolution multispectral imagery. International Journal of Applied Earth Observation and Geoinformation 7, 61–77. Li, X., He, H.S., Bu, R., Wen, Q., Chang, Y., Hu, Y., Li, Y., 2005. The adequacy of different landscape metrics for various landscape patterns. Pattern Recognition 38, 2626–2638. Lindenmayer, D.B., Franklin, J.F., 2008. Salvage Logging and its Ecological Consequences. Island Press, Washington, DC, USA, 227. McGarigal, K., Marks, B.J., 1995. FRAGSTATS: Spatial Pattern Analysis Program for Quantifying Landscape Structure. U.S. Department of Agriculture Forest Service. Pacific Northwest Research Station, Portland, OR, p. 122. Report nr Gen. Tech. Re. PNW-GTR-351. Mladenoff, D.J., White, M.A., Pastor, J., Crow, T.R., 1993. Comparing spatial pattern in unaltered old-growth and disturbed forest landscapes. Ecological Applications 3, 294–306. Nelson, T., Boots, B., White, K.J., Smith, A.A.C., 2006. The impact of treatment on mountain pine beetle infestation rates. BC Journal of Ecosystems and Management 7, 20–36. Noss, R.F., Lindenmayer, D.B., 2006. Special section: the ecological effects of salvage logging after natural disturbance. Conservation Biology 20, 946–948. Perry, G.L.W., 2002. Landscapes, space and equilibrium: shifting viewpoints. Progress in Physical Geography 26, 339–359. Pickett, S.T., White, P.S., 1985. The Ecology of Natural Disturbance and Patch Dynamics. Academic Press, New York, NY, USA, p. 472. Price, K.P., Jakubauskas, M.E., 1998. Spectral retrogression and insect damage in lodgepole pine forests. International Journal of Remote Sensing 19 (8), 1627–1632. Raffa, K.F., Aukema, B.H., Bentz, B.J., Carroll, A.L., Hicke, J.A., Turner, M.G., Romme, W.H., 2008. Cross-scale drivers of natural disturbance prone to anthropogenic amplification: the dynamics of bark beetle eruptions. BioScience 58 (6), 501–517. Roberts, D.A., Numata, I., Holmes, K.W., Batista, G., Krug, T., Monteiro, A., Powell, B., Chadwick, O., 2002. Large area mapping of land-cover change in Rondônia using multitemporal spectral mixture analysis and decision tree classifiers. Journal of Geophysical Research Atmospheres 107 (D20), 8073, LBA 40-1 to 40-18. 2365 Saunders, D.A., Hobbs, R.J., Margules, C.R., 1991. Biological consequences of ecosystem fragmentation: a review. Conservation Biology 5, 18–32. Schroeder, T.A., Cohen, W.B., Song, C., Canty, M.J., Yang, Z., 2006. Radiometric correction of multi-temporal Landsat data for characterization of early successional forest patterns in western Oregon. Remote Sensing of Environment 103, 16–26. Schott, J.R., Salvaggio, C., Volchok, W.J., 1988. Radiometric scene normalization using pseudoinvariant features. Remote Sensing of Environment 26, 1–16. Skakun, R.S., Wulder, M.A., Franklin, S.E., 2003. Sensitivity of the thematic mapper enhanced wetness difference index to detect mountain pine beetle red-attack damage. Remote Sensing of Environment 86, 433–443. StatSoft, Inc., (2005). STATISTICA (data analysis software system), version 7.1. www.statsoft.com. Tulsa OK. USA. Tesera Systems Inc., 2006. Morice TSA Timber Supply Analysis Report for the Current Status, Decision & Mitigation Scenarios. Morice and Lakes timber supply area innovative forest practices agreement, Cochrane, Alberta, p. 216. Tinker, D.B., Resor, C.A.C., Beauvais, G.P., Kipfmueller, K.F., Fernandes, C.I., Baker, W.L., 1998. Watershed analysis of forest fragmentation by clearcuts and roads in a Wyoming forest. Landscape Ecology 13, 149–165. Turner, M.G., O’Neill, R.V., Gardner, R.H., Milne, B.T., 1989. Effects of changing spatial scale on the analysis of landscape pattern. Landscape Ecology 3, 153–162. Turner, M.G., Collins, S.L., Lugo, A.L., Magnuson, J.J., Rupp, T.S., Swanson, F.J., 2003. Disturbance dynamics and ecological response: the contribution of long-term ecological research. Bioscience 53, 46–56. Turner, M.G., Gardner, R.H., O’Neill, R.V., 2001. Landscape Ecology in Theory and Practice. Springer, New York, 401 pp. Wickham, J., Riitters, K., Wade, T., Homer, C., 2008. Temporal change in fragmentation of continental US forests. Landscape Ecology 23, 891–898. Wilcove, D.S., McLellan, C.H., Dobson, A.P., 1986. Habitat fragmentation in the temperate zone. In: Soulé, M. (Ed.), Conservation Biology: Science of Scarcity and Diversity. Sinauer Associates, Sunderland, MA, pp. 237–256. Woodcock, C.E., Allen, R., Anderson, M., Belward, A., Bindschadler, R., Cohen, W., Gao, F., Goward, S.N., Helder, D., Helmer, E., Nemani, R., Oreopoulos, L., Schott, J., Thenkabail, P.S., Vermote, E.F., Vogelmann, J., Wulder, M.A., Wynne, R., 2008. Free access to landsat imagery. Science 320, 1011–11011. Wulder, M.A., Dymond, C.C., White, J.C., Leckie, D.G., Carroll, A.L., 2005a. Surveying mountain pine beetle damage of forests: a review of remote sensing opportunities. Forest Ecology and Management 221, 27–41. Wulder, M.A., White, J.C., Coops, N.C., Han, T., Alvarez, M.F., 2005b. A Protocol For Detecting And Mapping Mountain Pine Beetle Damage From A Time-series of Landsat TM or ETM+ Data. Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Victoria, British Columbia, p. 70. Report nr Version 1.0. Wulder, M.A., White, J.C., Bentz, B.J., Ebata, T., 2006a. Augmenting the existing survey hierarchy for mountain pine beetle red-attack damage with satellite remotely sensed data. The Forestry Chronicle 82, 187–202. Wulder, M.A., White, J.C., Bentz, B., Alvarez, M.F., Coops, N.C., 2006b. Estimating the probability of mountain pine beetle red-attack damage. Remote Sensing of Environment 101, 150–166.