# [R] Calculating volume under polygons

Micha Silver t@v|b@r @end|ng |rom gm@||@com
Tue Nov 21 17:53:44 CET 2023

```Again, please keep communication on the maillist, to help others as well.

The R spatial maillist is R-sig-Geo in case you want to continue there.

Here's my solution. I read in your layers, then looped over all
polygons. Inside the loop I extract the minimum value of the DEM, then
subtract that from the original values, to get the height above minimum.
And finally, using extract_exact() I get the sum of all these height
pixels within the polygon = volume.

After the loop, rbind to merge the new polygons with volume column.

Here's the code I tried:

# Load packages and Setup directories
pkgs <- c('terra', 'exactextractr', 'dplyr', 'sf')
invisible(lapply(pkgs, require, character.only = TRUE))
GIS_dir <- "./GIS"
dem_file <- file.path(GIS_dir, "Base1.tif")
poly_path <- file.path(GIS_dir, "Sites_Sarch.shp")

dem <- rast(dem_file)

# Start a loop over all polygons, and extract data from each
volumes_list <- lapply(1:length(polys\$geometry), function(x) {
poly <- polys\$geometry[x] # A single sfc object
poly <- st_as_sf(poly)    # coerced to an sf object
poly_id <- polys\$Name[x]
# Get minimum value in this polygon
dem_min <- min(values(dem_crop), na.rm = TRUE)
# Make a new cropped raster
# that is the height above minimum
dem_height <- dem_crop - dem_min

# Get the sum of the values in this dem_height
volume <- exactextractr::exact_extract(dem_height,
poly, 'sum') # in m^3
# Make a data.frame of all necessary values:
# minimum and sum for this poly
# and the id of the poly
poly_sf <- data.frame("Name" = poly_id,
"poly_min" = dem_min,
"volume" = volume)
poly_sf <- st_set_geometry(poly_sf, st_geometry(poly))

return(poly_sf)
})

# Bind the list output into a merged sf object
result_polys <- do.call(rbind, volumes_list)
print(st_drop_geometry(result_polys))

HTH

On 21/11/2023 10:33, javad bayat wrote:
> Micha;
> I finally calculated the volumes. But it needs too many codes. Is
> there any way to reduce the amount of codes?
> Do you have any idea regarding how to calculate the exact volume under
> the polygon?
>
>
> x1 = as.data.frame(x[1])
> x2 = as.data.frame(x[2])
> x3 = as.data.frame(x[3])
> x4 = as.data.frame(x[4])
> x5 = as.data.frame(x[5])
> x6 = as.data.frame(x[6])
> x7 = as.data.frame(x[7])
> x8 = as.data.frame(x[8])
> x9 = as.data.frame(x[9])
> x10 = as.data.frame(x[10])
> x11 = as.data.frame(x[11])
> x12 = as.data.frame(x[12])
>
> x1\$Depth = x1[,1] - min(x1[,1])
> x2\$Depth = x2[,1] - min(x2[,1])
> x3\$Depth = x3[,1] - min(x3[,1])
> x4\$Depth = x4[,1] - min(x4[,1])
> x5\$Depth = x5[,1] - min(x5[,1])
> x6\$Depth = x6[,1] - min(x6[,1])
> x7\$Depth = x7[,1] - min(x7[,1])
> x8\$Depth = x8[,1] - min(x8[,1])
> x9\$Depth = x9[,1] - min(x9[,1])
> x10\$Depth = x10[,1] - min(x10[,1])
> x11\$Depth = x11[,1] - min(x11[,1])
> x12\$Depth = x12[,1] - min(x12[,1])
>
>
> x1\$Vol = x1[,2] * x1[,3]
> x2\$Vol = x2[,2] * x2[,3]
> x3\$Vol = x3[,2] * x3[,3]
> x4\$Vol = x4[,2] * x4[,3]
> x5\$Vol = x5[,2] * x5[,3]
> x6\$Vol = x6[,2] * x6[,3]
> x7\$Vol = x7[,2] * x7[,3]
> x8\$Vol = x8[,2] * x8[,3]
> x9\$Vol = x9[,2] * x9[,3]
> x10\$Vol = x10[,2] * x10[,3]
> x11\$Vol = x11[,2] * x11[,3]
> x12\$Vol = x12[,2] * x12[,3]
>
> Volume = data.frame(ID = c(1:12), Vol = c(sum(x1\$Vol),
>                                   sum(x2\$Vol),
>                                   sum(x3\$Vol),
>                                   sum(x4\$Vol),
>                                   sum(x5\$Vol),
>                                   sum(x6\$Vol),
>                                   sum(x7\$Vol),
>                                   sum(x8\$Vol),
>                                   sum(x9\$Vol),
>                                   sum(x10\$Vol),
>                                   sum(x11\$Vol),
>                                   sum(x12\$Vol)))
> Volume
>
> On Tue, Nov 21, 2023 at 11:14 AM javad bayat <j.bayat194 using gmail.com> wrote:
>
>     Dear Micha, I have sent my question to the R-sig-Geosite too. To
>     clarify more, I am trying to calculate the volume of polygons by a
>     DEM raster. The shapefile has several polygons (12) and I want to
>     calculate the volume of the polygons based on their minimum
>     elevations. I tried these codes too, but I don't know how to
>     calculate the volumes at the end. The function I wrote uses the
>     minimum elevation of all polygons, not each polygon minimum.
>     You know I want to calculate the volume of the soil to be removed
>     on these sites. Maybe using the minimum elevation is not
>     appropriate for this work, maybe calculating the relation between
>     elevation and slope of the area under the polygons will act
>     better. I dont know how to do it.
>     I would be more than happy if you help me.
>     The coordinate system is UTM zone 40N.
>     I have uploaded the files.
>
>     Sincerely yours
>
>     ################################################################
>     library(terra)
>     library(exactextractr)
>     library(dplyr)
>     library(sf)
>
>     # Read the DEM raster file
>     r <- rast("E:/Base1.tif")
>     # Read the polygon shapefile that contains several polygons from
>     different sites
>     crop_r <- crop(r, extent(p))
>     # Extract the area of each cell that is contained within each polygon
>     x <- exact_extract(r, p, coverage_area = TRUE)
>     # Add polygon names that the results will be grouped by
>     names(x) <- p\$Name
>     # Bind the list output into a data frame and calculate the
>     proportion cover for each category
>     result <- bind_rows(x, .id = "Name") %>%
>       group_by(Name) %>% summarize(Area = sum(coverage_area)) %>%
>       group_by(Name) %>% mutate(Volume = Area * min(na.omit(a)))
>     ############################################################################################
>
>
>     On Tue, Nov 21, 2023 at 10:53 AM Micha Silver <tsvibar using gmail.com>
>     wrote:
>
>         (Keeping on the list. And consider my suggestion to move this
>         to R-sig-Geo)
>
>
>         What is the coordinate reference system of the polygon layer?
>         Can you
>         share these layers?
>
>
>         On 21/11/2023 4:58, javad bayat wrote:
>         > Dear Micha;
>         > The R version is 4.3.2, the raster package is
>         "raster_3.6-26", the sf
>         > package is "sf_1.0-14".
>         > The str of the raster and polygon is as follow:
>         > > str(r)
>         > Formal class 'RasterLayer' [package "raster"] with 13 slots
>         >   ..@ file    :Formal class '.RasterFile' [package "raster"]
>         with 13 slots
>         >   .. .. ..@ name        : chr "E:\Base1.tif"
>         >   .. .. ..@ datanotation: chr "INT2S"
>         >   .. .. ..@ byteorder   : chr "little"
>         >   .. .. ..@ nodatavalue : num -Inf
>         >   .. .. ..@ NAchanged   : logi FALSE
>         >   .. .. ..@ nbands      : int 1
>         >   .. .. ..@ bandorder   : chr "BIL"
>         >   .. .. ..@ offset      : int 0
>         >   .. .. ..@ toptobottom : logi TRUE
>         >   .. .. ..@ blockrows   : Named int 128
>         >   .. .. .. ..- attr(*, "names")= chr "rows"
>         >   .. .. ..@ blockcols   : Named int 128
>         >   .. .. .. ..- attr(*, "names")= chr "cols"
>         >   .. .. ..@ driver      : chr "gdal"
>         >   .. .. ..@ open        : logi FALSE
>         >   ..@ data    :Formal class '.SingleLayerData' [package
>         "raster"] with
>         > 13 slots
>         >   .. .. ..@ values    : logi(0)
>         >   .. .. ..@ offset    : num 0
>         >   .. .. ..@ gain      : num 1
>         >   .. .. ..@ inmemory  : logi FALSE
>         >   .. .. ..@ fromdisk  : logi TRUE
>         >   .. .. ..@ isfactor  : logi FALSE
>         >   .. .. ..@ attributes: list()
>         >   .. .. ..@ haveminmax: logi TRUE
>         >   .. .. ..@ min       : num 1801
>         >   .. .. ..@ max       : num 3289
>         >   .. .. ..@ band      : int 1
>         >   .. .. ..@ unit      : chr ""
>         >   .. .. ..@ names     : chr "Base1"
>         >   ..@ legend  :Formal class '.RasterLegend' [package
>         "raster"] with 5
>         > slots
>         >   .. .. ..@ type      : chr(0)
>         >   .. .. ..@ values    : logi(0)
>         >   .. .. ..@ color     : logi(0)
>         >   .. .. ..@ names     : logi(0)
>         >   .. .. ..@ colortable: logi(0)
>         >   ..@ title   : chr(0)
>         >   ..@ extent  :Formal class 'Extent' [package "raster"] with
>         4 slots
>         >   .. .. ..@ xmin: num 330533
>         >   .. .. ..@ xmax: num 402745
>         >   .. .. ..@ ymin: num 3307321
>         >   .. .. ..@ ymax: num 3345646
>         >   ..@ rotated : logi FALSE
>         >   ..@ rotation:Formal class '.Rotation' [package "raster"]
>         with 2 slots
>         >   .. .. ..@ geotrans: num(0)
>         >   .. .. ..@ transfun:function ()
>         >   ..@ ncols   : int 5777
>         >   ..@ nrows   : int 3066
>         >   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
>         >   .. .. ..@ projargs: chr NA
>         >   ..@ srs     : chr "+proj=utm +zone=40 +datum=WGS84
>         +units=m +no_defs"
>         >   ..@ history : list()
>         >   ..@ z       : list()
>         >
>         > > str(p)
>         > sf [12 × 10] (S3: sf/tbl_df/tbl/data.frame)
>         >  \$ Name      : chr [1:12] "01" "02" "03" "04" ...
>         >  \$ FolderPath: chr [1:12]  ...
>         >  \$ SymbolID  : num [1:12] 0 0 0 0 0 0 0 0 0 0 ...
>         >  \$ AltMode   : int [1:12] 0 0 0 0 0 0 0 0 0 0 ...
>         >  \$ Base      : num [1:12] 0 0 0 0 0 0 0 0 0 0 ...
>         >  \$ Clamped   : int [1:12] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
>         >  \$ Extruded  : int [1:12] 0 0 0 0 0 0 0 0 0 0 ...
>         >  \$ Shape_Area: num [1:12] 63.51 15.28 9.23 39.19 642.38 ...
>         >  \$ Area      : num [1:12] 64 15 9 39 642 56 36 16 53 39 ...
>         >  \$ geometry  :sfc_POLYGON of length 12; first list element:
>         List of 1
>         >   ..\$ : num [1:18, 1:3] 55.6 55.6 55.6 55.6 55.6 ...
>         >   ..- attr(*, "class")= chr [1:3] "XYZ" "POLYGON" "sfg"
>         >  - attr(*, "sf_column")= chr "geometry"
>         >  - attr(*, "agr")= Factor w/ 3 levels
>         "constant","aggregate",..: NA NA
>         > NA NA NA NA NA NA NA
>         >   ..- attr(*, "names")= chr [1:9] "Name" "FolderPath"
>         "SymbolID"
>         > "AltMode" ...
>         > Sincerely
>         >
>         >
>         >
>         >
>         >
>         >
>         >
>         > On Mon, Nov 20, 2023 at 5:15 PM Micha Silver
>         <tsvibar using gmail.com> wrote:
>         >
>         >     Hi:
>         >
>         >     Some preliminary comments first
>         >
>         >     1- You might get a better response posting to R-sig-Geo
>         >
>         >     2- Probably better to switch to the `terra` package
>         >     older
>         >     `raster`
>         >
>         >     3- Furthermore, there is the package `exactextractr`
>         that ensures
>         >     that
>         >     edge pixels are taken into account, weighted by the
>         actual, partial
>         >     coverage of each pixel by the polygon.
>         >
>         >     4- Can you post: your R version, and version of the relevant
>         >     packages,
>         >     and information about the polygons and DEM (i.e.
>         `str(p)` and
>         >     `str(r)` )
>         >
>         >
>         >     On 20/11/2023 14:30, javad bayat wrote:
>         >     > Dear all;
>         >     > I am trying to calculate volume under each polygon of a
>         >     shapefile according
>         >     > to a DEM.
>         >     > when I run the code, it gives me an error as follows.
>         >     > "
>         >     > Error in h(simpleError(msg, call)) :
>         >     >    error in evaluating the argument 'x' in selecting a
>         method
>         >     for function
>         >     > 'addAttrToGeom': sp supports Z dimension only for
>         POINT and
>         >     MULTIPOINT.
>         >     > use `st_zm(...)` to coerce to XY dimensions
>         >     > "
>         >     > I want to have a table that contains one column
>         corresponding to the
>         >     > polygon rank and another that contains the volume on
>         that polyon.
>         >     > I would be more than happy if anyone could help me.
>         >     > I provided codes at the end of this email.
>         >     > Sincerely
>         >     >
>         >     >
>         >     >
>         >
>          ##########################################################################################
>         >     > library(raster);
>         >     > library(sf)
>         >     > # Load the DEM raster and shapefile
>         >     > r <- raster("E:/Base1.tif")
>         >     > p <- read_sf(dsn = "E:/Sites.shp", layer = " Sites")
>         >
>         >
>         >     Now to the question:
>         >
>         >     In `read_sf()`, when the source is a shapefile, you do
>         not need the
>         >     'layer' option. It doesn't hurt, but in this case I see
>         that you
>         >     have an
>         >     extra space before the layer name. Could that be causing
>         the problem?
>         >
>         >
>         >
>         >     > # Extract the values of the DEM raster for each polygon
>         >     > values <- extract(r, p)
>         >     >    Error in h(simpleError(msg, call)) :
>         >     >    error in evaluating the argument 'x' in selecting a
>         method
>         >     for function
>         >     > 'addAttrToGeom': sp supports Z dimension only for
>         POINT and
>         >     MULTIPOINT.
>         >     > use `st_zm(...)` to coerce to XY dimensions
>         >     >
>         >     > # Calculate the volume of each polygon
>         >     > volumes <- sapply(values, function(x) rasterVolume(x, r))
>         >
>         >     What is this "rasterVolume" function that you call here?
>         >
>         >
>         >     > # Print the results
>         >     > for (i in 1:length(volumes)) {
>         >     >    cat(sprintf("Volume under polygon %d: %f\n", i,
>         volumes[i]))
>         >     > }
>         >     >
>         >
>          ##########################################################################################
>         >     >
>         >     >
>         >     >
>         >     >
>         >     >
>         >     --
>         >     Micha Silver
>         >     Ben Gurion Univ.
>         >     Sde Boker, Remote Sensing Lab
>         >     cell: +972-523-665918
>         >
>         >
>         >
>         > --
>         > Best Regards
>         > M.Sc. Environment Engineering
>         > Alternative Mail: bayat194 using yahoo.com
>
>         --
>         Micha Silver
>         Ben Gurion Univ.
>         Sde Boker, Remote Sensing Lab
>         cell: +972-523-665918
>
>
>
>     --
>     Best Regards
>     M.Sc. Environment Engineering
>     Alternative Mail: bayat194 using yahoo.com
>
>
>
> --
> Best Regards