# Running a function on subset of a raster image instead of whole image in R

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:

``````library(raster)
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

## 1 Answer

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

``````library(raster)
raster_input <- raster("raster.tif")

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

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

plot(raster_focal)
``````

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.

``````library(spatial.tools)
library(raster)

raster_input <- raster("raster.tif")

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

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))
sfQuickStop()

plot(raster_focal_spatial_tools)
``````

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

plot(raster_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