将波段描述添加到 rioxarray to_raster()

如何解决将波段描述添加到 rioxarray to_raster()

我已经看到可以使用光栅 [1] 向 geotiff 图像添加波段描述。使用 rioxarray 将数组保存到栅格时，我将如何做同样的事情？

我尝试将名称添加为坐标，但是当我保存并重新打开光栅时，波段被命名为 [1,2,3,4] 而不是 ['R','G','B','近红外']。

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

bands = ['R','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,'band': bands})
path = 'test.tiff'
im_save.rio.to_raster(path)
im_load = rioxa.open_Rasterio(path)
print(im_load)

[dtype=int32 的 640000 个值] 坐标：

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 属性：比例因子：1.0 添加偏移量：0.0 网格映射：spatial_ref

解决方法

您应该考虑从 3d DataArray 切换到具有 4 个变量的 Dataset，每个变量代表一个单独的波段。

如果你正确命名变量，它应该被写入 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,'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','NIR')
#     grid_mapping:  spatial_ref

您可以看到波段名称现在作为 long_name 包含在属性中。

运行gdalinfo，可以看到band描述已经设置：

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