(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 at 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 at yahoo.com
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918