Calculating the volume between two interpolation layers in QGIS?

Question

Given data

From a shapefile with ca. 80 points (displaying borewells) I extracted two interpolations:

1. interpolation values of the different borewells’s elevations;
2. interpolation values for the elevation of the bottom of the first geological layer.

Desired outcome

I need to know the volume of the first geological layer in the area of interest.

I am a beginner and not familiar with GRASS. It would be great if you could explain everything step by step, maybe even providing some screenshots.

Answer 1

The first step is to define a grid resolution (g.region) and then convert your points to rasters with v.surf.idw. Check that the rasters look OK. Sometimes interpolating introduced artifacts. Then, subtract the lower from the upper layer to get thickness. This could be done with r.mapcalc "thick = top - bot". Finally, the volume is the thickness times the x resolution times the y resolution.

This answer doesn't explain everything nor use screenshots (I use the CLI not the GUI, so don't have any screenshots to provide), but it might help you get started. If you can use this to break down the problem to smaller tractable issues, you should be able to figure them out, or post more specific questions here that will be easier to answer.

Edit: I just realized that your title says QGIS. I noticed the grass tag and answered about how to do this in GRASS.

Answer 2

Assuming that borewells’s elevations (in my example a constant raster with value 2000) and elevation of the bottom of the first geological layer (in my example a DEM raster) are two rasters with the same resolution and extension (both projected in meters), next PyQGIS code was tested to get the volume in km³.

registry  = QgsMapLayerRegistry.instance()

bore_elev = registry.mapLayersByName('borewells_elevations')
first_geol = registry.mapLayersByName('first geological layer')

extent = bore_elev[0].extent()

rows = bore_elev[0].height()
columns = bore_elev[0].width()

xsize = bore_elev[0].rasterUnitsPerPixelX()
ysize = bore_elev[0].rasterUnitsPerPixelY()

cell_area = xsize*ysize

block1 = bore_elev[0].dataProvider().block(1, extent, columns, rows)
block2 = first_geol[0].dataProvider().block(1, extent, columns, rows)

sum = 0

print "Processing..."

for i in range(rows):
    for j in range(columns):
        sum += (block1.value(i,j) - block2.value(i,j))

print "volume:{:.2f} km3".format(sum*cell_area/1e9)

After running it at the Python Console of QGIS, I got:

>>>execfile(u'/home/zeito/pyqgis_scripts/calculating_volume.py'.encode('UTF-8'))
Processing...
volume:1275.13 km3