-
Notifications
You must be signed in to change notification settings - Fork 105
Open
Description
GIven:
library(terra)
r <- rast(ncols=10, nrows=10)
values(r) <- 1:ncell(r)
z <- rast(r)
values(z) <- rep(c(1:2, NA, 3:4), each=20)
names(z) <- "zone"
The following does not work:
zonal(r, z, c("mean","sd"))
The following works, but the output is weird:
a <- zonal(r, z, function(x) {cbind(mean(x, na.rm=TRUE),sd(x, na.rm=TRUE))})
> str(a)
'data.frame': 4 obs. of 2 variables:
$ zone : int 1 2 3 4
$ lyr.1: num [1:4, 1:2] 10.5 30.5 70.5 90.5 5.92 ...
The output is quite impractical for a multi-layer r:
r <- c(r, r*2)
names(r) <- c("r1", "r2")
a <- zonal(r, z, function(x) {c(mean(x, na.rm=TRUE),sd(x, na.rm=TRUE))})
str(a)
> str(a)
'data.frame': 4 obs. of 3 variables:
$ zone : int 1 2 3 4
$ lyr.1.x: num [1:4, 1:2] 10.5 30.5 70.5 90.5 5.92 ...
$ lyr.1.y: num [1:4, 1:2] 21 61 141 181 11.8 ...
This is what I do:
df <- a$zone
for(i in 2:length(a)){
colnames(a[,i]) <- c(paste0(names(a)[i],".mean"), paste0(names(a)[i],".sd"))
df <- cbind(df, a[,i])
}
> df
df r1.mean r1.sd r2.mean r2.sd
[1,] 1 10.5 5.91608 21 11.83216
[2,] 2 30.5 5.91608 61 11.83216
[3,] 3 70.5 5.91608 141 11.83216
[4,] 4 90.5 5.91608 181 11.83216
But a behaviour such as the one of extract() would be a lot more convenient
(and/or perhaps a more straightforward solution exists)
Metadata
Metadata
Assignees
Labels
No labels