USDAForestService / gdalraster

R Bindings to GDAL (Geospatial Data Abstraction Library)
https://usdaforestservice.github.io/gdalraster/
Other
41 stars 7 forks source link

plot_raster rgb scaling is not absolute #429

Closed mdsumner closed 2 months ago

mdsumner commented 3 months ago

I think each band is being scaled individually to its own range.

It should look like this (this plotted with terra::plotRGB and checked externally).

image

At this line, it needs to not stretch to 0,1 - or, do I have to specify that these data are 0,255 intended? I could not see how to do that.

https://github.com/USDAForestService/gdalraster/blob/c5ef54c37be515ed3921ef1b35687c7b5904aee8/R/display.R#L82

I guess it's arguable, but the bands are being stretched individually in this case afaict - I'm not sure what to do.

sq <- seq(2400, 3100, by = 100)
col <- col2rgb(hcl.colors(length(sq)))
writeLines(paste(sq, col[1, ], col[2, ], col[3, ]), pal <- tempfile(fileext = ".txt"))
library(gdalraster)
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
dem_proc("color-relief", lcp_file, tf <- tempfile(fileext = ".tif"), color_file = pal)

ds <- new(GDALRaster, tf)
plot_raster(ds)  ## gives the wrong scaling

Reading the values and using base graphics gives me the same correct result as terra, as above

v <- read_ds(ds)
plot(as.raster(array(v, attr(v, "gis")$dim)/255))

Scaling each band by stretching gives the wrong pic:

## what if we scale individual bands (each is stretched independently)
a <- array(v, attr(v, "gis")$dim)/255
for (i in 1:3) {
  a[,,i] <- scales::rescale(a[,,i])
}
plot(as.raster(a))

image

ctoney commented 3 months ago

specify that these data are 0,255 intended?

I believe that could be done with minmax_def = c(0, 0, 0, 255, 255, 255).

Maybe a better alternative for that specific case should use the existing maxColorValue argument? But it would need to be brought through to https://github.com/USDAForestService/gdalraster/blob/c5ef54c37be515ed3921ef1b35687c7b5904aee8/R/display.R#L110

when using the default color map function grDevices::rgb() (i.e., add the maxColorValue arg if col_map_fn is the default rgb()).

If we do that, it would become normalize = FALSE, maxColorValue = 255, which maybe is more intuitive once maxColorValue is documented accordingly.

Otherwise, I think normalizing by band would be preferred default in the general case. I'm guessing that is the difference noted here: https://github.com/hypertidy/ximage/issues/11#issuecomment-1942964350

But let me know if you disagree with that.

mdsumner commented 3 months ago

oh yes I remember we actually discussed this - it's much clearer to me now, the bands can be anything so you do really need that explicit metadata for the interpretation

I didn't get close enough to realize how minmax_def is structured, thanks!

I don't really get why per-band scaling is right though, for when you want an RGB? I always clip and stretch Sentinel equally, for example. What's an example where it's per band? (maybe I'm doing it wrong )

mdsumner commented 3 months ago

oh heck, I see your example over there is exactly that case. I can reproduce what ximage() does there with

 plot_raster(r, minmax_def = c(6000, 6000, 6000, 41396, 41396, 41396))

thanks for bringing me along in this. Hmm, I'm not sure what's sensible vs convenient now. I'll be paying more attention, mostly I'm consider actual byte-imagery prior to this kind of thinking (I think of the clip and stretch as data-processing, not vis - but I'll explore more).

ctoney commented 2 months ago

FWIW, the arguments related to min/max in plot_raster() were done to mimic the min/max settings in QGIS which are per-band:, so: plot_raster() default -> QGIS "Min/max" minmax_def -> "User defined" minmax_pct_cut -> "Cumulative count cut" (2 - 98 default)

and I did not implement QGIS "Mean +/- standard deviation x".

QGIS defaults to "Cumulative count cut" (2 - 98) if I add a multi-band image. This generally works well as a default, e.g., displaying a Landsat scene, typically there are outlier pixels and we don't want to normalize to the full range so 2 - 98 percentile tends to be a reasonable default. If I bring in a small clip out of a full scene, then the same setting may not be ideal. So in some gdalraster examples that display a small clip, I may use the 2 - 98 percentile ranges that were derived from the full scene, by manually specifying the ranges with minmax_def (e.g., https://usdaforestservice.github.io/gdalraster/articles/raster-display.html#rgb).

There is a of course a lot more to be considered, but keeping it simple, these settings at least provide flexibility to enable reasonable looking display for various band combinations / types of images.

All that said, I expect there is huge room for improvement.

ctoney commented 2 months ago

It should look like this (this plotted with terra::plotRGB and checked externally).

That is with no stretch to min/max applied (i.e., no contrast enhancement). plot_raster() gives that pic with

plot_raster(ds, nbands = 3, minmax_def = c(0, 0, 0, 255, 255, 255))

Scaling each band by stretching gives the wrong pic

It's not wrong, it just has contrast enhancement applied by default (by stretch to min/max). Personally I think it looks better with contrast enhancement. But either way these are 0-255 RGB. It's different if we're looking at band combinations in 16-bit imagery. In that case, not stretching to min/max tends to give a dull picture and we probably want contrast enhancement by default. So it's arguable what the defaults should be. But we have the options with minmax_def and minmax_pct_cut (plus normalize = TRUE|FALSE and col_map_fn) to render in different ways as needed (and matching a subset of the QGIS rendering options).

I may be missing something, let me know!

mdsumner commented 2 months ago

Oh it's just interesting to stretch by default, obviously I don't want stretch when I'm setting the colour map (and arguably is where I should bake the colours into a palette).

I'm not sure how to think about this yet, fundamentally different scenarios afaiu

ctoney commented 2 months ago

That's a good point. This case is a color map already applied. The default behavior should depend on the raster data type. A 3-band byte raster probably should not have stretch by default. But single-band probably should, and 16 bit or larger RGB basically always needs stretch for a decent picture. Turns out QGIS is doing exactly that. It is also configurable in settings but those are the defaults, so I lean toward matching those. Does that seem reasonable?

image

mdsumner commented 2 months ago

Absolutely, I just have to update my expectations and remember the context, thanks!

ctoney commented 2 months ago

The changes are described in comments at #435. I merged so that testing is easier and the online documentation can be checked, but I'm happy to revisit if you think it needs more attention. Thanks!

ctoney commented 2 months ago

@mdsumner, do you think there should be additional code changes or documentation for this, or good to close? The added documentation is mainly under Details. I agree with your original assessment that the previous behavior should be viewed as giving the wrong output. Although we could override the default normalization, it wasn't the default that users would probably expect for 8-bit RGB. It at least wasn't clearly documented that minmax_def may be needed in that case. Thanks again for the feedback.

mdsumner commented 2 months ago

yes, very good thanks!