Adding band description to rioxarray to_raster()

Question

I've seen that one can add band descriptions to a geotiff image using rasterio [1]. How would I do the same thing when saving an array to a raster with rioxarray?

I tried adding the names as coords, but when I save an re-open the raster, the bands are named [1, 2, 3, 4] instead of ['R', 'G', 'B', 'NIR'].

import numpy as np
import xarray as xa
import rioxarray as rioxa

bands = ['R', 'G', 'B', 'NIR']
im_arr = np.random.randint(0, 255, size=(4, 400, 400))
im_save = xa.DataArray(im_arr, dims=('band', 'y', 'x'), 
        coords={'x': np.arange(0, 400), 'y': np.arange(0, 400), 
                'band': bands})
path = 'test.tiff'
im_save.rio.to_raster(path)
im_load = rioxa.open_rasterio(path)
print(im_load)

<xarray.DataArray (band: 4, y: 400, x: 400)> [640000 values with dtype=int32] Coordinates:

• band (band) int32 1 2 3 4
• y (y) float64 0.0 1.0 2.0 3.0 4.0 ... 396.0 397.0 398.0 399.0
• x (x) float64 0.0 1.0 2.0 3.0 4.0 ... 396.0 397.0 398.0 399.0 spatial_ref int32 0 Attributes: scale_factor: 1.0 add_offset: 0.0 grid_mapping: spatial_ref

Answer 1

You should consider switching from a 3d DataArray to a Dataset with 4 variables, each representing a separate band.

If you name the variables correctly, it should get written to the tiff:

import numpy as np
import xarray as xa
import rioxarray as rioxa

bands = ['R', 'G', 'B', 'NIR']

xa_dataset = xa.Dataset()

for band in bands:

    xa_dataset[band] = xa.DataArray(np.random.randint(0, 255, (400, 400), dtype="uint8"), dims=('y', 'x'), 
            coords={'x': np.arange(0, 400), 'y': np.arange(0, 400)})

# see the structure
print(xa_dataset)

# <xarray.Dataset>
# Dimensions:  (x: 400, y: 400)
# Coordinates:
#   * x        (x) int64 0 1 2 3 4 5 6 7 8 ... 391 392 393 394 395 396 397 398 399
#   * y        (y) int64 0 1 2 3 4 5 6 7 8 ... 391 392 393 394 395 396 397 398 399
# Data variables:
#     R        (y, x) uint8 18 41 126 79 64 215 105 ... 29 137 243 23 150 23 224
#     G        (y, x) uint8 1 18 90 195 45 8 150 68 ... 96 194 22 58 118 210 198
#     B        (y, x) uint8 125 90 165 226 153 253 212 ... 162 217 221 162 18 17
#     NIR      (y, x) uint8 161 195 149 168 40 182 146 ... 18 114 38 119 23 110 26


# write to disk
xa_dataset.rio.to_raster("test.tiff")

# load

im_load = rioxa.open_rasterio('test.tiff')
print(im_load)


# <xarray.DataArray (band: 4, y: 400, x: 400)>
# [640000 values with dtype=uint8]
# Coordinates:
#   * band         (band) int64 1 2 3 4
#   * y            (y) float64 0.0 1.0 2.0 3.0 4.0 ... 396.0 397.0 398.0 399.0
#   * x            (x) float64 0.0 1.0 2.0 3.0 4.0 ... 396.0 397.0 398.0 399.0
#     spatial_ref  int64 0
# Attributes:
#     scale_factor:  1.0
#     add_offset:    0.0
#     long_name:     ('R', 'G', 'B', 'NIR')
#     grid_mapping:  spatial_ref

You can see the band names are now included in the attributes as long_name.

Running gdalinfo, you can see the band description has been set:

Driver: GTiff/GeoTIFF
Files: test.tiff
Size is 400, 400
Origin = (-0.500000000000000,-0.500000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -0.5000000,  -0.5000000) 
Lower Left  (      -0.500,     399.500) 
Upper Right (     399.500,      -0.500) 
Lower Right (     399.500,     399.500) 
Center      (     199.500,     199.500) 
Band 1 Block=400x5 Type=Byte, ColorInterp=Red
  Description = R
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=400x5 Type=Byte, ColorInterp=Green
  Description = G
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=400x5 Type=Byte, ColorInterp=Blue
  Description = B
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=400x5 Type=Byte, ColorInterp=Alpha
  Description = NIR