Having endured much frustration trying to use Geographically Weighted Regression to determine relationships between the patch metrics I calculated using FRAGSTATS, I’ve decided to change things up a bit for my final blog post. I’ll continue working with Geographically Weighted Regression, but rather than applying it to FRAGSTATS-derived metrics (shape index, edge:area ratio etc.), I will experiment with a different LandTrendr output: forest biomass.

Dataset:

Explaining more about the dataset here will help me articulate the research question below. So, the biomass layer is a raster calculated using a Tassled Cap Transformation, which is a method of enhancing the spectral information content of Landsat imagery. It is essentially a measure of per-pixel brightness (soil), greenness (vegetation) and wetness (soil and canopy moisture). I will be using yearly time-series biomass layers and “time-stamped” clear-cut disturbance patches.

Research Question:

I’m still not sure how well I can articulate the question, but here goes: Is there a statistically significant relationship between the mean biomass value within a clear-cut patch at the timing of the clear-cut, and the mean biomass within that same patch, before and after the clear-cut?

 

Hypothesis:

Clear-cutting obviously results in a loss of biomass, and I expect that quantities of biomass within a clear-cut patch before, during and after a clear-cut will exhibit a significant relationship.

Approach:

While I had easy access to the biomass data, creating a dataset of disturbance patches with attributes for biomass before, during, and after the timing of each clear cut was a carpal-tunnel inducing task. I ought to have approached it programmatically, and I did try, but my Python skills were lacking. I ended up using a definition query to filter disturbance patches by year, and then ran three Zonal Statistic operations on the biomass layers (one for the year before a given set of clearcuts, one during, and one after). I then joined each biomass calculation back to the clear-cut patches. Below is an attribute table for one set of disturbance patches (note the three mean biomass values) and an example of a disturbance patch overlaid on a “before, during, and after” set of images. I did this for three sets of yearly clear-cuts, and then merged them into one dataset of roughly 700 features.

attrib befaftdur

I then ran Geographically Weighted Regression, with “mean biomass after clear-cut” as the dependent variable, and “mean biomass before” and “mean biomass at timing of clear-cut” as explanatory variables.

Results:

I experimented with multiple combinations of the three biomass mean variables, and also tried adjusting the kernel type. The most significant run was that on the left, which had parameters as described above.

r2         r2

gwr_r2gwr2gwr5

gwr3

 

gwr4

Significance:

While it was satisfying to finally produce some significant statistics, I recognize that the analysis is not groundbreaking. While change in biomass is certainly of interest to forest managers and ecologists, the way in which it was calculated here (as a mean of an annual snapshot within a patch) may not have significant implications.

What I learned:

If this course is nicknamed “Arc-aholics Anonymous” then you could say I had somewhat of a relapse, as most of my analyses throughout the quarter made use of tools I had used in ArcMap before. That said, I gained a much more thorough understanding of their functionality, and feel I have a better command over interpreting their sometimes perplexing results. I now have a much better idea of the types of datasets and variables that lend themselves to certain methods, and the experience of working with a large dataset of Landsat time-series-derived forest disturbance will be invaluable to my research moving forward. I learned a great deal from others in the course and am glad to have made some new contacts (especially you coding gurus). Some of the work produced in this course was truly outstanding and I feel inspired to hone my own skills further, particularly with open-source software.

For this final tutorial post, I’ll be describing a workflow in which FRAGSTATS is used to generate a number of continuous metrics of disturbance patches. Then, using these metrics, a Grouping Analysis is performed in ArcMap to identify groups of patches that have multiple similarities.

FRAGSTATS is a spatial pattern analysis program developed by landscape ecologists for quantifying the structure (i.e., composition and configuration) of landscapes at multiple scales (see this link to the documentation). Through the lens of a landscape ecologist, a landscape may be considered as an assemblage of patches whose spatial homo/heterogeneity characterize the landscape they comprise. While the patches of interest to landscape ecologists are often habitat types, they may represent any spatial phenomenon, including forest disturbance.

The installation of FRAGSTATS is very straightforward (link to downloads page), and the GUI is friendly! Below I outline the steps for adding data to FRAGSTATS and generating patch metrics:

1) Add a layer. FRAGSTATS accepts raster images in a variety of formats (again, see documentation here). I worked with GeoTIFF (.tif) files representing disturbance patches.frags3frag_input

2) Select metrics to calculate. Descriptions of each metric can be found in the documentation.

results

3) Generate Results. Simply click the run button on the main toolbar, and view the results.

results2

If your goal is to attach these tables back to your input data for mapping/analysis, in a GIS for example, then it is crucial to generate a “patch ID file”. To do this, simply check the box for “Generate patch ID file” under the Analysis parameters tab:

Capture

The patch ID file and associated tables will be saved to the directory you choose. Note that here I’ve checked only the box for Patch metrics. The patch ID file will have a suffix of “_id8” appended to whatever name your input file is, and it’s associated extension (“input”_id8.tif). The patch metrics file will have a .patch extension. Open the .patch file in Excel or the spreadsheet editor of your choice, delimit by comma, and save it as a file type that ArcMap will recognize, such as .csv or .txt. I suggest removing the “LID” field which contains the file path where your initial input raster resides.

4) Join output .patch file to patch ID file. In ArcMap, bring in both the patch ID and copy of the patch file in .csv or .txt format. Then, proceed with a tabular join:

joinjoin

Right-click the _id8 file, click Join…in the prompt window that appears choose “Value” for the field in the patch ID file, and PID for the field in the table of patch attributes. Click OK. The patch attributes are now joined to each patch in the patch ID file. Open the attribute table of the patch ID file to verify. Remember to export this file to make the join permanent.

*Note: if you don’t see “Value” as an option for the join field, you may have to run the “Build Raster Attribute Table” tool on your patch ID file.

5) Proceed with Mapping/Analysis!

*Note that the FRAGSTATS tabular outputs are useful in their own right. Below I’ve charted “class” level metrics (as opposed to patch level) that show trends in clear-cut patches over time. In this case, I filtered the disturbance data for clear-cuts only, and had FRAGSTATS generate outputs by year (year is the “class”).

charts

 

 

 

 

Moving forward with the disturbance centroids, I ran two other analysis tools in ArcMap that identify not only spatial clustering, but clustering of similar numeric (not necessarily continuous) variables attributed to the centroids. These tools are Cluster and Outlier Analysis (Anselin Local Moran’s I) and Hot Spot Analysis (Getis-Ord Gi*). While these tools essentially answer the same question (where are clusters of similar values?) they calculate this in slightly different ways. See this link for more information. For both tools I used year of disturbance as the input value to be evaluated:

yr_moransihotspot

The map on the left shows the results of the Local Moran’s I tool. High-High and Low-Low Clusters represent statistically significant clusters of high and low values, whereas High-Low and Low-High clusters are outliers representing either high values surrounded by low values, or vice versa. The results of the Hot Spot Analysis show similar patterns of significant clustering by year. These patterns are indicative of geographic shifts in timber harvest concentration, which could reflect management decisions. Using “Select by Location”, the land use designations associated with significant Moran’s I clusters were tabulated. Low-Low Moran’s I clusters are on the left (1985-1989) and High-High clusters (1996-2012) on the right:

moransiLLmoransiHH

Here we see an expected increase in clustering of clear-cuts and partial harvests on Non-Forest Service (private) and Matrix lands, yet it is important to keep in mind that these land-use designations did not exist until 1994.

Following this analysis of clustering by year, I was interested in attaching continuous variables to the actual disturbance patches, rather than their centroids. Toward this end, I brought in an another output from the LandTrendr algorithm: disturbance magnitude, which is a function of the slope of change in spectral values for each pixel over a yearly time step. Since magnitude is measured on a per-pixel basis, I used the Zonal Statistics tool in ArcMap to calculate the mean disturbance magnitude within each patch, and attach it to the patch attributes. I then ran Hot Spot Analysis again, with mean disturbance magnitude as the input value.

magn_hotspotsmagn

The distribution of high and low magnitude disturbances is interesting, especially when overlaid on the land use designations. As shown in the map on the right, a cluster of high magnitude disturbance is associated with a section of adjacent private lands (in grey) and a cluster of low magnitude disturbance is associated with Matrix land (in orange). This may indicate a higher proportion of clear-cutting (high magnitude disturbance) on private lands, and more partial harvesting (lower magnitude disturbance) on Matrix lands.

To begin my analysis of forest disturbance patterns, I wanted to get a broad scale understanding of how clear-cuts and partial harvests are distributed throughout Willamette National Forest, as well as the land use designations and forest districts they spatially coincide with. Using ArcMap I ran the Kernel Density tool on cumulative disturbances from 1985 to 2012. This first required converting my disturbance patch data (rasters) into polygons (vectors). Additionally, Kernel Density requires point or polyline inputs, so I also generated disturbance centroids from the polygons. Below I outline the function, results and my interpretations:

Kernel Density calculates a magnitude per unit area using a kernel function to fit a smoothly tapered surface to each input feature. Here are the results for two runs of this tool:

disturb_density3disturb_density2

Below are the parameters used for each kernel density output above. In the first run on the left, the output cell size of 250 meters results in a surface that is less smooth than that of the second run on the right, which has a much smaller output cell size. Additionally, in the second run, the optional parameter for a search radius was set to 1600 meters, a value that is roughly double that of the average nearest-neighbor distance between disturbance centroids. This creates a density surface that is more appropriate for the scale of the map (on the right). Both maps give an indication that higher densities of clear-cuts and partial harvests between 1984 and 2012 occurred on Matrix Lands and Non-Forest Service Lands, which are mostly composed of private and industrial landowners. These results are not surprising, given that timber harvest is the primary function of Matrix Lands, and that private and industrial landowners are inclined to produce timber.

CaptureCapture3

*It is important to note that after converting the disturbance patches to polygons, I used the Normal QQ Plot and Histogram functions (built in to the Geostatistical Analyst extension for ArcMap) to remove outliers using the “Shape_Area” attribute field. Polygons of extremely low area (less than eleven 30-meter Landsat pixels) were removed because they likely represent incorrectly classified pixels. Polygons of extremely high area were also removed, because they would not be accurately represented by a single centroid point.

Mapping cumulative clear-cuts and partial harvests over the 27 year period between 1985 and 2012 paints an interesting picture, but in order to better assess the effects of forest governance on landscape patterns, it is probably more interesting to map disturbances temporally. Using the Definition Query function in ArcMap, I filtered the disturbance data at 5-year intervals, and ran the Kernel Density tool again for each interval:

dens_time

Bringing in the temporal dimension reveals some interesting changes in the density of forest disturbance that may correspond with significant events in the history of forest governance. For example, 1990 shows very high density of clear-cuts and partial harvests, which indicates a peak in timber harvest activity prior to the May 29, 1991 Dwyer Injunction, which banned timber sales. Following the 1994 Record of Decision for the Northwest Forest Plan, the maps for 2000, 2005 and 2010 show the expected drop in overall disturbance density, but interesting peaks associated with certain forest districts. For one final run of the Kernel Density tool, I mapped cumulative disturbance during the period of the Dwyer Injunction (1990-1994). The result shows higher disturbance density than expected, but I assume that there is lag between the timing of timber sales and when sold forest land is actually harvested. Thus, this map likely shows clear-cuts and harvests on forest lands sold prior to the injunction, but may also be indicative of which forest districts are more inclined toward timber production.

91_94_dens

1. Background and Research Question:

The goal of my work in this course is to assess the influence of forest governance on spatial patterns of forest disturbance. Forest governance can be understood as forest-related decisions, their implementation and resulting effects within a given institutional setting, whereas forest disturbances are events that cause change in the structure and composition of a forest ecosystem. For this course, I’ll be focusing my analysis on disturbances associated with timber production, e.g., clear-cutting and partial harvests. My study area is Willamette National Forest, which encompasses roughly 6,800 square km in the central portion of Oregon’s West Cascades.

Governance of the Willamette and other federally managed forests of the Pacific Northwest is shaped largely by the Northwest Forest Plan (NWFP). A key aspect of the NWFP is a system of land use designation (LUD) in which spatially explicit zones are managed according to a single or dominant management priority (Charnley, 2006). For example, wilderness areas are designated as “Congressionally Withdrawn” and are thus protected from timber harvest, while “Matrix Lands” are those on which timber harvest is concentrated. The various LUDs, along with interspersed state, private and other lands, form a mosaic of ownership and management priorities that are manifested as disturbance patterns in the forest landscape. And so, the research question I’ll be addressing is: How do land use designations and ownership influence patterns of disturbance in Willamette National Forest?

2. Datasets:

I’ll be relying primarily on a Landsat imagery time-series (30 meter spatial resolution, 1985-2012) which has been processed using a change-detection algorithm called LandTrendr (Kennedy et al., 2010). Outputs from LandTrendr include disturbance patches categorized by agent of change (e.g., clear-cut, partial harvest, etc.) as well as their timing, duration and magnitude. This data will be clipped to the boundary of Willamette National Forest. To structure my analysis, I’ll be using vector data representing the administrative forest boundaries, and the LUDs and ownerships within them.

wnf_ludsdisturbs

3. Hypotheses:

My general hypothesis is that patterns of disturbance will vary according to land use designation across space and time. For example, on Matrix Lands, I expect to see spatial clustering of clear-cuts that will increase during the period of the NWFP’s implementation (from 1994 onward), while on Adaptive Management Areas, I expect to see significantly fewer, more dispersed clearcuts, but increases in partial harvests. 

4. Approach:

I will analyze spatial patterns of clear-cut and partial harvest disturbance patches within each LUD (e.g., clustering), as well as the spatial characteristics of these disturbance patches (e.g., edge-to-area ratio). My analysis will be done primarily in ArcGIS, but may also include Python scripting, statistical analysis in R, or “patch analysis” using FRAGSTATS.

5. Expected outcome:

I will produce maps and graphs of disturbance patterns by LUD over time. Hopefully this will help visualize whether or not the NWFP has been implemented as intended. 

 6. Significance:

In the Pacific Northwest, forest governance is shaped largely by the federally mandated Northwest Forest Plan (NWFP), which initiated a momentous shift in forest management priorities; from the provision of sustained timber harvest to protection of ecosystems (Thomas, 2006). The NWFP is implemented through an adaptive management strategy that must identify high priority inventory and monitoring objectives needed to assess the plan’s effectiveness over time (FEMAT, 1993). Ideally, this project and my overall research will contribute to this ongoing assessment.

7. Experience levels with…

ArcGIS = high

Modelbuilder, Python = medium

R = low

References:

Charnley, S., & Pacific Northwest Research Station. (2006). Northwest Forest Plan, the first 10 years (1994-2003) : Socioeconomic monitoring results (General technical report PNW ; 649). Portland, OR: U.S. Dept. of Agriculture, Forest Service, Pacific Northwest Research Station.

Kennedy, R. E., Z. Yang, and W. B. Cohen. 2010. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr – Temporal segmentation algorithms. Remote Sensing of Environment 114:2897-2910.

Thomas, J., Franklin, J., Gordon, J., & Johnson, K. (2006). The Northwest Forest Plan: Origins, Components, Implementation Experience, and Suggestions for Change. Conservation Biology, 20(2), 277-287.

Report of the Forest Ecosystem Management Assessment Team (FEMAT, 1993)