terra
terra copied to clipboard
rspatial
plotRGB: Consistent colors to compare several SpatRasters
I have a set of raster bricks that I want to display using plotRGB()
The problem is that plotRGB(x, stretch="lin") performs
ad-hoc stretching to each object, which is good in some cases
but not if you want to visually compare the plots.
Is there a way to enforce one given scale to be applied
in all plotRGB() displays?
Can argument zlim be used here?
library(terra)
b <- rast(system.file("ex/logo.tif", package="terra"))
plotRGB(b, zlim=c(100,255))
If zlim refers to the input object and the output range is always the same, yes. In that case, please clarify in the doc: zlim: "numeric vector of length 2. Range of values to plot (optional)" should be "Range of input values to plot (optional)" Is the max in zlim the same as scale? Also, this is an RGB plot, thus 3 bands are involved. Perhaps a 3x2 matrix would be better? Finally, for the user, it would be a lot simpler to just extract the scaling from one plotRGB() and apply it to another.
This is a little real example:
https://www.dropbox.com/s/p9tzrak01w3vowh/a.rds?dl=0
a <- readRDS("a.rda")
close.screen(all=TRUE)
split.screen(c(4,1))
#[1] 1 2 3 4
screen(1)
plotRGB(a)
#Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale) :
# color intensity 1.00587, not in [0,1]
screen(2)
plotRGB(a, zlim=c(0,4096))
#Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale) :
# color intensity 1.00587, not in [0,1]
screen(3)
plotRGB(a, stretch="lin")
screen(4)
plotRGB(a, zlim=c(0,4096), stretch="lin")
close.screen(all=TRUE)
I was not expecting the first 2 errors and was expecting screen 3 and 4 to be different (as the max of a bands are <<4096). I have images that do go up to 4096, and those should be displayed as much brighter.
The example does not work. You can save/read a SpatRaster with saveRDS/readRDS if you want, but not with save//load
Sorry about that. Here: https://www.dropbox.com/s/p9tzrak01w3vowh/a.rds?dl=0
> a <- rast(readRDS("a.rds"))
> global(a,range,na.rm=TRUE)
X1 X2
B170 950.0561 1865.897
B14 903.6075 1768.794
B86 1214.1776 2435.019
These are different issues. If your data are not within [0, 255] you need to provide argument scale. For example
plotRGB(a, scale=2500)
That should probably solve your problem.
You could also use the stretch method before using plotRGB.
I certainly use stretch() a lot and can achieve what I need by generating previous stretched 3-band versions of the images, but this is complicated for fast visualizaitions. My suggestion here is a faster method that would mimic what can be done in QGIS (and other RSGIS software): given 4 multi-spectral images a1, a2, a3, a4, you display a stretched color composite of a1 (or even just a stretched grey scale display of one given band), copy the style and paste it to the display of a2, a3 and a4 at once. I suggest being able to do so as easily in terra. For the grey-scale display, a solution is to create an intermediate object compiling the band of interest from each multi-spectral image and use tidyterra::geom_spatraster() with facet_wrap(), which ensures that the same scale is being applied. But again, we have to create an intermediate object.
And it could be really helpful if you made an example (with code, not files) of two SpatRasters that have different ranges, but for which you would like to colors to match.
All I suggest is something like: plotRGB(a1, r=7, g=5, b=2, stretch="lin") #assume this is a nice display estilo <- plotRGB(a1, r=7, g=5, b=2, stretch="lin")$style plotRGB(a2, style=estilo) plotRGB(a3, style=estilo) plotRGB(a4, style=estilo)
which can be done using stretch, but not as easily, and visualization is meant to be quick and easy. Then, for a final figure, a more complicated procedure can be followed.
When argument zlim is used with stretch=lim, the stretching now accounts for zlim. So you can simply set the limits you want/need to get consistent coloring, even with stretching different original value ranges.
library(terra)
a <- d <- b <- rast(system.file("ex/logo.tif", package="terra"))
d[d>250 | d < 100]=NA
a[a>150]=NA
par(mfrow=c(1,3))
plotRGB(b, stretch="lin", zlim=c(50,200))
plotRGB(d, stretch="lin", zlim=c(50,200))
plotRGB(a, stretch="lin", zlim=c(50,200))
The problem is that we do not necessarily know a priori that zlim=c(50,200) is the appropriate range. I know this can be found, but a procedure such as the one I suggest (extracting it from a useful visualization) would be much easier and faster, much better for data exploration, and equivalent to what is done in GIS packages. Anyway, I think my point is clear and you obviously have the right to close this if you prefer investing your time on other issues.
My approach is not because of time constraints. I do not want to add functionality now that I will remove later. I think that plot(RGB) should return an object that could recreate the plot. From that, you could possibly extract what you are after. This will take some work and time. So for now, something, like what I suggest would have to do.
@aloboa, for data visualization, I suggest using dedicated packages like "ggplot2" or "tmap" (or something else) -- they are much more advanced; you can save plots as objects and copy their attributes.
@t-wojciech I often use ggplot-based graphics, in particular, for this case, tidyterra geom_spatraster() and geom_spatraster_rgb(). The point here is a fast and simple visualization, which is the goal of terra::plotRGB(). My suggestion is making plotRGB() not only fast and simple, but also consistent across a set of raster objects in an equally simple manner. I really think my point is clear enough now and that this issue does not need further discussion.