使用函数在 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'
请问有什么帮助吗?
请找到下面的代表,其中显示了一种可能的解决方案
两条初步说明:
第一个问题:一个函数只能return一个对象。所以 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 创建
您的示例数据
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)))
预先计算常量,通过对整个向量进行较少的计算来节省一些时间。
我想将一个函数应用于 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'
请问有什么帮助吗?
请找到下面的代表,其中显示了一种可能的解决方案
两条初步说明:
第一个问题:一个函数只能return一个对象。所以
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 创建您的示例数据
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)))
预先计算常量,通过对整个向量进行较少的计算来节省一些时间。