R script to implement true color transformation for Sentinel-2 Quarterly Mosaic raw data, based on https://custom-scripts.sentinel-hub.com/sentinel-2/l2a_optimized/ and @[email protected]'s R implementation

sentinel_true_color.R

#    input: ["B04", "B03", "B02"],

DN <- 10000 # unit scaling. See page 40 in https://sentinel.esa.int/documents/247904/685211/Sentinel-2-MSI-L2A-Product-Format-Specifications.pdf

# make list of bands. any order, we'll use names
# adjust aggregation factor for quick testing
smp <- list(B02 = rast("B02.tif")  |> aggregate(2) / DN,
           B03 = rast("B03.tif")  |> aggregate(2) /DN,
           B04 = rast("B04.tif") |> aggregate(2) /DN
)

# Parameters for contrast and gamma
maxR <- 3.0 # max reflectance
midR <- 0.13
gamma <-1.8

evaluatePixel = function(smp) {
  ## R is B04, G is B03, B is B02
  rgbLin <-  satEnh(sAdj(smp$B04), sAdj(smp$B03), sAdj(smp$B02))
  return (list(sRGB(rgbLin[[1]]), sRGB(rgbLin[[2]]), sRGB(rgbLin[[3]])))
}

# pre-compute some parameters for gamma correction
gOff = 0.01;
gOffPow = gOff**gamma;
gOffRange = ((1 + gOff)**gamma) - gOffPow;

adjust_gamma <- function(b) {
  return ((((b + gOff)**gamma) - gOffPow) / gOffRange)
}

sAdj <-  function(a) {
  return (adjust_gamma(adj(a, midR, 1, maxR)))
}


#// Saturation enhancement 
sat <-  1.2

satEnh = function(r, g, b) {
  avgS = (r + g + b) / 3.0 * (1 - sat)
  return ( list(clamp(avgS + r * sat, lower = 0, upper = 1),
                clamp(avgS + g * sat, lower = 0, upper = 1), 
                clamp(avgS + b * sat, lower = 0, upper = 1)) )
}


# //contrast enhancement with highlight compression

adj <- function(a, tx, ty, maxC) {
  ar <- clamp(a / maxC, lower = 0, upper = 1)
  return (ar * (ar * (tx / maxC + ty - 1) - ty) / (ar * (2 * tx / maxC - 1) - tx / maxC))
}

##const sRGB = (c) => c <= 0.0031308 ? (12.92 * c) : (1.055 * Math.pow(c, 0.41666666666) - 0.055);

sRGB <- function(c) {
  thresh <-  0.0031308
  ifel(c <= thresh, #condition 1
       12.92*c, #TRUE
       (1.055 * c**0.4166666666) - 0.055 #true
       ) #false
}

esmp <-  evaluatePixel(smp)
rgbE <- rast(esmp[c(1,2,3)])

plotRGB(rgbE, scale = 1)

https://gist.github.com/vgXhc/f6fcffb37e4778b44e1379353edf78b0#file-sentinel_true_color-r-L61-L63

Is this nested ifel required? It seems to me that it can be removed because the code is applied directly to values that do not meet the first condition?

@kadyb You're right. I had adapted that from Barry's script. The nesting doesn't do anything. I'll edit.