How to measure volume in raster?

The best way I've found to do so is to use a combination of the raster calculator and zonal statistics.

You first need to calculate the difference between the two rasters, using a filter to get either the positive or the negative difference, and no values elsewhere. In Qgis calculator this would give something like :

((A-B)>0)*(A-B) To get the positive change

((A-B)<0)*(A-B)*-1 To get the negative change (multiply by -1 if you're interested in the absolute value of the change)

Once this is done, create a new vector layer, showing the zone of interest (blue + red, just avoid the places where you are not interested. As this is construction, you can probably select only the zones where you know there has been changes).

Then use the tool in Raster -> Zonal Statistics to read the sum of the values of all pixels within the vector layer. Do it twice, both for positive and negative change.

Multiply then the value obtained for each case with the size of your pixel squared, to get the total volume.

A very fast way to get the volume with some python code and gdal2xyz. You have to calculate your difference raster only. Don't care about positiv or negativ value, we will distinguish between them later.

1. Convert your difference raster file to a csv file with the OSGeo4W Shell:

   gdal2xyz your_raster.tif your_raster.csv
   

2. Create a numpy array from this csv file and create two subsets:

   import numpy as np
   x = np.genfromtxt (r'your_raster.csv', delimiter=" ", usecols=range(2,3))
   pos = x[x>0.2]
   neg = x[x<-0.2]
   pos.sum(axis=0)
   neg.sum(axis=0)
   

Explanation of the Python code:

The csv file usually has 3 columns x,y,z. With usecols= range(2,3) you import only the 3rd column (values of difference).

With pos = x[x>0.2] you save all positive values >0.2 to a new array.

Afterwards you can sum the two subsets. Use pos.sum(axis=0) to get the sum of all positive values >0.2.

If your raster has NODATA (-99999 for example) values you can remove them from your array with:

index = np.argwhere(x==-99999)
x = np.delete(x, index)

Then you can go further with creating the subsets.

There is also a way to create a numpy array from raster files directly. I get slightly different results with this approach:

import numpy as np
from osgeo import gdal
ds = gdal.Open("your_raster.tif")
x = np.array(ds.GetRasterBand(1).ReadAsArray())
pos = x[x>0]
neg = x[x<0]
pos.sum(axis=0)
neg.sum(axis=0)

Finally you have to multiply your sum with the cellsize you have.