Is there an option to maintain factor levels when save a `stars` object to a tif file?

Question

I have a raster in R as a stars object with an attribute that is a factor. When I save the raster to a tif file using write_stars, the factor levels are lost and only the integers are stored. I read somewhere that what I need is to create a Raster Attribute Table. Can this be done using stars?

Here is a reprex:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

m <- matrix(letters[1:9], nrow = 3, ncol = 3)
dim(m) <- c(x = 3, y = 3) # named dim
s <- st_as_stars(m)

s[["A1"]] <- factor(s[["A1"]])

plot(s)

write_stars(s, "reprex.tif")

^{Created on 2024-01-31 with reprex v2.1.0}

And this is how the legend of the tif looks in QGIS:

Answer 1

At the end it was not an R issue but a QGIS issue. The information is read by QGIS as it is stored in the xml:

<PAMDataset>
  <PAMRasterBand band="1">
    <CategoryNames>
      <Category></Category>
      <Category>a</Category>
      <Category>b</Category>
      <Category>c</Category>
      <Category>d</Category>
      <Category>e</Category>
      <Category>f</Category>
      <Category>g</Category>
      <Category>h</Category>
      <Category>i</Category>
    </CategoryNames>
  </PAMRasterBand>
</PAMDataset>

The only problem is that QGIS does not have support to visualize this in the legend, see this issue. Here a screenshot of this information in the properties of the raster in QGIS:

Currently, it seems that the only solution is using some plugins or some workarounds with python...

Answer 2

I'm assuming you have QGIS to test how this goes for your colleagues before you share with them, but I think I made things seem more complicated than needed in my comments. You want to look at ?terra::coltab for 'color table', that I take a detour through

s_rast = rast(s) 
plot(s_rast)
coltb = data.frame(values=1:9, col = rainbow(9, end = .9))
coltb
  values     col
1      1 #FF0000
2      2 #FFAC00
3      3 #A6FF00
4      4 #00FF06
5      5 #00FFB3
6      6 #009FFF
7      7 #0D00FF
8      8 #B900FF
9      9 #FF0099
has.colors(s_rast)
[1] FALSE
coltab(s_rast) <- coltb
plot(s_rast)
writeRaster(s_rast, 'sof77912406.tif') # this SO ? #
s_rast2 = rast('sof77912406.tif')
has.colors(s_rast2)
[1] TRUE
plot(s_rast2)
s_rast2_stars <- st_as_stars(s_rast2)
plot(s_rast2_stars)

s_rast2_stars
stars object with 2 dimensions and 1 attribute
attribute(s):
 sof77912406.tif 
 a      :1       
 b      :1       
 c      :1       
 d      :1       
 e      :1       
 f      :1       
 (Other):3       
dimension(s):
  from to offset delta         refsys x/y
x    1  3      0     1 WGS 84 (CRS84) [x]
y    1  3      3    -1 WGS 84 (CRS84) [y]

s_rast2
class       : SpatRaster 
dimensions  : 3, 3, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 3, 0, 3  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
source      : sof77912406.tif 
color table : 1 
categories  : A1 
name        : A1 
min value   :  a 
max value   :  i

write_stars(s_rast2_stars, 'sof77912406_stars.tif')
s_rast2_stars_2 = read_stars('sof77912406_stars.tif')
plot(s_rast2_stars_2)

This isn't satisfactory, yet, for your purposes as I've lost the factor values associated with the plot colors in the stars object... So, deeper digging in stars, but perhaps allows you to test colors as appear in QGIS.

Answer 3

You can do

r <- terra::rast(s)
writeRaster(r, "test.tif")