使用函数的 RGB 栅格中的归一化值答案

【问题标题】：Normalized values in RGB raster using a function使用函数的 RGB 栅格中的归一化值
【发布时间】：2021-12-30 13:58:02
【问题描述】：

我想将一个函数应用于 RGB 栅格，首先使用公式 DNstd = (DN-mean(DN)/sd(DN) 标准化原始值，然后是标准化过程，但现在使用公式 DNnor = 255*(DNstd-min(DNstd)/(max(DNstd)-min(DNstdmin)) 没有成功。我不想逐层应用，而是在单个函数中应用于所有栅格 (r)。在我的例子中：

 # Create a RBG raster
library(raster)
r1 <- raster(ncol=10, nrow=10) 
values(r1) <- rpois(100, 150)
r2 <- raster(ncol=10, nrow=10) 
values(r2) <- rpois(100, 200)
r3 <- raster(ncol=10, nrow=10) 
values(r3) <- rpois(100, 150)
r <- stack(r1, r2, r3)
names(r)=c("R","G","B")
r
# class      : RasterStack 
# dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
# resolution : 36, 18  (x, y)
# extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# names      :   R,   G,   B 
# min values : 124, 170, 126 
# max values : 183, 235, 181 

# Standardized the values using the first DNstd = (DN-DNmean)/DNstd 
#and second Normalized the values using now DNnor = 255*(DNstd-DNstdmin)/(DNstdmax-DNstdmin)
DNnor <- function(img, i, j, k) {
  bi <- img[[i]]
  bj <- img[[j]]
  bk <- img[[k]]
  bi_std <- (bi-mean(bi))/sd(bi) 
  bj_std <- (bj-mean(bj))/sd(bj)
  bk_std <- (bk-mean(bk))/sd(bk)
  bi_nor <- 255*(bi_std-min(bi_std))/(max(bi_std)-min(bi_std))
  bj_nor <- 255*(bj_std-min(bj_std))/(max(bj_std)-min(bj_std))
  bk_nor <- 255*(bk_std-min(bk_std))/(max(bk_std)-min(bk_std))
  return(bi_nor, bj_nor, bk_nor)
}
normalized_image <- DNnor(r, 1, 2, 3)
# Error in as.double(x) : 
# cannot coerce type 'S4' to vector of type 'double'

请问有什么帮助吗？

【问题讨论】：

标签： r raster r-raster

【解决方案1】：

请在下面找到显示一种可能解决方案的代表

两个初步说明：

第一个问题：一个函数只能返回一个对象。所以return(obj1, obj2, obj3) 是不可能的。因此，您需要在函数内部构建生成的stack 栅格，并在函数末尾构建return。
第二个问题：光栅堆栈是S4 object 类型。因此，无法使用img[[i]] 访问值。您必须使用以下语法img[[i]]@data@values。注意...之间的区别...

r[[1]]
#> class      : RasterLayer 
#> dimensions : 10, 10, 100  (nrow, ncol, ncell)
#> resolution : 36, 18  (x, y)
#> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : memory
#> names      : R 
#> values     : 113, 185  (min, max)

...和

r[[1]]@data@values
#>   [1] 128 150 153 151 148 143 150 147 155 147 140 140 148 153 162 152 142 156
#>  [19] 163 159 140 140 138 151 173 116 142 136 164 148 148 168 141 149 139 138
#>  [37] 164 144 160 135 147 145 132 152 149 123 124 167 177 144 165 141 156 185
#>  [55] 138 128 149 163 147 157 168 158 150 149 166 127 149 163 161 169 147 170
#>  [73] 160 137 160 129 166 140 135 162 165 145 150 155 152 129 139 136 143 156
#>  [91] 149 145 145 157 150 147 156 144 156 113

所以，请在下面找到我建议的解决方案

Reprex

函数代码

library(raster)

DNnor <- function(img, i, j, k) {
  bi <- img[[i]]@data@values
  bj <- img[[j]]@data@values
  bk <- img[[k]]@data@values
  bi_std <- (bi-mean(bi))/sd(bi) 
  bj_std <- (bj-mean(bj))/sd(bj)
  bk_std <- (bk-mean(bk))/sd(bk)
  bi_nor <- 255*(bi_std-min(bi_std))/(max(bi_std)-min(bi_std))
  bj_nor <- 255*(bj_std-min(bj_std))/(max(bj_std)-min(bj_std))
  bk_nor <- 255*(bk_std-min(bk_std))/(max(bk_std)-min(bk_std))
  
  r1n <- raster(ncol=10, nrow=10)
  values(r1n) <- bi_nor
  r2n <- raster(ncol=10, nrow=10)
  values(r2n) <- bj_nor
  r3n <- raster(ncol=10, nrow=10)
  values(r3n) <- bk_nor
  
  rn <- stack(r1n, r2n, r3n)
  
  return(rn)
}

函数的输出

normalized_image <- DNnor(r, 1, 2, 3)

normalized_image
#> class      : RasterStack 
#> dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
#> resolution : 36, 18  (x, y)
#> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> names      : layer.1, layer.2, layer.3 
#> min values :       0,       0,       0 
#> max values :     255,     255,     255

^{由reprex package (v2.0.1) 于 2021 年 12 月 30 日创建}

【讨论】：

你应该永远使用img[[i]]@data@values这样的语法。首先是因为那里可能没有值（它们通常在文件中）；第二，因为这可能不是内存安全的。如果您不担心内存限制，可以使用values(img[[1]]) 提取单元格值
非常感谢 @Robert Hijmans 的具体反馈，更一般地说，感谢您花时间分享您在 raster 和 terra 库上的知识。特别是，您的评论将使我能够纠正一种非常有问题的获取价值的方法！俗话说，一个人每天都在学习；-) 干杯。
不客气，当然，我们都在这里学习。在我看来，直接读取插槽 @ 的值要么是不好的使用，要么是不好的设计（也有例外）。写入插槽是一个真正的禁忌，因为设计可能是一个插槽的更改也会影响其他插槽。
感谢您的友好回复和指导。我确实会修改旧脚本以避免意外！我利用此评论向您询问更多信息。阅读您给出的答案，我顺便了解到 stretch() 函数根本不会将栅格加载到内存中，这与手动设计的函数在内部调用 values() 不同。我理解正确吗？如果是这样，如果一个人想要手动设计一个尽可能避免（或最小化）内存问题的函数，那么“策略”是什么？
raster 和 terra 都会在必要时以块的形式读取和处理大型数据集。在用于此问题的小示例中，没有区别。但通常安全总比后悔好。为了内存安全，请将您的自定义函数与软件包附带的函数一起使用，例如 terra 中的 _app 系列。手头的情况更复杂，因为 global 统计信息不能以块的形式计算。但是，您可以预先计算它们。

【解决方案2】：

您的示例数据

library(raster)
r <- raster(ncol=10, nrow=10) 
set.seed(1)
r1 <- setValues(r, rpois(100, 150))
r2 <- setValues(r, rpois(100, 200))
r3 <- setValues(r, rpois(100, 150))
r <- stack(r1, r2, r3)
names(r)=c("R","G","B")

有了这些示例数据，简单的解决方案是使用stretch

a <- stretch(r)

另见：

plotRGB(r, stretch="lin")

更一般地说，您的算法可能与使用stretch 后使用scale 相同

a <- scale(r) |> stretch()

如果您想编写自己的函数，您需要识别全局计算（所有单元格的单一统计，使用cellStats）和本地计算（每个单元格的新值）之间的区别。例如：

norfun <- function(x) {
    rmean <- cellStats(x, "mean")
    rsd <- cellStats(x, "sd")
    bstd <- (x-rmean)/rsd
    brng <- cellStats(bstd, "range")
    as.vector((255 / diff(brng))) * (bstd - brng[1,])
}

b <- norfun(r)

# same values
plot(a, b)

另请注意，cellStats(x, "mean") 是内存安全的等价于 mean(values(x))

您也可以使用terra 进行此操作

library(terra)
R <- rast(r)
A <- scale(R) |> stretch()

norfun2 <- function(x) {
    rmean <- unlist(global(x, "mean"))
    rsd <- unlist(global(x, "sd"))
    bstd <- (x-rmean)/rsd
    brng <- t(as.matrix(global(bstd, "range")))
    as.vector((255 / diff(brng))) * (bstd - brng[1,])
}
B <- norfun2(R)

plot(c(A,B), rast(stack(a, b)))

请注意，norfun2 不会单独处理每一层。如果您想要或需要这样做，您可以为单个层编写一个函数并将其应用于所有层。这是一个简化函数的示例，您不关心内存限制，因此您可以使用values() 创建单元格值向量。

norf <- function(x, ...) {
  v <- values(x)
  v <- (v - mean(v)) / sd(v)
  setValues(x, (255 / (max(v)-min(v))) * (v-min(v)))
}

# `sapp` loops over layers
z <- sapp(R, norf)

我愿意

(255 / (max(v)-min(v))) * (v-min(v)))

然后

255 * (v-min(v))) / (max(v)-min(v)))

预先计算常量，通过对整个向量进行更少的计算，这样可以节省一些时间。

【讨论】：