addRasterRGB without stretching ?

[VicenteYago]

Hi, i have the following RGB image:

The visualization its the result of this code:

plotRGB(brick(r,g,b))

Note plotRGB its not performing stretching.

And if i run this code, with default stretching:

m <- leaflet() %>%
  setView(lng=IPEX12.latlong$longitud, lat=IPEX12.latlong$latitud, zoom = 20)  %>% 
  addProviderTiles(providers$Esri.WorldImagery) %>% 
  addRasterRGB(brick(b,g,r))
m

The result is:

my question is if addRasterRGB can be used without stretching so that the result is the same as plotRGB

I tried the following without succes:

  addRasterRGB(brick(b,g,r),quantiles = NULL, domain = NULL)

Thanks!.

[tim-salabim]

Hi @VicenteYago and sorry foe the late reply.

The issue here is a bit complicated and not immediately obvious.

plotRGB stretches values to be plotted between 0 and the maximum value of the supplied stack/brick. In the case below, it does use c(0, max_all) whereas addRasterRGB stretches between the minimum and the maximum of the supplied layers (as opposed to the full stack/brick), hence uses max_432. So, in order to mimic the behaviour of plotRGB you need to explicitly set the domain (and set quantiles = NULL).

In 72730138e6fd1ae6b0c37ce413d1ac2ff9fc7618 I have now changed the default for quantiles from c(0.02, 0.98) to c(0, 1) but that does not necessarily reproduce the output of plotRGB. I am reluctant to change the remaining default behaviour as I think that the shown values should be stretched between their respective min/max values rather than zero and the maximum of the whole stack/brick.

I'd love to hear your thoughts if you think plotRGB is doing the right thing and why. Also pinging @rhijmans to hear his thoughts on why plotRGB behaves the way it currently does. Maybe I am missing something obvious?

library(raster)
#> Loading required package: raster
#> Loading required package: sp
library(plainview, include.only = "poppendorf")
#> Loading required package: plainview
library(leaflet)
#> Loading required package: leaflet
library(leafem)

# supress all these proj related warnings for this example
options(warn = -1)

(max_all = max(plainview::poppendorf@data@max))
#> [1] 24972
(max_432 = max(plainview::poppendorf[[2:4]]@data@max))
#> [1] 15080
(min_432 = min(plainview::poppendorf[[2:4]]@data@min))
#> [1] 6154

# plotRGB default -> using max_all
plotRGB(plainview::poppendorf, 4, 3, 2)

# addRasterRGB default -> using max_432
leaflet() %>%
  addRasterRGB(
    plainview::poppendorf
    , 4, 3, 2
  )

# plotRGB -> using max_432
plotRGB(plainview::poppendorf, 4, 3, 2, scale =  max_432)

# addRasterRGB -> using domain = c(0, max_432)
leaflet() %>%
  addRasterRGB(
    plainview::poppendorf
    , 4, 3, 2
    , quantiles = NULL
    , domain = c(0, max_432)
  )

# plotRGB default -> using max_all
plotRGB(plainview::poppendorf, 4, 3, 2)

# addRasterRGB -> using max_all
leaflet() %>%
  addRasterRGB(
    plainview::poppendorf
    , 4, 3, 2
    , quantiles = NULL
    , domain = c(0, max_all)
  )

^{Created on 2021-05-24 by the reprex package (v2.0.0)}

[rhijmans]

The default in raster::plotRGB is to not stretch (in my terminology). But it needs to know the scale of the RGB values and assumes that these are between 0 and 255, but there is some logic that checks if there are higher values. If these are found, the maximum value found is used as maxColorValue in grDevices::rgb. However, it uses the max values of all the layers, not just the three RGB layers. In most cases that is probably OK, and sometimes it could be better (if all layers have the same scale). But I do not think that this was intentional, so I have changed it to use the max values for the layers only. The lowest values is not affected, because that is always assumed to be zero; it would be odd to have RGB values that do not start at zero.

If you know the scale, you could provide it

r <- plainview::poppendorf
plotRGB(r, 4, 3, 2, scale=2^14)

If you want stretching, you can do

plotRGB(r, 4, 3, 2, stretch="lin")

[tim-salabim]

Thanks @rhijmans ! Can you please clarify why

The lowest values is not affected, because that is always assumed to be zero;

?

[rhijmans]

The RGB channels (bands, color components) are intensities in red, green and blue, between 0 and n, where n is the maximum color depth (range, scale, resolution). n is often 255, or some other power of 2 (minus 1). zero is a natural origin; with c(0,0,0) == "black", no reflection, and c(n,n,n) == "white"; all light reflected. We just need to know, or guess, what n is. With remote sensing data n should be equivalent to the radiometric resolution.

If addRasterRGB subtracts the minimum observed value of the three channels, that could make the image too dark, as it would push the values towards zero (black).

Now if I look at plainview::poppendorf

#class      : RasterBrick 
#names      : B001n, B002n, B003n, B004n, B005n 
#min values :  8921,  7983,  6894,  6154,  5954 
#max values : 14123, 14060, 14390, 15080, 24972 

You could argue (and it was raster did) that `n`` should be 24972, or the nearest power of 2 that is larger than that (2^14)-1; even if that band is not used as a RGB channel. I believe that in this case (Landsat data) the numbers can go up to (2^16)-1 = 65535 (unsigned 2 byte integer). But even if that were theoretically correct, it cannot be used here, as the entire image becomes black. And the layers in a Raster* object could of course come from different sources, in which case all bets are off, and it may be safer to just use the values in the R, G, B, channels; also to always have the same result if these three layers are used, no matter what larger object they may be part of.

So in practice, stretching is used to avoid dull, or even entirely black, images. By default, plotRGB effectively does some stretching toward white (if the highest observed value > 255, and by then using it as n). By subtracting the minimum value, addRasterRGB effectively stretches towards black.

[tim-salabim]

Thanks again @rhijmans , I think I'm gonna stick with the current implementation to "stretch" the RGB values between min/max of the supplied layers. The user can decide to change it to e.g. c(0, max) or any other combination if they want to. But having a default that effectively shows the most possible detail is IMO desirable.