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