Stars: Does stars have an eqivalent way to create RGB rasters from arrays like raster::brick?

[MilesMcBain]

Here's my reprex:

img <- structure(c(143, 143, 128, 49, 8, 5, 25, 89, 143, 143, 143, 143, 
143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 137, 143, 
143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 133, 133, 
133, 133, 133, 136, 137, 143, 143, 143, 143, 143, 143, 143, 143, 
143, 143, 143, 154, 172, 172, 172, 154, 137, 133, 133, 136, 143, 
143, 143, 143, 143, 143, 143, 143, 143, 145, 145, 145, 145, 145, 
145, 125, 100, 79, 100, 137, 143, 143, 143, 143, 143, 143, 143, 
82, 82, 82, 82, 82, 82, 82, 82, 74, 88, 133, 137, 143, 143, 143, 
143, 143, 143, 232, 232, 232, 232, 232, 232, 232, 232, 232, 230, 
163, 133, 137, 143, 143, 143, 143, 143, 232, 232, 230, 230, 230, 
230, 230, 230, 230, 230, 230, 158, 133, 137, 143, 143, 143, 143, 
74, 158, 220, 230, 230, 230, 230, 230, 230, 230, 230, 230, 145, 
133, 143, 143, 143, 143, 133, 49, 33, 100, 172, 189, 203, 203, 
220, 230, 230, 230, 211, 133, 136, 143, 143, 143, 163, 128, 89, 
53, 18, 65, 137, 189, 230, 247, 238, 230, 230, 158, 133, 143, 
143, 143, 128, 125, 124, 158, 178, 107, 33, 33, 121, 220, 238, 
230, 230, 203, 133, 137, 143, 143, 125, 124, 133, 195, 195, 195, 
158, 44, 21, 143, 247, 232, 230, 230, 136, 136, 143, 143, 133, 
145, 136, 189, 125, 53, 29, 100, 203, 247, 247, 232, 230, 230, 
152, 133, 143, 143, 178, 163, 59, 25, 59, 133, 189, 195, 232, 
247, 247, 232, 230, 230, 163, 133, 143, 143, 65, 33, 74, 158, 
195, 195, 195, 195, 238, 247, 247, 230, 230, 230, 172, 133, 143, 
143, 152, 220, 230, 203, 172, 178, 195, 195, 247, 238, 230, 230, 
230, 232, 172, 133, 143, 143, 232, 230, 220, 172, 158, 178, 195, 
211, 247, 238, 230, 230, 230, 232, 154, 133, 143, 143, 189, 203, 
121, 107, 107, 107, 114, 178, 203, 211, 230, 230, 232, 232, 137, 
136, 143, 143, 8, 14, 2, 2, 2, 2, 5, 21, 21, 107, 232, 232, 232, 
203, 133, 137, 143, 143, 121, 133, 121, 121, 121, 124, 195, 220, 
220, 220, 232, 232, 232, 158, 133, 143, 143, 143, 124, 124, 133, 
133, 137, 203, 232, 232, 232, 232, 232, 232, 211, 136, 136, 143, 
143, 143, 145, 143, 154, 189, 230, 232, 232, 232, 232, 232, 232, 
230, 145, 133, 143, 143, 143, 143, 230, 230, 232, 232, 232, 232, 
232, 232, 232, 232, 232, 163, 133, 137, 143, 143, 143, 143, 230, 
230, 230, 230, 230, 230, 232, 232, 232, 230, 163, 133, 136, 143, 
143, 143, 143, 143, 230, 230, 230, 230, 230, 230, 230, 232, 211, 
145, 133, 136, 143, 143, 143, 143, 143, 143, 142, 142, 127, 49, 
8, 5, 25, 89, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 
142, 142, 142, 142, 142, 137, 142, 142, 142, 142, 142, 142, 142, 
142, 142, 142, 142, 142, 132, 132, 132, 132, 132, 136, 137, 142, 
142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 154, 171, 171, 
171, 154, 137, 132, 132, 136, 142, 142, 142, 142, 142, 142, 142, 
142, 142, 145, 145, 145, 145, 145, 145, 125, 100, 79, 100, 137, 
142, 142, 142, 142, 142, 142, 142, 82, 82, 82, 82, 82, 82, 82, 
82, 74, 88, 132, 137, 142, 142, 142, 142, 142, 142, 232, 232, 
232, 232, 232, 232, 232, 232, 232, 230, 163, 132, 137, 142, 142, 
142, 142, 142, 232, 232, 230, 230, 230, 230, 230, 230, 230, 230, 
230, 158, 132, 137, 142, 142, 142, 142, 74, 158, 220, 230, 230, 
230, 230, 230, 230, 230, 230, 230, 145, 132, 142, 142, 142, 142, 
132, 49, 33, 100, 171, 189, 203, 203, 220, 230, 230, 230, 211, 
132, 136, 142, 142, 142, 163, 127, 89, 53, 18, 65, 137, 189, 
230, 247, 238, 230, 230, 158, 132, 142, 142, 142, 127, 125, 123, 
158, 178, 107, 33, 33, 121, 220, 238, 230, 230, 203, 132, 137, 
142, 142, 125, 123, 132, 195, 195, 195, 158, 44, 21, 142, 247, 
232, 230, 230, 136, 136, 142, 142, 132, 145, 136, 189, 125, 53, 
29, 100, 203, 247, 247, 232, 230, 230, 151, 132, 142, 142, 178, 
163, 59, 25, 59, 132, 189, 195, 232, 247, 247, 232, 230, 230, 
163, 132, 142, 142, 65, 33, 74, 158, 195, 195, 195, 195, 238, 
247, 247, 230, 230, 230, 171, 132, 142, 142, 151, 220, 230, 203, 
171, 178, 195, 195, 247, 238, 230, 230, 230, 232, 171, 132, 142, 
142, 232, 230, 220, 171, 158, 178, 195, 211, 247, 238, 230, 230, 
230, 232, 154, 132, 142, 142, 189, 203, 121, 107, 107, 107, 114, 
178, 203, 211, 230, 230, 232, 232, 137, 136, 142, 142, 8, 14, 
2, 2, 2, 2, 5, 21, 21, 107, 232, 232, 232, 203, 132, 137, 142, 
142, 121, 132, 121, 121, 121, 123, 195, 220, 220, 220, 232, 232, 
232, 158, 132, 142, 142, 142, 123, 123, 132, 132, 137, 203, 232, 
232, 232, 232, 232, 232, 211, 136, 136, 142, 142, 142, 145, 142, 
154, 189, 230, 232, 232, 232, 232, 232, 232, 230, 145, 132, 142, 
142, 142, 142, 230, 230, 232, 232, 232, 232, 232, 232, 232, 232, 
232, 163, 132, 137, 142, 142, 142, 142, 230, 230, 230, 230, 230, 
230, 232, 232, 232, 230, 163, 132, 136, 142, 142, 142, 142, 142, 
230, 230, 230, 230, 230, 230, 230, 232, 211, 145, 132, 136, 142, 
142, 142, 142, 142, 142, 142, 142, 127, 49, 8, 5, 25, 89, 142, 
142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 
142, 137, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 142, 
142, 132, 132, 132, 132, 132, 135, 137, 142, 142, 142, 142, 142, 
142, 142, 142, 142, 142, 142, 153, 171, 171, 171, 153, 137, 132, 
132, 135, 142, 142, 142, 142, 142, 142, 142, 142, 142, 145, 145, 
145, 145, 145, 145, 125, 100, 79, 100, 137, 142, 142, 142, 142, 
142, 142, 142, 82, 82, 82, 82, 82, 82, 82, 82, 73, 87, 132, 137, 
142, 142, 142, 142, 142, 142, 232, 232, 232, 231, 231, 231, 231, 
231, 231, 230, 163, 132, 137, 142, 142, 142, 142, 142, 231, 231, 
230, 230, 230, 230, 230, 230, 230, 230, 230, 158, 132, 137, 142, 
142, 142, 142, 73, 158, 220, 230, 230, 230, 230, 230, 230, 230, 
230, 230, 145, 132, 142, 142, 142, 142, 132, 49, 33, 100, 171, 
189, 203, 203, 220, 230, 230, 230, 211, 132, 135, 142, 142, 142, 
163, 127, 89, 53, 18, 65, 137, 189, 230, 247, 238, 230, 230, 
158, 132, 142, 142, 142, 127, 125, 123, 158, 178, 107, 33, 33, 
121, 220, 238, 230, 230, 203, 132, 137, 142, 142, 125, 123, 132, 
195, 195, 195, 158, 44, 21, 142, 247, 231, 230, 230, 135, 135, 
142, 142, 132, 145, 135, 189, 125, 53, 29, 100, 203, 247, 247, 
232, 230, 230, 151, 132, 142, 142, 178, 163, 59, 25, 59, 132, 
189, 195, 232, 247, 247, 232, 230, 230, 163, 132, 142, 142, 65, 
33, 73, 158, 195, 195, 195, 195, 238, 247, 247, 230, 230, 230, 
171, 132, 142, 142, 151, 220, 230, 203, 171, 178, 195, 195, 247, 
238, 230, 230, 230, 231, 171, 132, 142, 142, 231, 230, 220, 171, 
158, 178, 195, 211, 247, 238, 230, 230, 230, 231, 153, 132, 142, 
142, 189, 203, 121, 107, 107, 107, 114, 178, 203, 211, 230, 230, 
231, 231, 137, 135, 142, 142, 8, 14, 2, 2, 2, 2, 5, 21, 21, 107, 
231, 231, 231, 203, 132, 137, 142, 142, 121, 132, 121, 121, 121, 
123, 195, 220, 220, 220, 231, 231, 231, 158, 132, 142, 142, 142, 
123, 123, 132, 132, 137, 203, 231, 231, 231, 231, 231, 232, 211, 
135, 135, 142, 142, 142, 145, 142, 153, 189, 230, 232, 232, 232, 
232, 232, 232, 230, 145, 132, 142, 142, 142, 142, 230, 230, 231, 
232, 232, 232, 232, 232, 232, 232, 231, 163, 132, 137, 142, 142, 
142, 142, 230, 230, 230, 230, 230, 230, 232, 232, 232, 230, 163, 
132, 135, 142, 142, 142, 142, 142, 230, 230, 230, 230, 230, 230, 
230, 232, 211, 145, 132, 135, 142, 142, 142, 142, 142, 142), .Dim = c(18L, 
                                                                      26L, 3L))

library(raster)
library(stars)
library(ggspatial)
library(mapview)
library(ggplot2)
library(palr)

## raster
r_tile <- raster::brick(img)
raster::crs(r_tile) <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"

ggplot() +
  layer_spatial(r_tile)
## correct

mapview::viewRGB(r_tile)
## yes

palr::image_raster(r_tile)
## works(?)

## stars

s_tile <- st_as_stars(img)
st_crs(s_tile) <- 3857

image(s_tile)
## fails

ggplot() +
  layer_spatial(s_tile)
## applying a colour scale I don't want

mapview::viewRGB(s_tile)
## fails

palr::image_stars(s_tile)
## fails

[MilesMcBain]

edit: added ggplot2 and palr to reprex.

[MilesMcBain]

edit: accidently wrote image_raster instead of image_stars for the last one. Same result anyhow.

[edzer]

I don't think that raster::brick creates an RGB raster, just try plot(r_tile) or image(r_tile). It interprets your array as a set (stack? brick?) of raster layers, and plots these. That ggspatial::layer_spatial() makes a different assumption is outside the scope of raster, and in addition it interpolates your cells, which is usually a no-go when showing raster data (where colors often correspond to categories).

Since your data seem not to be spatially referenced, I'd suggest to try your luck with one of the R packages dedicated to image analysis or image handling. Otherwise, plot(s_tile, rgb = 1:3) seems to do what you want (not stated anywhere, but this started here: https://twitter.com/MilesMcBain/status/1278293526766211072 ). I'll make a note that stars::geom_stars still lacks an rgb argument.

[MilesMcBain]

The data are spatially referenced. I know the extent and CRS, I just didn't include that for the reprex. Full context is here:https://github.com/anthonynorth/snapbox/blob/master/R/static_map.R#L27

Ultimately what I want to do is ggplot over the top of this thing. So right now I can't use geom_stars with 3 channel arrays of RGB values?

[paleolimbot]

Looks like my autodetection for "RGB" for stars images is not working...I think I based it off of what one gets when reading using read_stars() on a very specific image.

https://github.com/paleolimbot/ggspatial/blob/master/R/layer-spatial-stars.R#L23-L24

If you set the dimension names, it works but is incorrect...

s_tile <- st_as_stars(img) %>% 
  st_set_dimensions(names = c("x", "y", "band"))

It's not standalone, but you could rig your own Geom based on what I use to plot OSM tiles:

https://github.com/paleolimbot/ggspatial/blob/master/R/annotation-map-tile.R#L157-L199

[edzer]

But band is way too generic! I can imagine that we make this default if the (or a) dimension is called rgb and has cardinality 3.

[MilesMcBain]

Yes so if I write my image to PNG and read it back with read_stars, it's in a format that works nicely with ggspatial::layer_spatial.

stars object with 3 dimensions and 1 attribute
attribute(s):
    foo.png     
 Min.   :  2.0  
 1st Qu.:137.0  
 Median :143.0  
 Mean   :162.6  
 3rd Qu.:230.0  
 Max.   :247.0  
dimension(s):
     from to offset delta    refsys point values    
x       1 26      0     1 EPSG:3857    NA   NULL [x]
y       1 18     18    -1 EPSG:3857    NA   NULL [y]
band    1  3     NA    NA        NA    NA   NULL    

Is this an RGB raster?

[edzer]

No, it is a three-band raster dataset. The semantics that band 1, 2 and 3 refer to R, G and B values in 0-255 is lost.

[edzer]

How do you btw get read_stars to assign a crs to a PNG?

[MilesMcBain]

Oh sorry. I assigned the crs after reading from png file. :)

[edzer]

Instead of moving the rgb=... stuff to geom_stars, I propose to introduce a function that creates color rasters from multi-band rasters, so that we can handle there things like alpha, maxColorValue, and leave other dimensions untouched e.g.

library(stars)
s = st_as_stars(img)
r = st_apply(s, 1:2, function(x) rgb(x[1], x[2], x[3], maxColorValue=255))

has the colors in the proper places; supporting this by plot.stars and geom_stars is pretty straightforward.

[edzer]

st_rgb()?

[edzer]

This already works:

ggplot() + geom_stars(data=r)+ scale_fill_identity()

[edzer]

Here is an example how things look for a 4-dimensional object, using plot:

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
x1 = x[,,,1:3]
x2 = x[,,,4:6]
st_dimensions(x1) = st_dimensions(x2)
(xx = c(x1, x2, along = "foo"))
# stars object with 4 dimensions and 1 attribute
# attribute(s):
#   L7_ETMs.tif    
#  Min.   :  1.00  
#  1st Qu.: 54.00  
#  Median : 69.00  
#  Mean   : 68.91  
#  3rd Qu.: 86.00  
#  Max.   :255.00  
# dimension(s):
#      from  to  offset delta                       refsys point values    
# x       1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
# y       1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
# band    4   6      NA    NA                           NA    NA   NULL    
# foo     1   2      NA    NA                           NA    NA   NULL    
r = st_rgb(xx, 3)
plot(r)

[MilesMcBain]

This already works:

ggplot() + geom_stars(data=r)+ scale_fill_identity()

Sort of. This appears transposed from what I expect and get when using layer_spatial on a raster brick:

It's an M. For either Mapbox or openstreetMap.

[edzer]

fortunes::fortune("illogical")

Indeed, it matters how your matrix (or array) should be interpreted in terms of a (spatially oriented) image. You didn't tell.

[edzer]

@MilesMcBain : swapping dimensions 1 and 2 and then flipping the second (y) gives you

r = st_as_stars(aperm(img, c(2,1,3))[,18:1,])
ggplot()+geom_stars(data=st_rgb(r))+scale_fill_identity()

[edzer]

Feel free to re-open for follow-up issues!