Brian McMullen

CEE 6440 Fall 2005 Term Paper

Using GIS to depict TEUI survey results from the Shoshone National Forest, Wyoming

Introduction

My M.S. thesis research is coupled with full-time summer survey work with the United States Forest Service (USFS) in the Wind River Mountains of Wyoming.  The survey project follows the protocols of the USFS’s Terrestrial Ecosystem Unit Inventory (TEUI) and couples soils, plants, and geology data collected in the field in order to produce maps and interpretations for land management decisions.  Two summers (2004 and 2005) of field data have been collected, culminating in 253 ecological site descriptions. Each plot contains detailed information on the soil taxonomy, parent material, and plant community.  Remotely-sensed (Landsat 7) imagery combined with existing geographic information system (GIS) maps and digital elevation models were used for site selection.  Maps for this project will help further refine the data collected and provide ground-truthing for previously prepared maps that were based on geology and air photo interpretation.

Research Area

The Wind River Mountain Range is one of the numerous mountain chains that comprise the Rocky Mountains. They are part of a large anticlinal uplift characteristic of the Central Rocky Mountain Province. Measuring 190 km (120 mi) in length and 60 km (40 mi) in width, the “Winds” extend from South Pass in the south central part of Wyoming to the Gros Ventre Mountains to the north (Svalberg et al., 1997). This aptly-named Range receives Chinook winds from the Pacific that descend the east slope towards Lander and Dubois (Maheny, 1977). The Range trends northwest to southeast, with central latitude and longitude of 42º52´ N and 109º52´ W. Located in a sparsely populated area of Wyoming, the Winds straddle the Continental Divide and rises sharply out of the surrounding Green River and Wind River Basins, huge expanses of flatland that exist in stark contrast to the numerous 4000+ meter peaks and alpine glaciers of the Range.

Figure 1: Landsat 7 image of the Wind River Range and surrounding basins

Lander, Wyoming is considered the gateway town to the eastern slope of the Wind River Range. It is a 250 mile drive from Logan, Utah to Lander; both locales are shown in Figure 2.

Figure 2: Project location obtained from Google Earth imagery

The eastern slope of the Winds are predominantly managed by federal land agencies, including the USFS, the Bureau of Land Management (BLM), and the Bureau of Indian Affairs (BIA). Figure 3 shows the extent of the USFS’s Shoshone National Forest, with the research area and survey points demarcated in the southern portions.

Figure 3- USFS land in the eastern portion of the Wind River Range

Geologic substrate serves as soil parent material; for the survey area of the Shoshone National Forest, this connection was paramount to understanding the spatial distribution of soil and vegetation. Figure 4 shows a generalized cross section of the geology of the area.

Figure 4: West (Green River Basin) to East (Wind River Basin) geologic cross section of the Wind River Range (modified from Dahms et al., 2003)

Methods and Materials

Data Collected

The information collected at each survey plot is highly suitable for map production using ArcGIS.  Slope, aspect, and elevation were recorded at each plot. Soil profile characteristics described include soil texture, drainage class, moisture/temperature state, rock (surface and profile) content, and pH. Base saturation analyses were made for selected soil pits from the 2004 field season to aid in taxonomic classification of the soils and to better discern spatial nutrient patterns. Vegetation measurements made in each plot include percent species coverage, diameter at breast height (dbh) of representative tree species, tree cores for dendrochronology, and habitat type classification. Bedrock and surficial geology were described and compared to preexisting geology maps for the area.

Plot data has been entered into an Excel spreadsheet; this tabular data includes UTM coordinates and can therefore be projected geographically using a series of basic steps in ESRI’s ArcCatalog and ArcMap.

Figure 5: Data entered in Excel format

Data from the Excel spreadsheet was copied and saved as a database (.dbf) file. From this file format, a personal geodatabase was created using ArcCatalog, as shown in Figure 6.

Figure 6: Geodatabase, shapefile, and feature class creation in ArcCatalog

The digitized map shown in Figure 7 was provided to our survey crew prior to the start of the initial field season. I projected the sampling points from coordinates contained in the previously described geodatabase.

Figure 7: 2004 sampling points

Prior to the 2005 field season, areas with sampling deficiencies were identified; these are depicted in Figure 8.

Figure 8: Areas identified for sampling in the 2005 field season

The map in Figure 8 served as a guide for developing a sampling plan for the 2005 field season. Figure 9 shows how these sampling gaps were addressed.

Figure 9: 2005 Sampling points

All but 12 of the 253 sampling points are shown in Figure 10 below; the missing 12 are in a remote northern part of the forest and are not included for the majority of GIS analyses included in this project.

Figure 10: Combined 2004 and 2005 sampling points

Embedded within each point is a wide variety of information contained within the attribute tables of each plot. Using the identify tool in ArcMap, all of the plot data contained in the geodatabase can be accessed through the click of a button.

Figure 11: Plot data contained for each point feature

Survey Analysis using ArcGIS Spatial Analyst and Spatial Interpolation

Interpolation is the development of new data points from an existing set of measured points and is useful for delineating soil, vegetation, and geologic trends across the landscape. Interpolation assumes that close samples tend to be more alike than distant samples and uses mathematical algorithms to create continuous surfaces from point data (Wackernagel, 2003). The different algorithms used in ArcGIS Spatial Analyst provide various ways of estimating measurements at unmeasured locations (e.g. base saturation status of the soil where no measurement was made).

Geographic information systems provide an ideal framework to graphically depict spatially interpolated data. A summary of the various algorithms used in ESRI’s ArcGIS follows:

Inverse Distance Weighting (IDW)

Data from an unknown point is estimated from the sum of all neighboring locations. The simplest weighting function is inverse power: w(d)= 1/d^p
where p is a user-specified value greater than 0. Inverse distance weights the data points by the inverse of their distance from the estimation point, thereby giving greater influence to nearby data points than those farther away.

Mathematically, this is expressed as

“…where v₀ is the estimated concentration at (x₀, y₀, z₀), v_i is a neighboring data value at (x_i, y_i, z_i), d_i is the distance between (x₀,y₀,z₀) and (x_i,y_i,z_i), P is the power, and N(v₀) is the number of data points in the neighborhood of v₀.” (Spatial Analysis and Decision Assistance (SADA) Home Page).

Kriging

The term “Kriging” is named after the 2^nd greatest South African engineer the world has ever known, D.G. Krige. Krige developed the geostatistical technique to better predict the location of ore reserves and mineral concentrations at unobserved points. This technique was further refined by the French mathematician Georges Matheron in the 1960s (Wackernagel, 2003). One improvement over IDW is that predictions have explicit uncertainty calculations associated with them. Kriging uses least squares to find a best fit to a data set by minimizing the sum of the squares of the ordinate differences between the data and the fitted function.

Natural Neighbor Algorithm (NN)

Natural Neighbor interpolation is a weighted moving average technique utilizing geometric relationships to select and weight nearby points. Circumcircles are used to determine the number of natural neighbors; points are natural neighbors if they lie on the same natural neighbor circle (SADA, 2005). Weights (wi) are interpolated using Delaunay triangulation and depend on the area about each of the data points rather than the distance between data points. The equation for the NN interpolation is:

where: G(x,y) is the NN estimation at (x,y), n is the number of nearest neighbors used for interpolation; f(xi,yi) is the observed value at (xi,yi); and wi is the weight associated with f(xi,yi). (SADA, 2005).

Base Saturation Interpolation

Base saturation is the percentage of soil cation exchange sites occupied by the basic cations calcium (Ca²⁺), magnesium (Mg²⁺), potassium (K⁺), and sodium (Na⁺). The relative ratio of base cations to acidic cations such as hydrogen (H⁺) and aluminum (Al³⁺) is important because plant toxicities and nutrient deficiencies occur when imbalances occur. Soil parent material, atmospheric acid deposition, and management practices all can impact the base saturation status of soils. Base saturation was determined for 56 of the 117 soils from the 2004 field season using ammonium acetate extraction at Utah State University’s Soil Genesis Laboratory. Data was used to more confidently classify the soil according to the guidelines of the Keys to Soil Taxonomy (USDA, 2003). Figure 12 shows this sampling subset.

Figure 12: 2004 sample points with base saturation data

A shapefile was created for this subset, and a mask was used to confine subsequent interpolations to the southern portion of the forest (REVISEDFINAL.shp).

Figure 13: Step 1 of mask procedure

Figure 14: Step 2 of mask procedure

The three previously described interpolation algorithms were used to predict base saturation status for the study area. All three methods captured the general trend of decreasing base saturation going east to west across the survey area. This corresponds with the geologic distribution of calcareous limestone and dolomite in the flanking rocks of the eastern portion of the range, transitioning to acidic granite and gneiss in the central Precambrian core. IDW was able to capture the smallest scale variability, as shown in Figure 15.

Figure 15: Base saturation predicted by IDW interpolation

The NN algorithm produced a similar distribution pattern, but with one area in the northeastern area of the forest with an broader distribution of higher base saturation soils than would be expected given the acidic nature of the granitic parent material. This area is highlighted in red in Figure 16.

Figure 16: Base saturation predicted by Natural Neighbor interpolation

Kriging provided the coarsest prediction of base saturation distribution due to the relative lack of data points for interpolation.

Figure 17: Base saturation predicted by Ordinary Kriging

For taxonomic purposes, 60% base saturation is the dividing line between Dystro- and Eutro- Inceptisols (e.g. Typic Dystrocryepts vs. Typic Eutrocryepts). Figure 18 shows the geographic divide that represents the partitioning of the forest into areas of high ( ≥ 60% base saturation, Eutro-) and low (< 60% base saturation, Dystro-) base saturation soils.

Figure 18: IDW Base Saturation interpolation using Dystro- and Eutro- breaks

Soil pH distribution using Spatial Interpolation

Soil pH was measured for all 253 soils using colorimetric indicator dyes; an overview of the trends observed shows a similar spatial distribution as observed with base saturation. Calcareous parent material weathers to produce more Ca and Mg in the soils of the eastern flank rocks (Paleozoic limestones and dolostone shown in Figure 4). Traveling east to west in the study area leads to an overall decrease in pH due to the lower concentration of base cations on the soil exchange sites relative to acidic cations. A notable exception to this general trend lies in the most southerly location of the study area (near Louis Lake); here, granodiorite and metasedimentary bedrock is intruded by a complex of mafic and felsic dikes that produces soils intermediate in pH and base saturation values between the flanking calcareous soils to the northeast and the acid granite soils of the high elevation core rocks to the west/northwest. Figures 19 and 20 depicts this trend and also demonstrate the effect of having more sampling points for the algorithms used in interpolation.

Figure 19: pH interpolated by IDW using only 2004 data points (n= 117)

Figure 20: pH interpolated by IDW using combined data from 2004 and 2005 (n= 241)

A lack of sampling points in the northwestern part of the study area (shown in red) led to lower-than-expected pH values for this area; for most of the forest, the NN algorithm reproduced the observed trend of decreasing pH concentration going east to west, with intermediate values in the southern portion of the forest.

Figure 21: Natural Neighbor pH prediction (2004/2005 combined)

Broader underestimations of pH occurred using NN interpolation with pH values from 2004 only.

Figure 22: Natural Neighbor pH prediction (2004)

Predicting clay distribution using Spatial Interpolation

Soil texture determinations were made in the field for all of 2004 and ~ 85% of 2005 soils. Percentages of sand and clay were determined by hand, leading to identification on the USDA’s Soil Textural Triangle. Clay imparts many important chemical and physical properties to the soil. The trends of higher pH and base saturation in the eastern flanking rocks of the forest correspond with generally higher levels of clay in this area as well. This makes sense as clay has greater surface area (and hence higher amounts of negatively charged exchange sites) to retain base cations, which in term raises the pH (lowers the amount of H⁺ in the soil). The calcareous parent materials has more impurities that weather to produce clay than the igneous core rocks. Figure 23 shows why I chose to leave out the 12 northern points for interpolation purposes. This figure also illustrates gaps in the 2005 clay percentages (the missing textures for 15% of the samples); some points in areas with red circles had a value of 0% for clay since texture determinations have not yet been made for these points.

Figure 23: Clay % interpolation for 2004/2005 (n= 253) points

Better clay predictions were made when 2005 clay percentages were left out altogether, as shown in Figures 24 and 25.

Figure 24: Clay distribution using NN interpolation of 2004 data points

Figure 25: Clay distribution using IDW interpolation of 2004 data points

Miscellaneous survey results depicted using Spatial Interpolation

To conclude this project, I experimented with interpolating fields other than base saturation, clay, and pH. Figure 26 shows a soil depth prediction map based on how deep each hand-dug soil pit was. Our goal was to excavate a pit to at least 1 meter where possible; this was usually achieved unless bedrock or densic material (e.g. dense basal till) was encountered. In general, the eastern flanks rocks weathered to produce deeper soils, with a few exceptions.

Figure 26: Soil depth interpolation using IDW

Elevation was recorded at each site using a handheld global positioning system device. Using all of the 2004 and 2005 elevation points from the survey for the southern portion of the forest, I was able to construct a digital elevation model using IDW interpolation.

Figure 27: IDW interpolated Digital Elevation Model

Slope percentage was measured at each site using an inclineometer. Using the combined 2004 and 2005 data set for the southern portion of the forest, slope maps were derived using NN and IDW interpolation (Figures 28 and 29, respectively).

Figure 28: Slope prediction using NN interpolation of 2004/2005 data points

Figure 29: Slope prediction using IDW interpolation of 2004/2005 data points

One measured parameter that will need much greater attention is the geographic distribution of soil temperature. A confounding variable in this is disparate weather patterns in the two field seasons. Overall, the weather was much cooler and wetter in the 2004 field season than the 2005 field season. Soil temperature was recorded using a thermometer inserted into the soil profile at a depth of 50 cm (when possible). The lack of identifiable trends with regards to elevation, parent material, and soil type is shown in Figure 30.

Figure 30: Soil temperature interpolation by IDW (2004/2005 southern points)

Conclusions

The geographic trends analyzed with spatial interpolation are summarized in Figure 31 below. Completion of all of the 2005 data points and concentration on areas with exceptions to the general spatial trends observed will further aid in delineating a complex landscape into more simple management units.

Figure 30: General geospatial trends for the southern unit of the Shoshone National Forest

Works Cited

Dahms, D.E., Hall, R.D., Shroba, R.R., Sorenson, C.J., Lynch, E.A., and Applegarth, M.T. 2003. The Rocky

Mountain glacial model: The Wind River Range, Wyoming: in Easterbrook, D.J., ed., 2003. Quaternary

Geology of the United States: INQUA 2003 Field Guide Volume, Desert Research Institute. Reno, NV, p. 345-364.

Mahaney, W.C. 1978. Late-Quaternary Stratigraphy and Soils in the Wind River Mountains, Western

Wyoming. in Quaternary Soils, York University Press, Toronto, Canada, p 225-264.

Spatial Analysis and Decision Assistance (SADA) Home Page http://www.tiem.utk.edu/~sada/help/TH_33.htm

Svalberg, T., Tart, D., Fallon, D., Ferwerda, M., Lindquist, E., and Fisk, H. 1997. Bridger-East Ecological Unit

Inventory, Volume 1, USDA NRCS and University of Wyoming.

United States Department of Agriculture. 2003. Keys to Soil Taxonomy, Ninth Ed. Soil Survey Staff of the Natural Resource Conservation Service.

Wackernagel, Hans. 2003. Multivariate Geostatistics: An Introduction with Applications. Springer, Berlin.