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.bayat...@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 <tsvi...@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
<tsvi...@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: bayat...@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: bayat...@yahoo.com
--
Best Regards
Javad Bayat
M.Sc. Environment Engineering
Alternative Mail: bayat...@yahoo.com
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
