ds.rio.to_raster() not preserving band level attributes

mahrsee1997 commented 1 year ago

I have a xr.Dataset that contains around 43 bands/data variables. Each of those bands/data variables have 34 attributes [assume keys are same but values are different] + it have dataset level attributes as well.

Now I'm creating a tiff out of it using -- ds.rio.ro_raster('filename.tiff'). Then I opened the created tiff file using -- rioxarray.open_rasterio('filename.tiff', band_as_variable=True) I observed that all of the dataset level attributes are present but all of the band_level attributes are not preserved.

Lets us take this file. Note that gh variable had 29 attributes in dataset vs 3 attributes when converted into tiff.

Code:

import rioxarray
import cfgrib

ds = cfgrib.open_datasets(file_path)

to_convert_dataset = ds[1]

print(to_convert_dataset)
print(to_convert_dataset.gh)

to_convert_dataset.rio.to_raster('converted_tiff.tiff')

converted_tiff_ds = rioxarray.open_rasterio('converted_tiff.tiff', band_as_variable=True)

print(converted_tiff_ds)
print(converted_tiff_ds.band_1)

Console logs:

<xarray.Dataset>
Dimensions:     (latitude: 278, longitude: 370)
Coordinates:
    time        datetime64[ns] 2022-12-08
    step        timedelta64[ns] 00:00:00
    cloudBase   float64 0.0
  * latitude    (latitude) float64 0.138 0.246 0.354 0.462 ... 29.84 29.95 30.05
  * longitude   (longitude) float64 260.0 260.1 260.2 ... 299.6 299.7 299.9
    valid_time  datetime64[ns] ...
Data variables:
    gh          (latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
<xarray.DataArray 'gh' (latitude: 278, longitude: 370)>
[102860 values with dtype=float32]
Coordinates:
    time        datetime64[ns] 2022-12-08
    step        timedelta64[ns] 00:00:00
    cloudBase   float64 0.0
  * latitude    (latitude) float64 0.138 0.246 0.354 0.462 ... 29.84 29.95 30.05
  * longitude   (longitude) float64 260.0 260.1 260.2 ... 299.6 299.7 299.9
    valid_time  datetime64[ns] ...
Attributes: (12/29)
    GRIB_paramId:                             156
    GRIB_dataType:                            fc
    GRIB_numberOfPoints:                      102860
    GRIB_typeOfLevel:                         cloudBase
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_name:                                Geopotential height
    GRIB_shortName:                           gh
    GRIB_units:                               gpm
    long_name:                                Geopotential height
    units:                                    gpm
    standard_name:                            geopotential_height
<xarray.Dataset>
Dimensions:      (x: 370, y: 278)
Coordinates:
  * x            (x) float64 260.0 260.1 260.2 260.3 ... 299.5 299.6 299.7 299.9
  * y            (y) float64 0.138 0.246 0.354 0.462 ... 29.73 29.84 29.95 30.05
    spatial_ref  int64 0
Data variables:
    band_1       (y, x) float32 ...
Attributes:
    Conventions:             CF-1.7
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_edition:            2
    GRIB_subCentre:          0
    institution:             US National Weather Service - NCEP
<xarray.DataArray 'band_1' (y: 278, x: 370)>
[102860 values with dtype=float32]
Coordinates:
  * x            (x) float64 260.0 260.1 260.2 260.3 ... 299.5 299.6 299.7 299.9
  * y            (y) float64 0.138 0.246 0.354 0.462 ... 29.73 29.84 29.95 30.05
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     gh

Lib versions:

rasterio -- 1.3.2 rioxarray -- 0.13.1 xarray -- 2022.11.0 gdal -- 3.5.1

snowman2 commented 1 year ago

If you specify it like so, it works:

ds[1].gh.rio.to_raster('converted_tiff.tiff')
rds = rioxarray.open_rasterio("converted_tiff.tiff")
rds.attrs

{'GRIB_cfName': 'geopotential_height',
 'GRIB_cfVarName': 'gh',
 'GRIB_dataType': 'fc',
 'GRIB_gridDefinitionDescription': 'Latitude/longitude. Also called equidistant cylindrical, or Plate Carree',
 'GRIB_gridType': 'regular_ll',
 'GRIB_iDirectionIncrementInDegrees': 0.108,
 'GRIB_iScansNegatively': 0,
 'GRIB_jDirectionIncrementInDegrees': 0.108,
 'GRIB_jPointsAreConsecutive': 0,
 'GRIB_jScansPositively': 1,
 'GRIB_latitudeOfFirstGridPointInDegrees': 0.138,
 'GRIB_latitudeOfLastGridPointInDegrees': 30.054,
 'GRIB_longitudeOfFirstGridPointInDegrees': 260.0,
 'GRIB_longitudeOfLastGridPointInDegrees': 299.852,
 'GRIB_missingValue': 3.4028234663852886e+38,
 'GRIB_name': 'Geopotential height',
 'GRIB_numberOfPoints': 102860,
 'GRIB_NV': 0,
 'GRIB_Nx': 370,
 'GRIB_Ny': 278,
 'GRIB_paramId': 156,
 'GRIB_shortName': 'gh',
 'GRIB_stepType': 'instant',
 'GRIB_stepUnits': 1,
 'GRIB_typeOfLevel': 'cloudBase',
 'GRIB_units': 'gpm',
 'long_name': 'Geopotential height',
 'standard_name': 'geopotential_height',
 'units': 'gpm',
 'scale_factor': 1.0,
 'add_offset': 0.0}

snowman2 commented 1 year ago

Thanks for the report. Fix in #616

mahrsee1997 commented 1 year ago

@snowman2 ds[1].gh.rio.to_raster('converted_tiff.tiff') -- this is good when you want to create a tiff for a single band. But if you want to create a tiff of multiple bands (convert ds containing multiple data variables/bands into tiff). How will you do that -- like this ds[0].rio.to_raster('converted_tiff.tiff'), correct ?

Eg: each of these data variables pwat, tcc have different attributes + ds level attributes are also different.

print(ds[0])

<xarray.Dataset>
Dimensions:                (latitude: 278, longitude: 370)
Coordinates:
    time                   datetime64[ns] 2022-12-08
    step                   timedelta64[ns] 00:00:00
    atmosphereSingleLayer  float64 0.0
  * latitude               (latitude) float64 0.138 0.246 0.354 ... 29.95 30.05
  * longitude              (longitude) float64 260.0 260.1 260.2 ... 299.7 299.9
    valid_time             datetime64[ns] 2022-12-08
Data variables:
    unknown                (latitude, longitude) float32 ...
    pwat                   (latitude, longitude) float32 ...
    tcc                    (latitude, longitude) float32 ...
    tcolr                  (latitude, longitude) float32 ...
    tcols                  (latitude, longitude) float32 ...
    tcolw                  (latitude, longitude) float32 ...
    tcoli                  (latitude, longitude) float32 ...
    tcolc                  (latitude, longitude) float32 ...
    refc                   (latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP

I'm assuming that #616 will also be handling multiple bands case as well i.e. attributes will be preserved of each bands + ds level. Any timeline for merging & releasing #616 ?

corteva / rioxarray