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")
# Read the data
dem <- rast(dem_file)
polys <- st_read(poly_path)
# 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_crop <- mask(crop(dem, poly), poly)
??? 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 at
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
> p <- st_read("E:/Sites.shp")
> crop_r <- crop(r, extent(p))
> mask_r <- mask(crop_r, 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
> a = values(mask_r)
> # 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 at
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;
> > Thank you for your reply.
> > 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 at 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
> instead of the
> >? ? ?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
> > Javad Bayat
> > M.Sc. Environment Engineering
> > Alternative Mail: bayat194 at yahoo.com
>
> --
> Micha Silver
> Ben Gurion Univ.
> Sde Boker, Remote Sensing Lab
> cell: +972-523-665918
>
>
>
> --
> Best Regards
> Javad Bayat
> M.Sc. Environment Engineering
> Alternative Mail: bayat194 at yahoo.com
>
>
>
> --
> Best Regards
> Javad Bayat
> M.Sc. Environment Engineering
> Alternative Mail: bayat194 at yahoo.com
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918