Sediment Volume Project -- 07/23/08

Detailed instructions for analysis

--> Download the latest code. <--

Here are instructions for processing the Skate Creek subset of the Knowles basin area:

Processing other subsets of Knowles or Golden Ridge will involve similar steps, but may involve 'tuning' some of the program parameters for best effect. These will be discussed below.

Step 1: Get subset coordinates using PlotDEM.py

The entire Knowles DEM is too large to display full-size on my 4 Gb machine, so it must be decimated. Try the following command (replacing '../Data/' with the location of your file):

> PlotDEM.py ../Data/knwl_grd_asc.txt 0 0 1 10
Reading:  ../Data/knwl_grd_asc.txt
Header:  {'ncols': 10883.0, 'cellsize': 1.0, 'nrows': 5312.0, 'xllcorner': 435672.5, 'yllcorner': 4865015.5, 'NODATA_value': -9999.0}
Array shape:  (5312, 10883)
Valid elements:  43194462
Min = 64.8579, max = 560.357, mean = 289.827, std = 98.0393
Decimating by: 10
Array shape:  (532, 1089)
Valid elements:  431982
Min = 66.9065, max = 559.631, mean = 289.81, std = 98.0291

The full command syntax for PlotDEM.py is:

> PlotDEM.py 
Usage: PlotDEM.py data-file [min max mode decimate smooth output-file]
       data-file = binary or text file of 2d array
                   set to "demo" for demo mode
       min = data minimum
       max = data maximum (set both to 0 for auto-scaling, default)
       mode = -1 (no display), 0 (b & w), 1 (color, default)
       decimate = skip rows/colums in array (default 1)
       smooth = smooth array by square odd width window (default none)
       output-file = file to save binary result

The first thing you may want to do is convert the DEM to binary, which is slightly larger, but which loads much faster. You can do so with something like the following command. However, make sure that 1) you specify mode = -1 so there is no display, and 2) you do not decimate or smooth the file (decimate = 0, smooth = 0):

> PlotDEM.py ../Data/knwl_grd_asc.txt 0 0 -1 0 0 ../Data/kc_dem.dat
Reading:  ../Data/knwl_grd_asc.txt
Header:  {'ncols': 10883.0, 'cellsize': 1.0, 'nrows': 5312.0, 'xllcorner': 435672.5, 'yllcorner': 4865015.5, 'NODATA_value': -9999.0}
Array shape:  (5312, 10883)
Valid elements:  43194462
Min = 64.8579, max = 560.357, mean = 289.827, std = 98.0393
Saving to: ../Data/kc_dem.dat

To extract coordinates for Skate Creek, display the new binary DEM using a lower decimation value:

> PlotDEM.py ../Data/kc_dem.dat 0 0 0 4
Reading:  ../Data/kc_dem.dat
Loading from: ../Data/kc_dem.dat
Header:  {'ncols': 10883.0, 'cellsize': 1.0, 'nrows': 5312.0, 'xllcorner': 435672.5, 'yllcorner': 4865015.5, 'NODATA_value': -9999.0}
Array shape:  (5312, 10883)
Valid elements:  43194462
Min = 64.8579, max = 560.357, mean = 289.827, std = 98.0393
Decimating by: 4
Array shape:  (1328, 2721)
Valid elements:  2699618
Min = 66.4797, max = 560.277, mean = 289.825, std = 98.0393

Note that in using PlotDEM.py:

You can scroll the plot with the arrow keys.
You can drag the plot with the mouse.
Clicking the mouse prints the {x, y, z} coordinates under the crosshair.

To get coordinates, click the mouse at the lower-left corner and then the upper-right corner (or upper-left, lower-right, etc...) of a rectangle that encloses the area of interest:

There will be no visual feedback, but the coordinates of the two points will be printed:

   439296.500000   4866891.500000       474.669500
   440460.500000   4869739.500000       235.711500

Note that in selecting subsets:

You should select reaches which are fairly straight. If a reach changes direction, you should extract more than one subset, each fairly straight.
You should make sure that you also extract some area (a border) around the reach.

Step 2: Extract subsets using ExtractRegion.py

> ExtractRegion.py 
Usage: ExtractRegion.py dem-file br-file x1 y1 x2 y2 output-file
       data-file = text/binary file of 2d array
       br-file = text/binary file of bedrock points
       x1 y1 x2 y2 = coords of corners of rectangle to extract
       output-file = desired binary file of dem and br pts

Use copy/paste to enter the corner coords for subset extraction. Also enter the name of the appropriate bedrock point DEM. Points within the subset region will be attached to the output DEM:

> ExtractRegion.py ../Data/kc_dem.dat ../Data/skt_br_pts.txt 439296.5 4866891.5 440460.5 4869739.5 ../Data/SkateCreek.dat
Reading:  ../Data/kc_dem.dat
Loading from: ../Data/kc_dem.dat
Header:  {'ncols': 10883.0, 'cellsize': 1.0, 'nrows': 5312.0, 'xllcorner': 435672.5, 'yllcorner': 4865015.5, 'NODATA_value': -9999.0}
Array shape:  (5312, 10883)
Valid elements:  43194462
Min = 64.8579, max = 560.357, mean = 289.827, std = 98.0393
Reading:  ../Data/skt_br_pts.txt
Header:  {'ncols': 211.0, 'cellsize': 1.0, 'nrows': 973.0, 'xllcorner': 439856.99143266003, 'yllcorner': 4868462.5335611999, 'NODATA_value': -9999.0}
Array shape:  (973, 211)
Valid elements:  25
Min = 0, max = 24, mean = 12, std = 7.2111
25 points extracted
Extracting 2849x1165 subset from:
   i1 =  3624 x1 =    439296.500000
   j1 =   587 y2 =   4869739.500000
   i2 =  4788 x2 =    440460.500000
   j2 =  3435 y1 =   4866891.500000
copied 25 points
Saving to: ../Data/SkateCreek.dat

The extracted subset and attached bedrock points can be viewed with PlotDEM.py:

> PlotDEM.py ../Data/SkateCreek.dat 0 0 0 3
Reading:  ../Data/SkateCreek.dat
Loading from: ../Data/SkateCreek.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.1021, max = 487.853, mean = 247.126, std = 87.3216
Decimating by: 3
Array shape:  (950, 389)
Valid elements:  369550
Min = 73.6712, max = 487.624, mean = 247.125, std = 87.3446

Note:

There is also a program called PlotBR.py which can be used to combine the knls_basin_br.txt, skt_br_pts.txt, and lknl_br_asc.txt files into a single binary file, although they can be used separately by ExtractRegion.py:

Step 3: Smooth the subset DEM using PlotDEM.py

Smooth the subset with a 5x5 m square filter, and save it in a file ending in '_s.dat'. This file will then be used, along with a de-trended file (see below) to extract valley wall coordinates. The smoothing operation may take a few minutes to complete:

> PlotDEM.py ../Data/SkateCreek.dat 0 0 -1 0 5 ../Data/SkateCreek_s.dat 
Reading:  ../Data/SkateCreek.dat
Loading from: ../Data/SkateCreek.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.1021, max = 487.853, mean = 247.126, std = 87.3216
Smoothing array by:  5
Min = 73.9242, max = 487.585, mean = 247.126, std = 87.3128
Saving to: ../Data/SkateCreek_s.dat

Step 4: Remove the trend from the subset DEM using RemoveTrend.py

> RemoveTrend.py
Usage: RemoveTrend.py dem-file x1 y1 x2 y2 output-file
       dem-file = text/binary file of 2d array
       x1 y1 x2 y2 = coords of points defining trend
       output-file = desired binary file

Get the coordinates of two points defining the trend using PlotDEM.py. The points should lie at opposite ends of the stream bed, one at low elevation, the other at high elevation. For example, two points in the Skate Creek subset might be indicated as below:

Save the de-trended file with an ending of '_t.dat'. Be sure to use the smoothed file ('_s.dat') as input:

> RemoveTrend.py ../Data/SkateCreek_s.dat 440019.5 4869318.5 439686.5 4867851.5 ../Data/SkateCreek_t.dat
Reading:  ../Data/SkateCreek_s.dat
Loading from: ../Data/SkateCreek_s.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.9242, max = 487.585, mean = 247.126, std = 87.3128
Removing trend...
Trend: z1=84.437996, z2=176.105908, distance=1504.319780, slope=0.060936, vector=[-0.22136251 -0.97519159]
Saving to: ../Data/SkateCreek_t.dat

The resulting file should make it much easier to see the valley at higher elevations:

> PlotDEM.py ../Data/SkateCreek_t.dat 0 0 0 3                              
Reading:  ../Data/SkateCreek_t.dat
Loading from: ../Data/SkateCreek_t.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 0, max = 271.363, mean = 120.397, std = 58.1127
Decimating by: 3
Array shape:  (950, 389)
Valid elements:  369550
Min = 0.184239, max = 271.307, mean = 120.426, std = 58.1233

Note two things:

The de-trended DEM is used only as a visual aid for entering an initial centerline in the program that follows. All other programs use as data input the smoothed '_s.dat' version of the DEM.
Although you can decimate the DEM files in most programs, a full-size DEM with no decimation should be used for best results in the programs that follow.

Step 5: Estimate the valley centerline and boundaries using GetValley.py

> GetValley.py
Usage: GetValley.py data-file-base-name [decimate]
       data-file-base-name = common part of dem file names
       decimate = skip rows/colums in array (default 1)

GetValley.py reads the '_s.dat' and '_t.dat' versions of DEM files. It uses the de-trended version for viewing, and the smoothed version for calculations. For best results, don't decimate these files unless you run out of memory using the full-size files.

> GetValley.py ../Data/SkateCreek
Reading:  ../Data/SkateCreek_t.dat
Loading from: ../Data/SkateCreek_t.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 0, max = 271.363, mean = 120.397, std = 58.1127
Reading:  ../Data/SkateCreek_s.dat
Loading from: ../Data/SkateCreek_s.dat
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.9242, max = 487.585, mean = 247.126, std = 87.3128

Procedure:

Click the mouse to add new points in sequential order from the upstream to the downstream end of the valley.
Use the arrow keys (not the mouse) to move the image.
Points cannot be removed or edited, but you can start over again by choosing the 'Operations->New curve' menu item.
Points can be entered at any spacing, but will be interpolated after the curve is ended. Use just enough detail to stay in mid-valley, but not enough to follow every 'wiggle'. The centerline is used only to calculate profiles and their orientations, and it is best not to make the centerline too sinuous, or profiles will intersect one another.
It is best to enter 3 closely spaced colinear points at the start and end of the centerline, so that neither end has significant curvature.
Do not proceed from a tributary into the main channel (e.g. from Skate Creek into Lower Knowles) or 'bad things' will happen to the end of the valley. It is best to stop slightly before the valley widens into an intersection, although you may need to experiment a few times before getting it right.

Once you have entered the points, choose the 'Operations->End curve' menu item. The centerline will be interpolated, profiles will be calculated perpendicular to the centerline, derivatives on those profiles will be calculated, and left and right valley boundaries will be calculated at distances from the centerline where the slope of each wall first reaches some fraction (default 67%) of the maximum slope for the profile (calculated separately for left and right walls):

You will also see plots of the profiles (including thalweg in yellow and wall boundaries in white) and the derivatives along each profile:

At this point you could stop and save the results, or you can iterate to generate a new centerline and set of profiles and boundary curves. To do so, choose the 'Operations->Iterate' menu item. This will generate a new centerline midway between the left and right boundary curves, and then a new set of profiles perpendicular to the centerline, and then new left and right boundary curves:

New profile and derivative plots will also be shown:

When you like what you see, choose the 'Operations->Write file' menu item. This will interpolate the left and right walls, attach all three curves to the DEM, and write an output file having the base name plus '_v.dat'. You can then view that file using PlotDEM.py:

> PlotDEM.py ../Data/SkateCreek_v.dat 0 0 0 3
Reading:  ../Data/SkateCreek_v.dat
Loading from: ../Data/SkateCreek_v.dat
center len = 195, left slope len = 195, right slope len = 195
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.9242, max = 487.585, mean = 247.126, std = 87.3128
Decimating by: 3
Array shape:  (950, 389)
Valid elements:  369550
Min = 74.0283, max = 487.576, mean = 247.125, std = 87.3356

Tips for use:

It may take several iterations for you to like the calculated centerline and boundary curves.
Unfortunately, you can't undo an iteration. You must start over by choosing 'Operations->New curve'.
Pay careful attention to both ends of the valley. Each end should be fairly straight, and all three curves should start and end together without any significant change in direction. Under no circumstances should any curve reverse direction.
Often, if after an iteration the curves don't behave well at the ends, performing another iteration may fix the problem. If several iterations do not improve the geometry, you may need to start over.

There are several 'tunable parameters' in the program which you can change to get different results. These are actually in the Plots.py file, in the __init__() routine of DEMPlotPanel2:

		# adjustable parameters
		self.profileHalfWidth = 100
		self.slopeFraction = 0.667
		self.interpolationLength = 10.0
		self.smoothLength = 5
		self.leftRightDamping = 1
		self.centerDamping = 1

profileHalfWidth sets the distance left and right of the centerline that profiles are calculated. Must be an integer. The valley can't be any wider than 2x this value.
slopeFraction is the (decimal) fraction of the maximum slope along the left or right side of the profile that the boundaries are placed. A smaller/larger value will make the resulting valley narrower/wider.
interpolationLength is the resolution (distance) to which all curves are interpolated, using circular arcs, and is the resulting distance between profiles.
smoothLength is the number of points (not distance) which are used to smooth each curve. Must be an integer. For example, smoothLength = 5 averages 2 points before and after every individual point to perform smoothing.
leftRightDamping is the decimation value for calculating the left and right boundaries from the profiles, at the currently set interpolationLength. Must be an integer. For example, leftRightDamping = 2 uses every other profile to calculate left and right curves. This has the effect of damping the variation in each curve, and a larger value of leftRightDamping will make the valley width more constant along its length. The combined effects of interpolationLength, smoothLength, and leftRightDamping are to act as a low-pass filter on all curves.
centerDamping controls the decimation of left and right points used to calculate a new centerline. Must be an integer. Since setting leftRightDamping will also result in damping of the centerline, centerDamping provides an alternate way to control centerline sinuosity independently of left and right boundary variation.

You will have to play with these settings to see how they affect the resulting valleys, and which combination of settings you prefer. At present there is no way to alter these settings without restarting the program, so their effect on a single set of input points cannot be observed. Two strategies for allowing this might be 1) create a dialog and corresponding menu for changing settings during a run, or 2) add functionality for reading a set of initial centerline points from a file so that the same set of points can be used for multiple runs of the program.

Here is the effect of redigitizing the Skate Creek valley using a leftRightDamping value of 3:

Step 6: Estimate volumes using PlotValley2.py

> PlotValley2.py 
Usage: PlotValley2.py valley-dem [output-flag: 0 = text DEM, 1 = binary DEM (default)]

To estimate volumes along the Skate Creek valley produced in step 5:

> PlotValley2.py ../Data/SkateCreek_v.dat
Reading:  ../Data/SkateCreek_v.dat
Loading from: ../Data/SkateCreek_v.dat
center len = 195, left slope len = 195, right slope len = 195
Header:  {'ncols': 1165.0, 'cellsize': 1.0, 'nrows': 2849.0, 'xllcorner': 439296.5, 'yllcorner': 4866891.5, 'NODATA_value': -9999.0}
Array shape:  (2849, 1165)
Valid elements:  3319085
Min = 73.9242, max = 487.585, mean = 247.126, std = 87.3128
Read 3 curves:
   left: 202 pts, 1993.36 m
   right: 199 pts, 1968.73 m
   center: 195 pts, 1934.42 m
24 bedrock pts
Calculating path vectors...
195 vectors
Calculating wall intersections...
Calculating bedrock coordinates and elevations...
...

Full output

Here is the valley geometry:

Here is the valley elevation:

(Click for larger image)

Here are the valley profiles:

Here is the resulting thickness DEM:

(Click for larger image)

Note the following about the volume approximation:

The surface elevation can vary in an arbitrary manner from point to point.
The bedrock elevation consists of a level, horizontal, base that varies only in s (the down-valley coordinate) near the middle of the valley, plus a set of beveled edges that can vary in both s and n (the cross-valley coordinate) in a linear way (yielding a plane in 3d). The slope of the base wrt n is averaged to a single value for each step (although separate values for left and right edges) based on the actual wall slopes for each profile. Therefore, while the thickness near the middle of the valley should be fairly continuous across steps, it will in general be discontinuous across steps at the edges (depending on whether the surface and base are moving toward or away from each other along s). This is a failure of the simplification used to model the beveled bedrock base at the edges of the valley. Ideally, a more continuous model of the bedrock elevation should be used which avoids the discontinuities at the edges of the steps.

Here is a valley and thickness DEM generated for Lower Knowles, using a leftRightDamping value of 3: