How to measure volume in raster?

There are many similar questions like this. But I haven't found a usefull answer yet. So I hope some of you have a suggestion to this.

I've done some flyovers with the drone at a contruction area, and then got a DSM out of these, using Pix4D. Now I have 2 different DSM, and then imported these into QGIS. Through rastercalculator i've withdrawn these from each other to get the height difference. See below: I then would like to find out the volume of where there have been a height difference.. So some of the area I have marked red (increased volume) and the blue (decrease in volume). Also using rastercalculator and polygonize. See below.  Is there a smart way to find the volume of the blue and red marks of these rasterlayers in QGIS, or do any of you have an alternative program that would do the job easier and faster?

• You could have a look at this example gis.stackexchange.com/questions/258482/…. In your case you could skip the steps related with produce the surface of the previous contidion Dec 15 '17 at 11:45
• I have an open ticket in GRASS to make it esier with r.volume: trac.osgeo.org/grass/ticket/3442#no1 Dec 15 '17 at 11:48
• Magnus -- I saw your partial answer (now gone) and if I understand what you're getting at: Each cell of your raster has a size (x meters^2) and each cell has a value (z meters). So can't you just multiply each raster cell value by its size and get meters^3, then sum the values of those calculations to derive total volume? -- I'd post this as an answer but it seems too simple so I'm not sure what you are really asking. Dec 15 '17 at 15:42

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.

• Hi Kantan, I've been told to leave a comment instead of an answer.. Thanks for your reply. I've got a sum-value of 11758 for an specific red area. I'm not sure what that value shows? My pixelsize for the raster is 0.1 x 0.1. Now, i've checked the same pile of gravel in pix4D, to see what it calculate in m3. And it says about 2300 m3.. So, according to you, it's the 11758 x 0.1 = 1175.8.... Thats not close to pix4D's 2300.. Or am I doing it wrong? Dec 20 '17 at 12:30
• Hi Magnus. The formula to compute the total volume is the sum of the volume difference for all cells. For each cell, the volume is the area of the cell times the difference in height. As the area is constant for all pixel, you can take it out and compute the total difference in the raster values, that you then multiply by the cell area. Here there is a small mistake, it should be 11758x0.01. That is not even close to the Pix4D value yet. How did you compute the volume in pix4D? Could you detail your process in detail to see if you might have forgotten something? Dec 20 '17 at 13:21
• Ahh... Hi again Kantan. I've been sitting and trying it back- and forwards now. The mistake was that, when I use rastercalculator with 2 rasterlayers, both with 0.1 x 0.1, the output somehow come out 0.4 x 0.4... don't know why.. It's not the exact same as Pix4D in this example above, but I've tried it with some other piles, and that's quite accurate. Thanks a million!! Dec 20 '17 at 20:26
• Just a quick one.. I made some vector layers, to calculate the piles by.. But when i'm editing, or adding new polygons in the same vector layer, the calculated zonal statistic and my m3-calculation (done by field calculator) won't update the numbers in the Attribute tabel. I do remember to save it before checking, alså tried to update in general. Is there a reason for that? Dec 20 '17 at 21:10
• Glad that you made it work! If you use QGis Raster calculator, you might want to select your first pix4D layer and click "current layer extend" on the right panel. That might help. About the field calculator question, I'm not really using them as I wasn't sure about the way the values were refreshed. I'd recommend you open a new question for this specific issue, where more knowledgeable people could answer. I'd be interested in the answer too! Finally, don't forget to accept the answer, so other people know it worked for you! Welcome to GIS SE! Dec 21 '17 at 8:16

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")