(Click for larger image)
--> Download the latest code. <--
Last weekend I did some initial volume calculations, using the PlotValley.py program. Most of the enhancements to the code have been made in that file. The program now sorts bedrock points wrt downstream distance, and calculates both dz and dz/ds, and min, max, mean and std drop in elevation:
Calculating bedrock coordinates and elevations...
x y z s dz dz/ds
---------------- ---------------- ---------------- ---------------- ---------------- ----------------
441482.107862 4868617.690507 -9999.000000 -23.134745
441493.107862 4868701.690507 110.551156 60.073256 10109.551156 121.497343
441449.107862 4868912.690507 106.580796 297.214788 -3.970360 -0.016743
441438.107862 4868942.690507 104.933784 327.084366 -1.647012 -0.055140
441430.107862 4868962.690507 106.026768 347.166017 1.092984 0.054427
441404.107862 4868968.690507 105.174436 369.301664 -0.852332 -0.038505
441240.107862 4869065.690507 102.137268 562.170643 -3.037168 -0.015747
441182.107862 4869119.690507 100.126252 642.149938 -2.011016 -0.025144
441089.107862 4869155.690507 97.796675 737.935844 -2.329577 -0.024321
440931.107862 4869282.690507 96.028834 941.271280 -1.767841 -0.008694
...
Calculating thalweg coordinates and elevations...
x y z s dz dz/ds
---------------- ---------------- ---------------- ---------------- ---------------- ----------------
441484.962117 4868641.849659 111.928576 0.000000
441499.696638 4868658.362687 111.523200 19.995710 -0.405376 -0.020273
441498.681231 4868679.803017 110.464328 39.993051 -1.058872 -0.052951
441490.145996 4868702.269767 110.327356 59.992682 -0.136972 -0.006849
441491.984597 4868723.002980 109.724272 79.991224 -0.603084 -0.030156
...
438909.430426 4870197.842656 68.162886 3298.062108 -0.174478 -0.008723
438911.514459 4870220.263315 67.703177 3318.044051 -0.459709 -0.023006
438906.286759 4870241.165055 68.308993 3338.039228 0.605816 0.030298
438906.097914 4870262.200873 68.007907 3358.037477 -0.301086 -0.015056
438895.743404 4870284.255978 68.336619 3378.037047 0.328712 0.016436
dz: 169 -2.9229632 1.3470756 -0.257940572781 0.734904083293
The bedrock base elevation is then calculated as z(base) ≤ z(thalweg) - 0.5 and z(base) ≤ z(bedrock) and either slope(base) ≤ 0 or slope(base) < 0. Here is output for the latter (strictly descending base slope):
Calculating composite base elevations...
s z dz dz/ds
---------------- ---------------- ---------------- ----------------
0.000000 111.428576
19.995710 111.023200 -0.405376 -0.020273
39.993051 109.964328 -1.058872 -0.052951
59.992682 109.827356 -0.136972 -0.006849
79.991224 109.224272 -0.603084 -0.030156
...
3198.107066 68.588928 -0.132324 -0.002206
3238.095112 68.027490 -0.561438 -0.014040
3278.060737 67.837365 -0.190126 -0.004757
3298.062108 67.662886 -0.174478 -0.008723
3318.044051 67.203177 -0.459709 -0.023006
dz: 85 -2.037336 -0.0135904 -0.520298809412 0.431279311578
Base elevations are assigned to each step (and can be assigned to any point based on s):
Getting step elevations...
s z
---------------- ----------------
0.000000 111.428576
19.995710 111.023200
39.993051 109.964328
59.992682 109.827356
79.991224 109.224272
...
3298.062108 67.662886
3318.044051 67.203177
3338.039228 67.203177
3358.037477 67.203177
3378.037047 67.203177
Cell coordinates are calculated from steps and sidewall segments. Cells are polygons that consist of the left and right boundary step intersections, and any boundary points between the steps (usually about 6 points total), accumulated in a CCW direction beginning with the left intersection of the 'upper' step (lower s value). Steps which cross are removed from consideration (in the case of LKC digitized at 20 m, only 1 step was removed, at the turnoff for Bear Creek). Here are printouts of the cell coordinates {x, y, z_top, z_bottom, s}, the area of the rectangular bound of the cell, the area of all prismatic columns within a cell (simple square cross sections for {x, y} points within the boundary of the cell), the ratio of the two, the average thickness of the cell (lengths of prisms), and the volume of all prisms. Note that this volume is 'crude' in the sense that it is a bundle of vertical square prisms within the cell, each having a length delimited by the DEM elevation and bedrock base at that point, and horizontal faces with area equal to the DEM cell size squared:
Getting cell coordinates and volumes... [[ 4.41455943e+05 4.86865275e+06 1.18272660e+02 1.11428576e+02 0.00000000e+00] [ 4.41510321e+05 4.86863232e+06 1.14715376e+02 1.11428576e+02 0.00000000e+00] [ 4.41516914e+05 4.86864845e+06 1.16437928e+02 1.11105764e+02 1.59231171e+01] [ 4.41517937e+05 4.86865280e+06 1.16807588e+02 1.11023200e+02 1.99957098e+01] [ 4.41462394e+05 4.86866974e+06 1.15914544e+02 1.11023200e+02 1.99957098e+01] [ 4.41456614e+05 4.86865521e+06 1.18189104e+02 1.11371849e+02 2.79814945e+00]] bounds area = 9279.102276, cell area = 4668.000000, % = 0.503066, av thickness = 1.754927 vol = 8191.998339 [[ 4.41462394e+05 4.86866974e+06 1.15914544e+02 1.11023200e+02 1.99957098e+01] [ 4.41517937e+05 4.86865280e+06 1.16807588e+02 1.11023200e+02 1.99957098e+01] [ 4.41521490e+05 4.86866790e+06 1.17789812e+02 1.10268239e+02 3.42535332e+01] [ 4.41521972e+05 4.86867419e+06 1.17197660e+02 1.09964328e+02 3.99930510e+01] [ 4.41468544e+05 4.86868707e+06 1.17465868e+02 1.09964328e+02 3.99930510e+01] [ 4.41464009e+05 4.86867380e+06 1.15164852e+02 1.10771311e+02 2.47527655e+01]] bounds area = 8166.369134, cell area = 4548.000000, % = 0.556918, av thickness = 1.928215 vol = 8769.523623 ... [[ 4.38885311e+05 4.87024015e+06 7.57350980e+01 6.72031772e+01 3.33803923e+03] [ 4.38962715e+05 4.87024389e+06 7.58341184e+01 6.72031772e+01 3.33803923e+03] [ 4.38962507e+05 4.87024913e+06 7.61742704e+01 6.72031772e+01 3.34395167e+03] [ 4.38962976e+05 4.87026163e+06 7.63792612e+01 6.72031772e+01 3.35803748e+03] [ 4.38885099e+05 4.87026241e+06 7.63684152e+01 6.72031772e+01 3.35803748e+03] [ 4.38885299e+05 4.87024786e+06 7.65568944e+01 6.72031772e+01 3.34496288e+03]] bounds area = 6934.426125, cell area = 6136.000000, % = 0.884861, av thickness = 3.684001 vol = 22605.027482 [[ 4.38885099e+05 4.87026241e+06 7.63684152e+01 6.72031772e+01 3.35803748e+03] [ 4.38962976e+05 4.87026163e+06 7.63792612e+01 6.72031772e+01 3.35803748e+03] [ 4.38963259e+05 4.87026915e+06 7.61208284e+01 6.72031772e+01 3.36671734e+03] [ 4.38965389e+05 4.87027882e+06 7.75899496e+01 6.72031772e+01 3.37803705e+03] [ 4.38882783e+05 4.87028527e+06 7.77915380e+01 6.72031772e+01 3.37803705e+03] [ 4.38885024e+05 4.87026785e+06 7.66422988e+01 6.72031772e+01 3.36284069e+03]] bounds area = 7810.716614, cell area = 6344.000000, % = 0.812217, av thickness = 3.620681 vol = 22969.601227 Total volume: 5.591071e+06 Average width: 82.787177 Average thickness: 4.834486
Here is the full output of this run.
Here is a plot of the valley elevation along its length. The orange curve at bottom is the bedrock base, and is 1) strictly descending, 2) at or below all designated bedrock points, 3) at least 0.5 m below the thalweg:
(Click for larger image)
Here are plots of the valley cross-section profiles. At left are the 'raw' elevations. At center are only those profiles which are descending. At right are descending profiles smoothed by a 5 m filter. The meandering thalweg is in white. The colors mean nothing, and are only used to help see the different profiles. The centerline is at center. Left directions are to left, and negative of center:
(Click for larger image)
It would seem from these plots that I am over-estimating the valley width. The slope of the left wall is clearly visible in the last plot, although the right wall is less distinguishable. The thalweg goes right up to both walls, and probably indicates the actual limits of the valley floor. Therefore, I would expect the volumes estimated above to be too high, as they include part of the valley walls (or at least part of the left wall).
Here is a plot of the descending steps in the valley, with colors the same as in the profiles above. All other steps initially ascended, although their base elevations have since been reset according to the orange plot above:
(Click for larger image)
Here are the outlines of all the cells:
(Click for larger image)
Here are the volumes wrt stream length:
The volume estimate calculated above (5.591071e+06 m^3) seems too large. An estimate based on stream length, average width, and average thickness yields 3387.05 x 82.79 x 4.83 = 1354398.99 m^3, or only about 1/4 of that above. One explanation for this might be that the wider cells were also thicker. To see if this was the case, I cross-plotted thickness and width:
Indeed, there does seem to be a positive correlation between thickness and width, possibly because the wider cells included too much of the steeper valley walls (but also maybe not).
Well, this is just an initial calculation. Refinements might include:
There is also the question of how to perform the analysis on the entire KC grid area at once, rather than stream by stream, and how to automatically locate the valley bounds without digitizing by hand. Finally, another question which arises is whether you want me to do more work on the project in July, or if you are content to take what I've done so far and continue on your own. Please let me know by this Friday!