I'd like to applied a function to a RGB raster, first standardized the original values using the formula DNstd = (DN-mean(DN)/sd(DN) and second the normalization process but now using the formula DNnor = 255*(DNstd-min(DNstd)/(max(DNstd)-min(DNstdmin)) without successes. I don't wanted to apply by layer, but to all the raster (r) in a single function. In my example:
# 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'
Please, any help with it?
CodePudding user response：
Your example data
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")
With these example data, the simple solution is to use stretch
a <- stretch(r)
Also see:
plotRGB(r, stretch="lin")
More generally, your algorithm may be the same as using stretch after using scale
a <- scale(r) |> stretch()
If you want to write your own function, you need to recognize the difference between global computation (a single statistic for all cells, using cellStats) and local computation (a new value for each cell). For example:
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)
Also note that cellStats(x, "mean") is the memory-safe equivalent to mean(values(x))
You can also do this with 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)))
Note that norfun2 does not treat each layer separately. If you wanted or needed to do that, you could write a function for a single layer and apply that to all. Here an example with a simplified function, where you are not concerned about memory limitations, so you can use values() to create a vector of cell 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)
I do
(255 / (max(v)-min(v))) * (v-min(v)))
Rather then
255 * (v-min(v))) / (max(v)-min(v)))
To pre-compute constants, and this safe some time by doing fewer computations on the entire vector.
CodePudding user response：
Please find the reprex below which shows one possible solution
Two preliminary remarks :
First problem: a function can only return one object. So
return(obj1, obj2, obj3)is not possible. So you need to build the resultingstackraster inside the function andreturnit at the end of the function.Second problem: a raster stack is a
S4 objecttype. So, accessing values withimg[[i]]is not possible. You have to use the following syntaximg[[i]]@data@values. Note the difference between...
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)
...and
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
So, please find below the solution I suggest
Reprex
- Code of your function
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)
}
- Output of the function
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
Created on 2021-12-30 by the reprex package (v2.0.1)