I have a wide swath raster image. I want to detect cloudy pixels based on criteria such that a pixel in raster layer only consider its 25*25 neighborhood not the entire image.

My problem is here that I don't want to consider whole image when checking a pixel for cloudy.

I want a R code to do the following:

raster_layer = raster("path_to_raster_layer")
for (i in 1: nrow(raster_layer))
     for (j in 1:ncol(raster_layer))
        if(raster_layer(i,j)  < ( ( Q1 - 1.5 * ( Q3- Q1 )) of a neighborhood 
         of 25*25 )
         {  remove this pixel   }
        else { keep this pixel } 

Q1 = first quartile of neighborhood dataset Q3 = third quartile of neighborhood dataset


The raster package provides the function focal() (https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/focal). You first have to create a RasterLayer object by calling raster([path_to_your_raster_file]). Then define a custom cloud detection function which then can passed to focal(). Take a look at the examples in the documentation for further hints.

UPDATE 2018-08-17

This is probably what you need:

First step: calculate new raster with quantiles

raster_input <- raster("raster.tif")

fun <- function(x) {
  quantiles <- quantile(x, names = FALSE, na.rm = TRUE)
  return(quantiles[2] - 1.5*(quantiles[4] - quantiles[2]))

window <- matrix(data = 1, nrow = 21, ncol = 21)
raster_focal <- focal(x = raster_input, w = window, fun = fun)


I have defined a 21x21 matrix, because focal() requires a matrix with odd numbered shape as window.

You can also give the package spatial.tools a try. It provides a simple interface for parallel processing of raster images. This may speed up the calculations significantly (155 seconds vs. 257 seconds in the case of my test data). The Performance gain will probably increase with the size of your input raster.


raster_input <- raster("raster.tif")

fun <- function(inraster) {
  quantiles <- quantile(inraster, names = FALSE, na.rm = TRUE)
      return(quantiles[2] - 1.5*(quantiles[4] - quantiles[2]))

sfQuickInit(cpus=3) # this line is important for parallel processing
raster_focal_spatial_tools <- rasterEngine(
  inraster = raster_input,
  fun = fun,
  window_dims = c(21,21))


Second step: overlay

fun_over <- function(x, y) {
  return(ifelse(x < y, NA, x))

raster_over <- overlay(x = raster_input, y = raster_focal_spatial_tools, fun = fun_over)

  • I want a fiction that returns a value, in our case "quartile_1 - 1.5*(quartile_3-quartile_1)", which is calculated from the values of a 20*20 window around the central pixel. – Ritesh Porwal Aug 14 '18 at 10:48
  • Thank you so much for your time you spent on my problem. – Ritesh Porwal Oct 1 '18 at 5:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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