Writing concatenated data to_netcdf file changes values

[mariomech]

Hej,

I try to combine a bunch of ERA5 datasets to a single one with open_mdfdataset() and write them to a single netcdf. While the combining seems to work, the method to_netcdf() seem sto change some of the values. Before writing the values in the dataset are still as in the origin file.

This is what I do:

def _remove_z(ds):
    if 'z' in ds.variables:
        ds = ds.drop('z')
    return ds
ds = xarray.open_mfdataset(['/tmp/first.nc','/tmp/second.nc'],preprocess=_remove_z,decode_times=False,combine='by_coords'); ds.to_netcdf('/tmp/result.nc')

The variable that changes is 't' the temperature. In the single file /tmp/second.nc I have

ds2 = xr.open_dataset('/tmp/second.nc')                                                                                                                                               
ds2['t'][:,16:17,:,:]      
                                                                                                                                                          
<xarray.DataArray 't' (time: 2, level: 1, latitude: 2, longitude: 2)>
array([[[[264.09262, 264.12222],
         [264.59088, 264.6254 ]]],
       [[[320.05222, 320.04236],
         [319.92395, 319.90915]]]], dtype=float32)
Coordinates:
  * longitude  (longitude) float32 11.85 11.95
  * latitude   (latitude) float32 78.95 78.85
  * level      (level) int32 17
  * time       (time) datetime64[ns] 2001-01-27T12:00:00 2001-01-28T12:00:00
Attributes:
    units:          K
    long_name:      Temperature
    standard_name:  air_temperature

But the result after ds.to_netcdf('/tmp/result.nc') is

rs = xr.open_dataset('/tmp/result.nc')                                                                                                                                               
rs['t'][:,16:17,:,:]                                                                                                                                                                  

<xarray.DataArray 't' (time: 3, level: 1, latitude: 2, longitude: 2)>
array([[[[234.8576 , 234.86235],
         [234.85287, 234.85524]]],
       [[[264.09366, 264.12207],
         [264.59097, 264.6265 ]]],
       [[[164.84744, 164.83795],
         [164.71954, 164.70297]]]], dtype=float32)
Coordinates:
  * latitude   (latitude) float32 78.95 78.85
  * longitude  (longitude) float32 11.85 11.95
  * level      (level) int32 17
  * time       (time) datetime64[ns] 2000-12-31T12:00:00 ... 2001-01-28T12:00:00
Attributes:
    units:          K
    long_name:      Temperature
    standard_name:  air_temperature

So

       [[[320.05222, 320.04236],
         [319.92395, 319.90915]]]

changes to

       [[[164.84744, 164.83795],
         [164.71954, 164.70297]]]

the moment I write to the file.

And for completness: python==3.8.5, xarray==0.15.0

Any help is very much appreciated.

Thank you.
first.nc.gz
second.nc.gz

[kmuehlbauer]

@mariomech I've reproduced this using xarray 0.16.2. The source of the problem is most likely that your source files have different scale_factor and add_offset for decoding the packed data (observe ds.t.encoding when opening the single datafiles).

Also your data isn't exported as packed data but as floating point data. I've no idea at which point the wrong decoding takes place but obviously that happens. It seems to work using `combine="nested" though.

That sounds like a serious issue.

[keewis]

the reason for this is encoding["dtype"]:

>>> ds.t.encoding["dtype"]
dtype('int16')
>>> rs.t.encoding["dtype"]
dtype('float32')

I'm not really sure why, though. Maybe someone with more knowledge on how our netcdf backend works can help with that?

You can work around this by adding del ds.t.encoding["dtype"] before calling to_netcdf:

In [2]: import xarray as xr
   ...: 
   ...: def _remove_z(ds):
   ...:     if 'z' in ds.variables:
   ...:         ds = ds.drop('z')
   ...:     return ds
   ...: 
   ...: 
   ...: ds = xr.open_mfdataset(
   ...:     ['first.nc', 'second.nc'],
   ...:     decode_times=False,
   ...:     preprocess=_remove_z,
   ...:     combine='by_coords',
   ...: )
   ...: del ds.t.encoding["dtype"]
   ...: ds.to_netcdf("result.nc")
   ...: 
   ...: with xr.open_dataset('result.nc', decode_times=False) as rs:
   ...:     xr.testing.assert_allclose(ds, rs, atol=1e-7)
   ...: 

In [3]:

[DaveOri]

@kmuehlbauer your idea is indeed correct. scale_factor and add_offset are different for the various datasets and xr.merge or xr.open_mfdataset keeps the encoding of only one of them.
When the values are encoded again as int16 some temperature loop over the int16 limit of 32767 and are written as negative integers.

The workaround suggested by @keewis works great, but it is still only a temporary patch for @mariomech that noticed this problem. Moreover, the drawback is that the variable will be written as float32 instead of int16.
The dataset can have a lot of variables, and this procedure should be repeated over all of them to avoid ambiguities.

It seems that this issue arises from a feature request
#1297
In that discussion, it was already pointed out that the encoding might be no longer valid. It was decided that it was not worth adding an additional argument that allows to drop the encoding as it was the case in previous versions of xarray.
My feeling is that the issue should be reopened. At least the problem should be part of the documentation, or?

[dcherian]

the values are encoded again as int16 some temperature loop over the int16 limit of 32767 and are written as negative integers.

Should we raise at compute time if the data cannot be represented by encoding["dtype"]? Is this easy to check?

[DaveOri]

As far as I can understand the open_mfdataset decodes everything into floats and keeps the encoding instructions of one (the first?) dataset.
The hard check would be looking if there are float values larger than the maximum or smaller than the minimum number that can be encoded given the type, the scale_factor and the add_offset.
Maybe it is easier to check if the datasets to merge have different encodings (perhaps check only if dtype is integer and only scale_factor and add_offset). In that case, at least, a warning should be raised so that the user knows that calling to_netcdf can cause data corruption.

[DanielAdriaansen]

Hi @dcherian @DaveOri. Thanks for this discussion. I ran into a similar issue following the documentation for combine_nested() here: https://docs.xarray.dev/en/stable/user-guide/combining.html?highlight=combine#combining-along-multiple-dimensions.

Each of my 16 subsets of the hypercube (each in their own netCDF) had their own add_offset and scale_factor and after using combine_nested() and then to_netcdf() the resulting data were incorrect because the scale_factor and add_offset for the full hypercube should be different than whichever of the 16 subsets it was using.

I was able to overcome my problem like:
ds.to_netcdf('out.nc',encoding={"my_var_1":{"dtype":"float32"}}

but I'd rather not have to do that for each my_var_n variable in my Dataset.

I was wondering if two things might be possible (in lieu of building full encoding support into to_netcdf() which perhaps is the long-term solution):

1. Set a flag or other metadata on the Dataset object to know whether or not it was created from other Datasets via any of the combining methods, and if so and to_netcdf() is called, default to float32 encoding and override the encoding stored when reading, and warn the user
2. Allow a global option for setting netcdf encoding, so my workaround above reduces to:
   ds.to_netcdf('out.nc',encoding={"global_var_dtype":"float32"}

[DanielAdriaansen]

I was able to use dictionary comprehension to fairly efficiently achieve item 2:

ds.to_netcdf('out.nc',encoding={v:{'dtype':'float32'} for v in ds.data_vars})

It might still be nice to see if there is a way to handle item 1 in my list.

I also discovered that if the input data have an encoding of byte or short, that the valid_min and valid_max NetCDF attributes are not re-computed to match the output encoding. Thus, I encountered a situation where the maximum value for my field was a floating value of 13482.13, but my valid_max was 127. Downstream use of the output NetCDF file was obstructed because the field max was > valid max. After stripping the valid_max and valid_min attributes before writing, I was able to workaround the downstream problems. This seems related to #1614.

[dcherian]

Hi @DanielAdriaansen

The long-term solution we're planning on is to remove encoding entirely when any operation akes place. Xarray does not explicitly update or track any "attributes", so there's always the possibility of encoding being out of date.

If you care about strict .encoding for the output file, specify it before calling to_netcdf