Section of the valley we want to compute volumes on (shown here is one of the levels)

This tutorial presents a method to calculate volume(s) of a valley.

Obviously that volume depends on the elevation we want to compute the volume from (see image above).

Let's say that valley bottom is at 100 m, top of surrounding hills is at 250 m and, more importantly, the lowest elevation between hills surrounding the valley is at 190 m. Now, obviously again, volumes can only be calculated from 100 m to 190 m. There is no point to compute a volume above 190 m because it will encompass a volume much larger than the true volume of the valley (think of filling the valley with water, once you reach 190 m, the water starts spilling over to surrounding valleys!).

WARNING! **The method presented below is not that much accurate**. Accuracy depends on the DEM cell size.

On the software side we need:

- QGIS 2.x
- A spreadsheet program (Excel here, but could be anything modern)
- Ability to write batch/script files (DOS is used here, but easy to adapt to UNIX)

On the data side we need (image below):

- a DTM/DEM raster (or any raster with elevations value) -
`DEM.tif`

- a polygon vector file that defines an area of interest (the region around the valley) -
`AOI.shp`

- a polygon vector file that defines the region of interest (the valley and surrounding hills, the ROI is a subset of the AOI) -
`ROI.shp`

The following steps are used to calculate volumes:

- Clip
`DEM.tif`

to`ROI.shp`

[Processing -> GDAL -> Clip raster by mask layer] ->`ROI.tif`

*Note that ROI.tif grid cell size is the principal factor in the accuracy of volume computation. The smaller the cell size is the greater the accuracy of the volume will be*. - Save
`AOI.shp`

as`AOI.CSV`

file. - Open it in Excel and add columns, named
`Znn`

, where 'nn' is the level value, one for each level, fill each column with the level value, see below:

- Import the CSV file in QGIS and save as shapefile ->
`AOI_pg.shp`

- Create a raster for each isosurface, save as tiff:
gdal_rasterize -a Znn -ot Float32 -tr 5 5 -a_nodata -9999 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -l AOI_pg AOI_pg.shp nn.tif

Here`-tr 5 5`

is the cell size (5x5 m) of`ROI.tif`

.*Again, the accuracy of the volume calculation is very much dependent on this value*. - From there clip the isosurface to the ROI, save as tiff:

gdalwarp -ot Float32 -q -of GTiff -dstnodata -9999 -cutline ROI_pg.shp -crop_to_cutline -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 nn.tif area-nn.tif

- Compute height differences with raster calculator: (DEM - area-nn), save as tiff:

gdal_calc --calc "((A-B)<0)*(A-B)*-1.0" --format GTiff --type Float32 --NoDataValue -9999 -A ROI.tif --A_band 1 -B Area-nn.tif --B_band 1--outfile Diff-area-nn.tif

- Import all
`Diff-area-nn.tif`

files into QGIS.

- For each of them set the
*NaN*to 0 and compute the histogram from the Properties window. Close the Properties window.

- Open the Properties window again and in the Metadata Tab look for the
__Sum of all cells__and__Cell Count__values, copy and paste into the excel spreadsheet (`Volumes.xlsx`

)

- On the first level, also copy the
__Pixel Size__into Excel

- Do the volume computation in Excel. Take care of the sign, volumes are generally positive!!!

Steps 5 to 7 are automated in `dosur.bat`

and `surfaces.bat`

`dosur.bat`

@echo off for %%A in (89 87 85 83 81 79) DO call surfaces.bat %%A

replace the values within the parenthesis with your level values.

`surfaces.bat`

@echo off echo Computing isosurface for level %1 if exist %1.tif del %1.tif gdal_rasterize -a Z%1 -ot Float32 -tr 5 5 -a_nodata -9999 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -l AOI_pg AOI_pg.shp %1.tif echo Computing clipped isosurface for level %1 if exist area-%1.tif del area-%1.tif gdalwarp -ot Float32 -q -of GTiff -dstnodata -9999 -cutline ROI_pg.shp -crop_to_cutline -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 %1.tif area-%1.tif echo Computing height differences for level: %1 if exist Diff-area-%1.tif del Diff-area-%1.tif gdal_calc --calc "((A-B)<0)*(A-B)*-1.0" --format GTiff --type Float32 --NoDataValue -9999 -A ROI.tif --A_band 1 -B area-%1.tif --B_band 1 --outfile Diff-area-%1.tif

For each grid cell of the isosurface that is above the DEM we calculate the height difference (h - step 7 above). We save all those height differences in the `Diff-area-nn.tif`

and we get the sum of all the heights from the raster metadata (step 10 above).

Now we multiply the cell size (x.y) by that value to get to the volume.

*Note that the larger (x.y) is the less accurate the calculated volume is.*

MS Excel is used to do volume calculation, but any spreadsheet program can be used.

The Volumes.xlsx sheet

There are only two formulas to compute the volume:

- Compute the cell area: G3*H3

- Compute the volume: D5*$F$3

and propagate this formula down to the other levels.

The Excel spreadsheet can be downloaded here: Volumes.xlsx

** Footnote**:

This method is obviously very rough but it gives far more plausible results than volumes derived from other algorithms.

I've tried SAGA grid volume algo in Processing and Golden Software Surfer algo on **my** dataset but volume values were not what I was expecting!

After *Binging* the web I've found this page that I've used to derive the method described above.

Published date: 4 April 2016.

No comments yet