[R] Calculating volume under polygons

[Micha Silver]

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