使用函数在 RGB 栅格中归一化值

Normalized values in RGB raster using a function

我想将一个函数应用于 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'

请问有什么帮助吗?

请找到下面的代表,其中显示了一种可能的解决方案

两条初步说明:

  1. 第一个问题:一个函数只能return一个对象。所以 return(obj1, obj2, obj3) 是不可能的。因此,您需要在函数内部构建生成的 stack 栅格,并在函数末尾构建 return 它。

  2. 第二个问题:栅格堆栈是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 创建

您的示例数据

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")

更一般地说,你的算法可能在使用scale

之后与使用stretch相同
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)))

预先计算常量,通过对整个向量进行较少的计算来节省一些时间。