Search code examples
rrastervolume

Calculating the volume of a rasterpixel in R

I am performing a time series analysis on raster data (DEMs) from coastal regions. I want to measure the difference in volume (sand) between 2 different time stamps. I have already calculated the difference in height trough DEM differencing, but I do not know how to calculate the volume from that.

I have resampled the two rasters, so they now have the same resolution.

dem1 = resample(dem1, dem18)

DEM1

class      : RasterLayer 
dimensions : 5076, 6722, 34120872  (nrow, ncol, ncell)
resolution : 0.0464469, 0.0464469  (x, y)
extent     : 49584.86, 49897.07, 215276.8, 215512.6  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs 
source     : memory
names      : X20210223_DUDE_Oostende_T1_DEM 
values     : 3.852022, 19.46622  (min, max)

DEM18

class      : RasterLayer 
dimensions : 5076, 6722, 34120872  (nrow, ncol, ncell)
resolution : 0.0464469, 0.0464469  (x, y)
extent     : 49584.86, 49897.07, 215276.8, 215512.6  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs 
source     : 20220912_0216419_DUDE Oostende_T18_DEM.tif 
names      : X20220912_0216419_DUDE_Oostende_T18_DEM

I have subtracted both to measure the elevation difference.

elevation_difference = dem18 - dem1

I receive the plot below:

enter image description here

I now want to calculate the changes in volume of sand. Is is just elevation_difference * the size of the pixel?

Solution 
  • dem1  = raster("20210223_DUDE_Oostende_T1_DEM.tif")
dem18 = raster("20220912_0216419_DUDE Oostende_T18_DEM.tif")

dem1 = resample(dem1, dem18)

stacked = stack(dem1, dem18)
plot(stacked, xlim = x_lim, ylim = y_lim, col = terrain.colors(100))

    enter image description here

    # elevation difference

elevation_difference = dem18 - dem1
plot(elevation_difference, main = "Changes in height in m", col=terrain.colors(100), xlim = x_lim, ylim = y_lim)

# volume difference: height is in meter, pixels are in cm

area_pixel = prod(res(dem_stacked)) # the area of a cell is the product of the resolution in x and y
area_pixel # in cm²!
volume_difference = (elevation_difference$layer*100) * area_pixel
plot(volume_difference, main = "Changes in volume sand", col=terrain.colors(100), xlim = x_lim, ylim = y_lim)

    enter image description here enter image description here

    Normally the units should be correct, but better check it out :)