Using if-function to apply SCL depending pixel mutation

Hi everyone,

First of all, thank you for developing the openEO platform and the very usable R client — it is extremely helpful for our work!

We are new to using the R client and are currently trying to process a Sentinel scene for a later classification task. Specifically, we want to load a collection and set each pixel value in every band to NA (or something similar) when the SCL band does not have a value of 4, 5, or 6.
After that, we plan to reduce the time dimension by applying a median (excluding the NA values).

I have gone through the cookbook and process documentation and tried to come up with the correct code (using if or apply), but unfortunately, I am not able to make it work.

I’m not sure if this is the right forum for this question, but I would really appreciate any help or hints on how to achieve this.

Thank you very much!

This is my code so far in R Studio:

library(openeo)

# connect to the backend and authenticate
connection = connect(host = "https://openeo.dataspace.copernicus.eu")
login()

# assign the graph-building helper object to "p" for easy access to all openEO processes
p <- processes()

# create variables for loading collection
aoi <- list(west=16.06, south=48.06, east=16.65, north=48.35)
t <- c("2017-03-01", "2017-06-01")

# load first datacube
cube_s2 <- p$load_collection(
  id = "SENTINEL2_L2A",
  spatial_extent = aoi,
  temporal_extent = t,
  bands=c("B02", "B03", "B04","B05","B06","B07","B08","B8A","B09","B09","B11","B12","SCL")
)

# apply SCL based data manipulation
cube_s2 <- p$apply(cube_s2, function(x, context) {
  p$if(p$or(p$eq(data["SCL"], 4), p$eq(data["SCL"], 5), p$eq(data["SCL"], 6)), x ,NA)})

# reduce time dimension by setting median value
cube_s2_median_t <- p$reduce_dimension(data = cube_s2_mask,
                                       reducer = function(data, context) { p$median(data, na.rm = TRUE) },
                                       dimension = "t")

res <- p$save_result(data = cube_s2_median_t, format = formats$output$GTiff)

# send job to back-end
job <- create_job(graph = res, title = "cube_s2_median")

start_job(job = job)

jobs = list_jobs()
jobs

# download all the files into a folder on the file system
download_results(job = job, folder = "\somefolder")```

Hi,

two quick notes:

You seem to be using openeo.dataspace.copernicus.eu (Copernicus Data Space Ecosystem) which has its own forum at openEO - Copernicus Data Space Ecosystem | Community forum

Instead of the if process, I think it’s more elegant to use the mask process for this use case. I’m not that familiar with R, but in Python it would roughly go like this:

cube_s2 = connection.load_collection(...)
scl = cube_s2.band("SCL")
mask = (scl == 4) | (scl == 5) | (scl == 6)
masked = cube_s2.mask(mask)

Also note that the mask approach allows to do some morphological manipulation on the mask before applying it. For example you can achieve dilation/erosion by convolution with a gaussian kernel. There must be better examples out there, but this old outdated notebook illustrates that (see section “cloud masking”).

I hope this make sense and is translatable to R

Thanks Stefaan for your fast reply and your ideas/links - I will see if I can translate them it into R.

Update (just in case someone looks for a similar case):

Thanks also to the discussion on (cloud masking)

# load data cube with previously defined area of interest (aoi) and temporal extent (t)
cube_s2 <- p$load_collection(
  id = "SENTINEL2_L2A",
  spatial_extent = aoi,
  temporal_extent = t,
  bands=c("B02", "B03", "B04","B05","B06","B07","B08","B8A","B09","B11","B12","SCL")
)

# masking all bands according prefered SCL values to exclude e.g. clouded cells for later classification
cube_s2_mask <- p$apply(
  data = cube_s2,
  process = function(x, context) {
    scl = x["SCL"]
    scl_res = scl == 4 | scl == 5 | scl == 6
    p$`if`(value = scl_res, accept = x, reject = NULL)
  }
)

# reduce time dimension by setting median value and compute a median image for individual classification process
cube_s2_median_t <- p$reduce_dimension(data = cube_s2_mask,
                                       reducer = function(data, context) { p$median(data) },
                                       dimension = "t")

Hi,

I’d like to point out that there seems to be something conceptually wrong here:

The apply process in openEO is a process that provides it callback (the function provided in the process argument) individual “pixel” values of the cube. So the x argument should be assumed to be a number/float/int value, while you seem to be using it (x["SCL"]) as an array/mapping.
The easiest fix here would be to switch to apply_dimension along the “bands” dimension, which does provides the cube data as (labeled) array to the callback.