加权平均数计算与分区命令

Question

我希望对于每个类别，能够计算一个变量的加权平均值（因为每个像素都有一个权重）。下面的代码无法实现加权平均。如果您有任何建议将不胜感激！

#sample rasters
    r <- raster(ncols=10, nrows=10)
    category <-  raster(ncols=10, nrows=10)
    weight <-  raster(ncols=10, nrows=10)
    r[] <- runif(ncell(r)) * 1:ncell(r)
    category[] <- runif(ncell(category)) * 1:ncell(category)
    weight[] <- 1:ncell(weight)

#mean for each category
    zonalstats <- zonal(r, z=category, fun='mean', digits=0, na.rm=TRUE, count=T)

- Cecile

Answer 1

你可以使用这个函数。它将“zonal”应用于权重和栅格值，并计算加权平均值。

weightedZonalMean <- function(r, z, w, na.rm = TRUE, ...){
  zonalvalues <- zonal(r, z=z, fun=function(x, na.rm){x}, na.rm, ...)
  zonalweights <- zonal(w, z=z, function(x, na.rm){x},  na.rm, ...)
  weigtedValues <- cbind(zonalvalues[,1], 
                         sapply(1:nrow(zonalvalues), function(x){
    stats::weighted.mean(x = zonalvalues[x, 2:ncol(zonalvalues)], 
                         w = zonalweights[x, 2:ncol(zonalweights)])
  }))  
  return(weigtedValues)
}
weightedZonalMean(r, z, w, digits = 0)
#       [,1]      [,2]
#  [1,]    0 12.254438
#  [2,]    1  0.487548
#  [3,]    2  1.246249
#  [4,]    3 17.283858
#  [5,]    4 22.545906
#  [6,]    5 29.388179
#  [7,]    6 13.494853
#  [8,]    8  1.592464
#  [9,]    9 15.510155
# [10,]   10 10.091313
# ...

- loki

Answer 2

很遗憾，您不能使用区域函数实现您想要的功能，但是您可以构建一个类似的函数来计算加权平均值。

如果您运行findMethods(zonal)$'RasterLayer#RasterLayer'，您可以看到区域函数的工作方式，并通过调整它来创建一个解决您问题的函数。

zonalWMean <- function(x, w , z, digits = 0, na.rm = TRUE) {

  if(!compareRaster(c(x, w, z))) error("Raster not comparable!")

  df <- data.frame(x = getValues(x),
                   w = getValues(w))
  zones <- round(getValues(z), digits = digits)

  dfs <- split(df, zones)
  out <- lapply(dfs, function(l) weighted.mean(l$x, l$w, na.rm = na.rm))
  out <- cbind(as.numeric(names(out)), out)
  colnames(out)[1] <- "zone"
  colnames(out)[2] <- "weighted.mean"

  return(out)
}
zonalWMean(r, weight, category)
   zone weighted.mean
0  0    4.685347     
1  1    3.251086     
2  2    55.33882     
3  3    4.384738     
4  4    25.26902     
5  5    28.34853     
...

- Juan Antonio Roldán Díaz

---

非常感谢。在样本数据集上运行良好。但是，它不能与我使用的“r”层一起工作。每个区域都会出现NAs...我找不到问题所在。如果您有更多时间帮助我，我会分享我正在处理的3个图层（x、w、z）作为示例。谢谢！（EDU链接）https://cloudstor.aarnet.edu.au/plus/s/cz43y0zdxg5Qdrn - Cecile

你的 x_raster_layer.tif 文件中有 NA 值，请更新函数以使其能够处理 NA。 - Juan Antonio Roldán Díaz