GEOG 566

         Advanced spatial statistics and GIScience

June 11, 2017

Beetle outbreaks and wildfire: overlapping in time and space

Filed under: Final Project @ 10:53 pm


Often landscape scale disturbances are singular events, but contemporary disturbance regimes are demonstrating that large-scale disturbances can overlap in time and space.  The spatial configuration an initial disturbance may influence the process and subsequent spatial configuration of a secondary disturbance, such as beetle outbreak followed by fire.  The cumulative effects of sequential disturbances could influence forest structure and successional trajectories through the next century for the landscapes in central interior British Columbia affect by beetle outbreaks and wildfire.

Figure 1: Post-fire photo of a fire that burned through beetle-killed forest in 2014 with various patches of fire severity.

Research Question

How is the landscape variability in fire disturbance related to the landscape variability in antecedent beetle outbreak disturbance?

Data Set

Landsat imagery was used for the analysis. Landsat data is 30-meter resolution multispectral data.  A single Landsat scene contained the wildfire event of interested (Row/Path 50/23).  Pre- and post-fire images were acquired with minimal cloud coverage and anniversary dates to match phenology.  These images were used to generate a burn severity map through the differenced Normalized Burn Ratio (dNBR) (Key and Benson 2006). The pre-fire NBR image was used to characterize the beetle outbreak on the pre-fire landscape. Pixel indexing formulas:


dNBR = NBRprefire – NBRpostfire


My ecological null hypothesis is that burning conditions (i.e. weather during day of burn) exert strong control over patch size.  Fire conducive weather, dry and windy conditions, control the disturbance variability across the landscape by facilitating crown fire in lodgepole pine dominated landscapes.  The ecological alternative hypothesis is that beetle outbreaks change the landscape structure by killing trees, which results in decreased continuity of canopy fuels. The lack of canopy fuels may constrain fire spread to the surface and may influence the disturbance variability of fire across that landscape.



A semivariogram was used to assess the spatial variability between pixel values in the dNBR image. A semivariogram is a means of evaluating spatial correlation between locations, and the shape of the semivariogram may indicate the level of patchiness and average patch size. The semivariogram analysis was conducted in R version 3.3.2 (Team 2016) with the USDM Package (Naimi 2015).


The semivariogram analysis was followed by a cross-variogram to assess the spatial relationship between fire and beetle outbreak patterns. The cross-variogram analysis was conducted with the gstat package in R version 3.3.2 (Pebesma and Graeler 2017).  It should be noted that some additional data extraction steps were needed to conduct the cross variogram.  The Gstat package does not work with raster data; therefore, a csv file was created with the cell value and UTM coordinates for each cell in the raster. The analysis was conducted in R version 3.3.2 (Team 2016).

Patch Metrics (Fragstats – Species Distribution Models)

Patch metrics were calculated with with the SDMtools Package (VanDerWal 2015) in R version 3.3.2 (Team 2016).   The SDMtools package is based on the work by McGarigal et al. (2012) FRAGSTATS software package and calculates many of the same patch metrics. In order to calculate patch metrics, I converted the continuous raster into a classified raster.  The burn severity classification was based on break points in the dNBR values defined by Reilly et al. (2017), and the beetle severity classification was based on cross-validation with orthorectified aerial photographs.  Burn classes included a low (plus unburned), moderate, and high, and beetle class included none, low-moderate, and moderate-high. Patch metrics were calculated for each disturbance severity class.

Geographically Weighted Regression

 Geographically Weighted Regression (GWR) was used as a descriptive analysis to understand the relationship between the variability in beetle outbreak and wildfire on the landscape.  The purpose of GWR is to determine areas of non-stationarity, or inconsistencies in the spatial structure over the landscape. I used the continuous raster layers for beetle outbreak and burn severity.  While there are packages in R to conduct GWR, I ended up conducting the analysis in ArcGIS (ESRI 2011).  I had difficulty trouble shooting the R packages. The raster layers were converted to point layers for the analysis.

 Spatial Overlay Analysis

 Spatial overlay analysis has been used in previous studies examining disturbances that overlap in time and space (Hart et al. 2015).  It is simply overlaying maps for each disturbance and determining where the disturbances overlap. In this case, I am looking for overlap between beetle outbreak and wildfire.  Both disturbance maps were classified as binary either burned or no burn and beetle outbreak or no beetle outbreak. Classification was based cross-referencing location with higher resolution orthorectified aerial photographs or higher resolution satellite imagery. I overlaid the two disturbance maps and multiplied them together to end up with four distinct classes: No Beetle – No Burn, No Beetle – Burn, Beetle – No Burn, and Beetle – Burn. Analysis was conducted in R version 3.3.2 (Team 2016).



The semivariogram analysis indicated a spherical relationship in the dNBR image, which suggests some level of patchiness (Figure 2).  The sill of the semivariogram suggests that patches averaged about 1500m, which seemed a bit high (Figure 2).  The dNBR image contained a couple large bodies of water, which may influence the semivariogram.  I removed the water bodies from the dNBR image, and ran the semivariogram on the water free image (Figure 2).  The water free analysis shows a shift in the sill location and indicates that average patch size might be around 500m (Figure 2).  The semivariogram without water also seems to have a double sill, which might indicate that there are both smaller and larger patches.


Figure 2: Semivariogram for the Entiako dNBR image with water included (a) and excluded (b). The red circles indicate the sill, which provide a general indication of patch size. With increasing distance (Lag) there is decreasing variability. The dNBR images (c and d)


The semivariogram analysis indicated a spherical relationship in the wildfire pattern and the beetle outbreak pattern, which suggests some level of patchiness (Figure 3).  The sill of the semivariogram suggests that patches averaged about 500-750m for the wildfire patterns and 250-500m for beetle outbreak patterns (Red circle in figure 3).  The cross-variogram in Figure 4c shows an undulating pattern, which indicates auto-correlation at multiple lag distances. Auto-correlation at multiple lag distances suggests that there are similar patches of fire severity that had similar patches of beetle outbreak severity, which may support this idea of the pattern of beetle outbreak influencing the pattern of wildfire.

The cross-variogram illustrates auto-correlation at multiple lag distances with patches of multiple sizes or patches that are regularly spaced. The low cross-variance at lag 300 (the nugget or y-intercept) might indicate that low/high burn severity is associated with areas of low/high beetle outbreak, which seems like it would correspond well to certain areas of the maps. The cross-variogram shows a second dip at a lag distance of about 2000m, which could describe patterns of low/high burn severity that are 2000m from low/high beetle outbreak.  The landscape patterns shown in the maps could also support this relationship.

Figure 3: Semivariogram for the Entiako wildfire (A) and beetle outbreak (B). The red circles indicate the sill, which provide a general indication of patch size. With increasing distance (Lag) there is increasing variability. The dNBR image (C) and NBR image (D)

Figure 4: Cross-variogram analysis for beetle outbreak and wildfire overlap. Graph a is a semivariogram for the wildfire raster. Graph b is a variogram for the beetle outbreak raster. Graph c is the cross-variogram analysis that overlays the raster used in graph a with the raster used in graph b.

Patch Metrics (Fragstats with Species Distribution Models)

In the SDMtools package, I ran two analyses, one on the whole landscape and one on each disturbance severity class. For the first analysis, I found that high burn severity was the smallest component for this fire event only accounting for about 15 percent of the landscape (see Table 1). This was interesting, because lodgepole pine is often characterized by ‘stand replacing’ high severity fire, however, moderate burn severity can also account for ‘stand replacing’ fire in lodgepole pine dominated landscapes, which was about 45 percent of the landscape (see Table 1). The landscape was dominated by low-moderate and moderate-high beetle severity with only a small portion unaffected by beetle outbreak (see Table 2).

Table 1: Burn severity class summary for patches, area, and proportion of landscape

Burn Severity Class Number of Patches Total Area in Hectares Proportion of Landscape
Low 382 2767.14 0.394
Moderate 448 3179.79 0.452
High 422 1082.70 0.154

 Table 2: Beetle severity class summary for patches, area, and proportion of landscape

Beetle Severity Class Number of Patches Total Area in Hectares Proportion of Landscape
None 985 780.93 0.111
Low-Moderate 656 2968.56 0.423
Moderate-High 709 3273.30 0.466

 The second analysis examined the patch metrics for each disturbance severity class. Table 3 and 4 includes patch area metrics by burn severity and beetle severity respectively. While the output includes many other variables, I wanted to look at the variability in patch size for disturbance type and class. The high burn severity maximum and mean patch size are the most interesting part of Table 3, because they are quite different from the low and moderate burn severity.  Further statistical analysis is needed to determine if they numbers are statistically different.  The summary of patch size metrics for both disturbance types indicate that there are many small patches for all class types.

Table 3: Burn severity classes patch size summary

Burn Class Number of patches Minimum Patch size Maximum Patch size Mean Patch Size Median Patch size SD
Low 382 0.09 1449.45 7.24 0.18 77.7
Moderate 448 0.09 1461.69 7.10 0.18 75.3
High 422 0.09 445.95 2.57 0.18 22.5

 Table 4: Beetle severity classes patch size summary

Burn Class Number of patches Minimum Patch size Maximum Patch size Mean Patch Size Median Patch size SD
None 985 0.09 77.58 0.79 0.18 4.04
Low-Moderate 656 0.09 2380.86 4.53 0.18 93.07
Moderate-High 709 0.09 1461.60 4.62 0.18 57.68

The SDM tools packages offers much more then is being demonstrated here. The Patchstat function can provides the following outputs: number of cells for each patch, number of cells in the core area, number of outer perimeter cell edges of the patch, number of internal cell edges of the patch, the area of each patch within the landscape mosaic, the interior core area of a patch, perimeter of patch, perimeter area ratio, shape index, fractal dimension index, and the core area index. The ClassStat Function also offers 28 metrics that are included in the output.  The documentation for the SDMTools package offers documentation on all of these metrics (VanDerWal 2015).

Geographically Weighted Regression

The GWR analysis is visualized (Figure 5) with the coefficients.  Reds and oranges indicate areas of positive association while blues indicate negative associations.  These associations visually make sense when going back and looking at the beetle outbreak map and burn severity map.  High beetle outbreak would have lower raster values, while high fire severity would have higher raster values. While the GWR seems to indicate some interesting patterns, I explored some of the variability with scatter plots, which is discussed below.

Figure 5: Geographically weighted regression showing the coefficients from the output. Positive values indicate a positive association, while negative values indicate a negative association.

Spatial Overlay Analysis

The spatial overlay analysis is shown in Figure 6.  Most of the landscape is covered by red, which is a combination of beetle outbreak and wildfire. There are a few small patches of live trees that were unaffected by the beetle outbreak that burned and are shown in orange. A number of these live patches that burned seemed to have burned at high severity.  There are number of areas that have beetle outbreak but did not experience fire, shown in green.  While areas with no beetle and no fire are very few and appear to have relatively small areas when looking at the map.

Figure 6: Spatial Overlay Analysis using binary disturbance classifications to determine where beetle outbreak and fire overlap in the landscape.

Inspection of GWR with spatial overlay

The final piece of this project was revisiting all the maps and thinking critically of the GWR Analysis.  Scatter plots where generated from six selected spots at the extreme ends of the GWR coefficient range in Figure 8 and 10. NBR and dNBR values were extracted and plotted in Figure 10. The purpose of generating scatter plots was to assess the relationship between these two variables in relation to the GWR coefficients. Maps for beetle outbreak (Figure 7), fire disturbance (Figure 8), GWR (Figure9), and Spatial overlay analysis (Figure 9) are shown below with rectangles outlining were values were extracted from for the scatter plots in Figure 10.

Plot 3 and 8 showed a strong linear relationship (see Figure 10).  The other plots visually did not show strong relationships.  The strong and weak linear relationships are supported by the R-squared values for those areas, where Plot 3 ranged from 0.59-0.83, Plot 8 ranged from 0.29-0.55, and Plot 4 ranged from .11-.25.  In looking at the scatter plots, the values associated with the spatial overlay analysis, beetle outbreak and fire severity have quite a range of variability across the scatter plots.

Fire is a complex process and this project has been primarily examining it through a lens of the pre-fire landscape, which overlooks weather and topography as controls of fire severity. Wildfire is a dynamic process that is influenced by fuels, weather, and topography on both macro- and micro-scales. The focus here has been solely on fuels, which likely to does not capture all the landscape variability generated by the wildfire event. Future analyses will account for topography and fire weather.  Additionally, since the focus is on forested landscapes and the interaction of overlapping disturbances in forested landscapes, future data sets may exclude non-forested areas.

Figure 7: The landscape distribution of Beetle outbreak disturbance for Entiako burn perimeter and areas associated with scatter plots in Figure 10.


Figure 8: The landscape distribution of burn severity for the Entiako fire that burned in 2012 and areas associated with scatter plots in Figure 10.

Figure 8: The landscape distribution of coefficients from the geographically weighted regression analysis and areas associated with scatter plots in Figure 10.


Figure 9: The landscape distribution of spatial overlay analysis with areas associated with scatter plot in Figure 10.


Figure 10: Scatter plots of NBR values for Beetle outbreak on the x-axis and dNBR values for burn severity on the y-axis. Point colors correspond to the spatial overlay analysis.


Overlapping disturbances of beetle-outbreak and wildfire will influence landscape structure and function for decades in central interior British Columbia, Canada.  It is important to characterize the distribution of disturbance severity across a landscape for both individual disturbances and the cumulative effects of overlapping disturbances. There is still much to explore within the analyses here and further analyses that include additional fires as well as topography and weather as explanatory variables of wildfire processes.

The spatial overlay analysis is visually the most dramatic showing much of the landscape dominated by the cumulative effects of overlapping disturbances, but also additional spots of fire mortality. This is important for future forest structure and successional trajectories.  It will be important for scientists and managers to monitor the landscapes response to the cumulative disturbance effect.  Both disturbances affect forest that offer late seral forest habitat utilized in the winter by caribou. The fire in this project plus an additional fire have affected most of the wintering habitat for this caribou herd, which is of great concern for managers.

Wildfire following beetle outbreak may facilitate or impeded lodgepole pine regeneration by preparing the surface for seedling establishment and releasing seeds or by killing or consuming seeds during the fire process. It will be important for scientist and managers to monitor regeneration on this landscape to understand forest response to overlapping disturbances.  This landscape is managed as wilderness, so it provides insight to natural processes unfolding with minimal human intervention under current climate conditions, which may be useful for predicting how surround areas may respond under similar disturbance interactions.

Software learning

Time was invested in learning many of the spatial analysis I performed in R. I thought that this might be a useful invest of time and learning.  R offers a number of spatial tools similar to ARC GIS.  I think R has a number of benefits over ARC GIS, in particular, the ability to track the workflow that could then be recycled for future analyses.  Additionally, once the code is setup, the process is semi-automated.

I was working with a particularly large fire, tens of thousands of pixels.  I found R to process the data quickly. I used only one fire event for this project, but have two other fires I want to look at, which I can now easily run through the r-scripts I have generated.

During this class, I also learned how to use R-Markdown.  R-Markdown offers an interface for coding and commenting, essentially a scientific notebook. R-Markdown allows for easy documentation of code chunks, and the ability to generate PDFs or HTML files for sharing. This documentation will allow me to more easily go back to the code later and understand what I did or share the code with someone else. It offers the ability to reproduce my work.

I did end up using ARC GIS for visualizing my raster files, data processing such as clipping, and for geographically Weighted Regression.  R does offer a couple of packages for conducting geographically weighted regression, but I was having difficulty trouble shooting the packages in R.

What did you learn about statistics?

I learned ‘oh so many things’. I learned a number of spatial statistical methods throughout this course including semi-variograms, cross-variograms, patch metrics (FRAGSTATs), geographically weighted regression, and spatial overlay analysis. Semivariograms and cross-variograms are useful descriptive analysis to assess variance as a function of scale and generally describe landscape patterns. Patch metrics (FRAGSTATS) are also a useful tool for generating metrics and describing the distribution of disturbance severity across a landscape. I think there are many more patch metrics to explore in the SDM package that may be useful for examining some of the variability in disturbance severity across the landscape, overlapping disturbances, cumulative disturbance effects, and future successional pathways for disturbed forest landscapes. The geographically weighted regression was likely the most challenging analysis to make work and draw conclusions from. I will probably be thinking about this beyond this course. It was useful for seeing where the overlapping disturbances were positively and negatively correlated, but there was still a great deal of variability in the scatter plots. The spatial overlay analysis was useful to see were areas of beetle outbreak and/or no beetle outbreak overlap with areas burned and/or unburned.  The last piece tied in here was scatter plots of the geographically weighted regression base on the beetle outbreak, fire severity and spatial overlay, which showed a great deal of variability.

Response to Comments

I received a number of useful comments on my exercises and tutorials, some of which I have been able to address in this final blog post.

  1. Exercise 2: Julia’s Comment – regarding additional interpretation of the cross-variogram was addressed in the cross-variogram section above.
  2. Exercise 3: Julia’s Comment regarding SDMtools for future FRAGSTATs users – I tried to list some of the additional metrics that are provided in the output. The package was relatively easy to use, and I provided commented code with tutorial  2 that will hopefully be useful to a future user.
  3. Exercise 3: Julia’s Comment – “Useful to know about proportions of h,m.l severity. yes mean patch sizes differ by burn severity but not median patch size – what do the full distributions of patch sizes look like for the various severity classes? Note the high st dev for L and M – indicating that likely there are a couple of very large patches affecting the mean (not the median), hence not really a difference in patch sizes. Why not? What would control patch size?” – I have attempted to generate boxplots, but have not been able to make ones that are visually worth sharing. There seem to be lots of small patches and a few very large outliers. This is worth additional exploration, but I am out of time in terms of doing this for the class.
  4. Tutorial 1 comments: Peer Comments – Create many variograms; NDVI, EVI; look at different band combinations. These comments were all considered.  I had generated an NDVI layer early one, and after reading through some remoting sensing pieces regarding Landsat data and beetle outbreak I opted to use the pre-fire NBR raster.  We talk quite a bit about the fact that the semi-variogram takes a sample of the data.  I opted not to run multiple semivariograms, since the analysis is considered a descriptive tool.  The model run provided a general understanding of the spatial variance on the landscape that seem useful. I am still considering the use of additional bands, but chose to hold off on that for this project/class.
  5. Tutorial 2: Peer Comment – Identifying any connections between layers and considering other factors. I think that the next steps for the project will be to bring in topography and weather layers, and run additional spatial statistical analysis including topography and weather. Due to time this was outside of the purview of this project.




ESRI. 2011. Arc GIS Desktop: Release 10. Environmental Systems Research Institute, Redlands, CA.

Hart, S. J., T. Schoennagel, T. T. Veblen, and T. B. Chapman. 2015. Area burned in the western United States is unaffected by recent mountain pine beetle outbreaks. Proceedings of the National Academy of Sciences of the United States of America 112:4375–80.

Key, C. H., and N. C. Benson. 2006. Landscape assessment (LA): Sampling and analysis methods. USDA Forest Service General Technical Report RMS-GTR-164-CD:1–55.

McGarigal, K., A. Cushman, and E. Ene. 2012. FRAGSTATS v4: Spatial patterns analysis program for categorical and continuous maps. Univeristy of Massachusetts, Amherst, Massachusetts.

Naimi, B. 2015. Package “ usdm .” R Topics Document.

Reilly, M. J., C. J. Dunn, G. W. Meigs, T. A. Spies, R. E. Kennedy, J. D. Bailey, and K. Briggs. 2017. Contemporary patterns of fire extent and severity in forests of the Pacific Northwest, USA (1985-2010). Ecosphere 8:e01695.

Team, R. C. 2016. R: A language and environment for statisitical computing. R Foundation for Statistical Computing, Vienna, Austria.

VanDerWal, J. 2015. Species Distribution Modelling Tools: Tools for processing data associated with species distribution modeling exercises.


Print Friendly, PDF & Email

No Comments

No comments yet.

RSS feed for comments on this post.

Sorry, the comment form is closed at this time.

© 2019 GEOG 566   Powered by WordPress MU    Hosted by