[R] Calculating volume under polygons

Tue Nov 21 08:23:44 CET 2023

(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 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 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 using yahoo.com

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918