(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
______________________________________________
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.
