Using an Energy Gradient Metric to Predict Habitat Selection

Eric McLeskey

GIS in Water Resources Term Project

12-05-03

Introduction

Drift feeding stream salmonids have been shown to use velocity gradients to conserve energy and maximize energy intake. By maintaining positions in slow moving water, a fish can minimize it’s energy expenditures while taking advantage of the higher drift-food offerings of faster moving water.

Probably the most widely used fish habitat model today is the Physical Habitat Simulation System (PHABSIM) component of the USFWS Instream Flow Incremental Methodology. PHABSIM is designed to predict the quantity and quality of habitat for single species and life stage of fish as it relates to flow. PHABSIM, in the simplest sense, uses depth, velocity and substrate suitability criteria as model input before computing the useable area for a particular flow. Because the velocity suitability criteria refer to velocity magnitudes only, there is currently no way within PHABSIM to assess the importance of different velocity gradients around each point. A habitat metric that can be used to measure the variability in the flow field around each point would be a useful addition to PHABSIM.

David Crowder and Panayiotis Diplas developed an energy gradient metric and showed that it could differentiate between points in a flow field with the same magnitude but with different gradients of adjacent flow (Crowder and Diplas 2000). Their metric:

Is a description of the spatial change in kinetic energy between two points in a flow field an arbitrary distance of Ds apart. Their purpose in developing the metric did not include determining a biological meaningful Ds for stream fishes and did not include a recommendation of which direction (lateral, longitudinal, vertical) the metric should be applied. They simply showed that it has utility in describing energy gradients for points in flow fields.

Objectives

My research had two main objectives:

1. To apply the energy gradient metric to a velocity data set.

2. To conduct a sensitivity analysis for different values of Ds by comparing fish observations with metric outputs.

To accomplish the first objective within GIS, it was necessary to run multiple grid calculations, each involving two values from the same grid. The second objective was met by correlating positive metric values with fish observation data.

Materials and Methods

Study Site and Data

The GIS data used for this project were original assembled by Greg Guensch as part of his Masters Thesis (Guensch 1999). The study site was a 250 m section of the Blacksmith Fork River located in northeastern Utah (Figure 1). The river at this location is third order and is at an elevation of approximately 1750 meters.

Figure 1. Study site.

The GIS data sets utilized for this project were a velocity grid of the study reach with a resolution of 0.3048 m (1ft) and a fish observation shapefile for 254 mm brown trout. Fish observation data were collected in the summer months by visual observations both in the stream (e.g., placing fish tags while snorkeling) and out of the stream (e.g., from an observation tower).

Project Development

Metric Application

The energy metric proposed by Crowder and Diplas can be applied in any direction. For the purpose of this study, I assumed that lateral velocity gradients were most important and considered <s as the direction that a fish would move perpendicular to flow (Figure 2). Another assumption I made was that all the flow directions were approximately parallel with the bank edges. While this is probably not an accurate representation of the flow dynamics of the reach, in absence of any flow direction data, this was a simple assumption and facilitated the metric computations.

Figure 2. Transformation required for proper application of the energy gradient metric.

There is a deviation from what is observed in a natural system in regard to fish feeding and how a grid calculation that loops through rows and columns in an ESRI grid would be performed. This problem was that the distance <s would be made in the rectangular x direction, not necessarily in the lateral direction of the stream especially if it meanders (Figure 2). To handle this problem, a mesh was created that allowed each point in the reach to be transformed into x’ and y’ coordinates (Figure 2).

Project Outline

To meet the objectives of this project it was necessary to do quite a bit of data manipulation. Although each step in the data processing was not terribly difficult, there were many file formats and folders that needed to be organized so that the processing could move efficiently. An overall outline of the processing steps from start to finish for this project shows that many software tools, some outside of GIS, were required to complete this project (Figure 3).

Figure 3. Flow chart showing data processing steps.

Project Implementation

Creating The (x’,y’) Mesh

A point mesh was created using Map Window 3.0 using a utility called Stream Builder. This utility allows a user to import a point shapefile of a river reach, draw a series of individual lines to represent the center line, smooth the joins of the individual lines, and builds a mesh perpendicular to the centerline (Figure 4). To use this utility it was first necessary to use ArcToolbox to create a velocity point shapefile from the velocity grid.

Figure 4. Building the mesh using Map Window’s Reach Builder.

An exploded view of the mesh created in Map Window shows that the points in the x’ direction are uniformly spaced (Figure 5). Obviously this cannot be true for the y’ coordinates, however, because only lateral velocity gradients were being considered important, the spacing of the points in the y’ direction was not important for the metric calculation.

Figure 5. Exploded view of the mesh created in Map Window.

Velocity Interpolation Using Terra Model

Terra Model was used to interpolate velocities from the velocity point shapefile created with ArcToolbox onto the mesh point shapefile within the boundary of the study reach (Figure 6). A boundary shapefile was manually created in ArcView 3.2 using the velocity point shapefile. Terra Model used a TIN interpolation method to output an interpolated velocity point shapefile and a table that included the x’, y’, and interpolated velocity values in it.

Figure 6. Terra Model inputs.

The purpose of using Terra Model was to assign velocity values to the mesh created in Map Window and then to take those values and make an ESRI grid using the Ascii to grid utility in ArcToolbox. This transformation made the y’ values uniform and had the overall effect of straightening the reach (Figure 7). This grid was then imported into ArcMap so that the metric calculations could be performed using a Visual Basic macro.

Figure 7. Straightened Reach.

Grid Calculations Using a Visual Basic Macro in ArcMap

To perform a calculation on a grid using multiple values from the same grid, access the individual grid cells was required. A Visual Basic Macro from ESRI’s ArcObject developer help was found that could accomplish this task with some code modification (for full source code see Appendix: Energy Metric Code).

Because the code I used in my project can be modified for many purposes, I will take a moment and explain how to use the code available in the appendix. Macros are very simple to run in ArcMap. Simply select Macros from the Tools menu in ArcMap, type in the name of your macro (Energy Metric in this case) and click Create (Figure 8).

Figure 8. Creating a Macro in ArcMap.

This will bring up a Visual Basic form in which the contents of the Appendix can be pasted (Figure 9). The code in the appendix includes the beginning and ending Subroutine commands so be sure to delete “Sub EnergyMetric()…End Sub” that was placed their by default in the form window before running the macro. The code in the appendix is documented in the places where the user is required to input information unique to them such as file names and file paths. Documentation within the code shows where a user can modify the equation performed on each cell so that the utility of the code is extendable beyond the energy metric used in this project.

Figure 9. Visual Basic Form window in ArcMap.

Once the code has been modified and the file paths and file names have been specified, select Run from the Visual Basic Project window (Figure 10).

Figure 10. Run Macro.

Results

The metric was run six times using <s values ranging from 1-6, which produced a grid for each value of <s. Once the metric calculations were completed, Ascii tables were created using the grid to Ascii utility in ArcToolbox for each metric grid. The resulting Ascii tables were manipulated with the PFE text editor individually and the x’ and y’ values from the mesh created in Map Window were added to each table. Each table was then imported into ArcView 3.2 and opened as an event theme after which the fish observation shapefile was added (Figure 11).

Figure 11. Metric results for <s=2 and fish observations for 254 mm brown trout.

Visual observations were made to assess a correlation between positive metric values and the observed fish locations. By counting the number of fish in positive locations for each metric calculation, a moderate correlation was found for <s values of 2 and 3 (Table 1.)

Table 1. Percent of brown trout observed in locations with positive metric value.

<s	1	2	3	4	5	6
	51	63	63	54	48	41

None of the values of <s seemed to give a high correlation between observed fish locations and positive metric values. Some reasons this may have occurred could be due to high productivity in the river. If fish are well fed, they may not be concerned about conserving energy. Also, the flow field itself was fairly uniform having many areas of low velocity. It may be that the trout were just not expending that much energy in the areas with negative metric values.

Conclusions and Potential Further Research

The scope of this project could be expanded to include multiple species and size classes or include multiple study sites. The energy gradient metric did not seem to be able to predict where fish were located in the Blacksmith Fork River in the context in which it was applied. However, further reach is required before any sound conclusions can be made about it’s utility in predicting fish locations in streams.

Acknowledgements

Much thanks and appreciation goes to Greg Guensch for providing me with the data so that this project was possible. Also, thanks to Craig Addley, Jennifer Ludlow and Mark Winkelaar at the Institute for Natural Systems Engineering (INSE), Utah State University for helping me in the development stages of this project to conceptualize how the project could be handled and for offering technical expertise and advice when I ran into problems. Last but not least thanks to Shannon Clemens also at INSE for her help interpolating the velocities for this project with Terra Model.

References

Crowder, D.W. and P. Diplas. 2000. Evaluating spatially explicit metrics of stream energy gradients using hydrodynamic model simulations. Canadian Journal of Fisheries and Aquatic Science 57: 1497-1507.

Guensch, G. 1999. Validation of an individual-based, mechanistic, habitat selection model for drift-feeding salmonids. M.S. Thesis, Utah State University, Logan, Utah.

Appendix: Energy Metric Code

The following code was modified to meet the purposes of this term project from a macro found from the ArcGis desktop help following the path: Help for developers/Arc Objects developer help/Object Model Overviews/Spatial Analyst Extension/Using the PixelBlock for Custom Analysis.

Public Sub EnergyMetric()

' This example calculates an energy gradient metric value

' for each cell in a grid. The metric equation is:((V1+V2)/2)*(V2-V1)/delta S)

' where delta S is the distance between grid cells.

' It demonstrates how one would use the PixelBlock

' object to perform a grid calculation.

' Declaring variables

' sInputPath: path of the input raster dataset

' sInputFileName: name of input raster dataset

' sOutputPath: path of the output raster dataset

' sOutputFileName: name of the output raster dataset

Dim sInputPath As String

Dim sInputFileName As String

Dim sOutputPath As String

Dim sOutputFileName As String

' Initializing variables

' This is where you insert the file names and extensions of the input

' and output grids for your project.

' Note that for this macro to work sInputPath and sOutputPath

' need to have the same extension.

' The sOutputFileName can be whatever you want it to be but you must

' rename it if you are going to do multiple calculations. You get a

' runtime error if you try to overwrite the file.

sInputPath = "C:\Metrics"

sInputFileName = "vel_grid_new"

sOutputPath = "C:\Metrics"

sOutputFileName = "MetricCalc2"

' Open input raster dataset

Dim pRWS As IRasterWorkspace2

Dim pWSF As IWorkspaceFactory

Set pWSF = New RasterWorkspaceFactory

If Not pWSF.IsWorkspace(sInputPath) Then Exit Sub

Set pRWS = pWSF.OpenFromFile(sInputPath, 0)

Dim pInputRasterDS As IRasterDataset

Set pInputRasterDS = pRWS.OpenRasterDataset(sInputFileName)

' Get RasterBand from the input raster and QI raster

' properties interface from input raster dataset

Dim pInputBand As IRasterBand

Dim pInputBandCol As IRasterBandCollection

Set pInputBandCol = pInputRasterDS

Set pInputBand = pInputBandCol.Item(0)

Dim pInputRasProps As IRasterProps

Set pInputRasProps = pInputBand

' Create a default raster from input raster dataset

Dim pInputRaster As IRaster

Set pInputRaster = pInputRasterDS.CreateDefaultRaster

' Create output raster dataset

If Not pWSF.IsWorkspace(sOutputPath) Then Exit Sub

Set pRWS = pWSF.OpenFromFile(sOutputPath, 0)

Dim pOrigin As IPoint

Set pOrigin = New Point

pOrigin.X = pInputRasProps.Extent.XMin

pOrigin.Y = pInputRasProps.Extent.YMin

Dim pOutputRasterDS As IRasterDataset

Dim outputPixelType As rstPixelType

Select Case pInputRasProps.PixelType

Case PT_CHAR

outputPixelType = PT_LONG

Case PT_COMPLEX

MsgBox "Operation not supported on complex data", _

vbOKOnly, "Neighborhood Notation"

Exit Sub

Case PT_DCOMPLEX

MsgBox "Operation not supported on complex data", _

vbOKOnly, "Neighborhood Notation"

Exit Sub

Case PT_DOUBLE

outputPixelType = PT_FLOAT

Case PT_FLOAT

outputPixelType = PT_FLOAT

Case PT_LONG

outputPixelType = PT_LONG

Case PT_SHORT

outputPixelType = PT_LONG

Case PT_U1

outputPixelType = PT_LONG

Case PT_U2

outputPixelType = PT_LONG

Case PT_U4

outputPixelType = PT_LONG

Case PT_UCHAR

outputPixelType = PT_LONG

Case PT_ULONG

outputPixelType = PT_LONG

Case PT_USHORT

outputPixelType = PT_LONG

End Select

Set pOutputRasterDS = pRWS.CreateRasterDataset(sOutputFileName, "GRID", pOrigin, _

pInputRasProps.Width, pInputRasProps.Height, _

pInputRasProps.MeanCellSize.X, _

pInputRasProps.MeanCellSize.Y, 1, outputPixelType, _

pInputRasProps.SpatialReference, True)

' Create a default raster from output raster dataset

Dim pOutputRaster As IRaster

Set pOutputRaster = pOutputRasterDS.CreateDefaultRaster

' Get RasterBand from the output raster

Dim pOutputBand As IRasterBand

Dim pOutputBandCol As IRasterBandCollection

Set pOutputBandCol = pOutputRasterDS

Set pOutputBand = pOutputBandCol.Item(0)

' QI RawPixel interface for input

Dim pInputRawPixel As IRawPixels

Set pInputRawPixel = pInputBand

' QI RawPixel interface for output

Dim pOutputRawPixel As IRawPixels

Set pOutputRawPixel = pOutputBand

' Create a DblPnt to hold the PixelBlock size

Dim pPnt As IPnt

Set pPnt = New DblPnt

pPnt.SetCoords pInputRasProps.Width, pInputRasProps.Height

' Creates PixelBlock of input with defined size

Dim pInputBlock As IPixelBlock

Set pInputBlock = pInputRawPixel.CreatePixelBlock(pPnt)

' Creates PixelBlock of output with defined size

Dim pOutputBlock As IPixelBlock

Set pOutputBlock = pOutputRawPixel.CreatePixelBlock(pPnt)

' Read input PixelBlock

pPnt.X = 0

pPnt.Y = 0

pInputRawPixel.Read pPnt, pInputBlock

' Get the SafeArray associated with the first band

' of output

Dim vSafeArray As Variant

vSafeArray = pOutputBlock.SafeArray(0)

' QI RasterProps for output for NoData handling

Dim pOutputRasProps As IRasterProps

Set pOutputRasProps = pOutputBand

' To make grid calculations other than the metric described at the

' beginning of this subroutine, modify the code for the equation below.

' Loop through the SafeArray and calculate each pixel value

' using equation given below

Dim I, J, S As Long

' Define a S value (delta s) for the metric calculation

S = 2

For I = 0 To pInputRasProps.Width - S

For J = 0 To pInputRasProps.Height - 1

If I - 1 >= 0 And I + 1 <= pInputRasProps.Width - 1 And J - 1 >= 0 _

And J + 1 <= pInputRasProps.Height - 1 Then

If Not CDbl(pInputBlock.GetVal(0, I, J)) = CDbl(pInputRasProps.NoDataValue) _

And Not CDbl(pInputBlock.GetVal(0, I + S, J)) = _

CDbl(pInputRasProps.NoDataValue) Then

vSafeArray(I, J) = 0.5 * (CDbl(pInputBlock.GetVal(0, I, J)) + _

CDbl(pInputBlock.GetVal(0, I + S, J))) * ((CDbl(pInputBlock.GetVal(0, I + S, J)) - _

CDbl(pInputBlock.GetVal(0, I, J))) / S)

Else

vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)

End If

Else

vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)

End If

Next J

Next I

' Write out the result

pOutputRawPixel.Write pPnt, pOutputBlock

End Sub