1

I have a multiband raster of the first day temperatures exceed 15 degrees Celsius from 2000 to 2020 for a given region. As these are inferred from satellite images, we have cloud noise etc., therefore some pixel values of given bands will be way too low or too high (e.g. data for 20 bands are around day 100, but band 21 is at day 149 due to cloud interference - see picture).

Band01 120
Band02 10
Band03 13
Band04 117
Band05 149
Band06 124
Band07 3
Band08 129
Band09 134
Band10 145
Band11 109
Band12 130
Band13 133
Band14 135
Band15 133
Band16 117
Band17 122
Band18 128
Band19 125
Band20 127
Band21 149

I would like to set the multiband pixels experiencing this phenomenon to no-data and thereby not included them in the time-series. This would be done by calculating the overall mean and std.dev for the 21 bands and hereafter set the pixels to no-data if even one band is below or above the two standard deviations so something like:

if(multibands@value < mean-std.dev*2, "no-data", multibands@value) OR (multibands@value > mean+std.dev*2, "no-data", multibands@value)

I started something in R, but I thought I could also do it in QGIS.

library(raster)
stacked <- stack("path/temp15.tif")

meanR <- calc(stacked, fun = mean)
sdR <- calc(stacked, fun = sd)

upper <- meanR + sdR * 2
lower <- meanR - sdR * 2

result <- stack()

for i in nlayers(1:i){
   layer = stacked
   layer[layer > upper] = NA
   layer[layer < lower] = NA
   result = layer
}

writeRaster(result, path, overwrite = TRUE)

How can I set pixels to no-data if they exceed two standard deviations?

3
  • Your text doesn't match up with the image - Band 21 is 149 - not 166. And it would be better if you could have posted the "Image" as a table of data and that would have saved me some typing.
    – Spacedman
    Commented Jul 13, 2021 at 16:05
  • For that data the sd is 44, and the mean is 112, so no values are > mean+2sd (200), and three are less mean-2sd (24). So your output for this cell of this stack would be the same as the input except those three are set to NA? Yes?
    – Spacedman
    Commented Jul 13, 2021 at 16:07
  • Does this help? Yes Band 2, 3 and 7 would eventually be left out as NA. But how can I do that in my r script?
    – Thomas
    Commented Jul 14, 2021 at 14:06

1 Answer 1

0

This is the answer:

library(raster)
stacked <- stack("path/temp15.tif")

meanR <- calc(stacked, fun = mean)
sdR <- calc(stacked, fun = sd)

upper <- meanR + sdR * 2
lower <- meanR - sdR * 2

result <- stack()

for (i in nlayers(stacked)){
  layers = stacked
  layers[layers > upper] = NA
  layers[layers < lower] = NA
  result = layers
}

writeRaster(result, path, overwrite = TRUE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.