Writing datarray/dataset to netcdf results in bad values

[etsmith14]

I am trying to save a regional subset of a netcdf file as netcdf file. I am first opening some data with dimensions of time, latitude, and longitude and then slicing that data by latitude and longitude to produce a smaller subset of the data. I save that smaller subset with the to_netcdf command. But when I go to open the new netcdf, the timeseries definitely wrong (see figures). The figure named 'correct' is what the temperature timeseries looks like when plotting directly from the original dataset. The figure named 'wrong' is what the temperature timeseries looks like when plotting from the newly saved netcdf (hopefully both figures attached properly). This happens when I select just a single point and save the data as a netcdf and it also happens when I save as a zarr file. However, when I load a single netcdf with open_dataset (instead of open_mfdataset) and save it as a new netcdf, everything is correct. So the issue seems to be coming from open_mfdataset. I've also noticed that not all grid points are incorrect, only some grid points have this issue. This doesn't happen when I convert to a series then save as a CSV, just happens when saving as a netcdf or zarr.

import xarray as xr
import matplotlib.pyplot as plt

lats = [33.35]
lons = [-112.86]

# Load data from original files
ERA5_t2m = xr.open_mfdataset(r'E:\ERA5\Temperature\T2m_*' + '.nc', parallel=True).sel(latitude = slice(37.25,31), longitude = slice(-115, -109))
# plot original data (looks good)
plt.plot((((ERA5_t2m.t2m.sel(latitude = lats[0], longitude = lons[0], method='nearest')- 273.15) * 9/5) + 32))

# save regional subset as new netcdf
ERA5_t2m.t2m.to_netcdf(r"E:\Arizona_test\Arizona_Temperature.nc")
# open new netcdf regional subset
ERA5_t2m_AZ = xr.open_dataset(r'E:\Arizona_test\Arizona_Temperature.nc')
# plot same point from new netcdf (looks bad)
plt.plot((((ERA5_t2m_AZ.t2m.sel(latitude = lats[0], longitude = lons[0], method='nearest')- 273.15) * 9/5) + 32))

[rabernat]

That definitely seems weird. I have a couple of hunches about what might be going wrong, but they're hard to confirm with more detail.

Can you start by sharing more detailed metadata about the two datasets?

print(ERA5_t2m)
print(ERA5_t2m_AZ)

and also

print(ERA5_t2m.t2m.attrs)
print(ERA5_t2m_AZ.t2m.attrs)
print(ERA5_t2m.t2m.encoding)
print(ERA5_t2m_AZ.t2m.encoding)

[etsmith14]

Thanks for responding to my issue so quickly! Here is the print outs:

print(ERA5_t2m)
`<xarray.Dataset>
Dimensions: (longitude: 25, latitude: 26, time: 631152)
Coordinates:

• longitude (longitude) float32 -115.0 -114.8 -114.5 ... -109.5 -109.2 -109.0
• latitude (latitude) float32 37.25 37.0 36.75 36.5 ... 31.5 31.25 31.0
• time (time) datetime64[ns] 1950-01-01 ... 2021-12-31T23:00:00
  Data variables:
  t2m (time, latitude, longitude) float32 dask.array<chunksize=(744, 26, 25), meta=np.ndarray>
  mx2t (time, latitude, longitude) float32 dask.array<chunksize=(744, 26, 25), meta=np.ndarray>
  mn2t (time, latitude, longitude) float32 dask.array<chunksize=(744, 26, 25), meta=np.ndarray>
  Attributes:
  Conventions: CF-1.6
  history: 2022-05-03 17:18:58 GMT by grib_to_netcdf-2.24.3: /opt/ecmw...`

print(ERA5_t2m_AZ)
`<xarray.Dataset>
Dimensions: (longitude: 25, latitude: 26, time: 631152)
Coordinates:

• longitude (longitude) float32 -115.0 -114.8 -114.5 ... -109.5 -109.2 -109.0
• latitude (latitude) float32 37.25 37.0 36.75 36.5 ... 31.5 31.25 31.0
• time (time) datetime64[ns] 1950-01-01 ... 2021-12-31T23:00:00
  Data variables:
  t2m (time, latitude, longitude) float32 ...`

print(ERA5_t2m.t2m.attrs)
{'units': 'K', 'long_name': '2 metre temperature'}

print(ERA5_t2m_AZ.t2m.attrs)
{'units': 'K', 'long_name': '2 metre temperature'}

print(ERA5_t2m.t2m.encoding)
{'source': 'E:\\ERA5\\Temperature\\T2m_19501.nc', 'original_shape': (744, 109, 245), 'dtype': dtype('int16'), 'missing_value': -32767, '_FillValue': -32767, 'scale_factor': 0.0011997040993123216, 'add_offset': 269.40377331689564}

print(ERA5_t2m_AZ.t2m.encoding)
{'zlib': False, 'szip': False, 'zstd': False, 'bzip2': False, 'blosc': False, 'shuffle': False, 'complevel': 0, 'fletcher32': False, 'contiguous': True, 'chunksizes': None, 'source': 'E:\\Arizona_test\\Arizona_Temperature.nc', 'original_shape': (631152, 26, 25), 'dtype': dtype('int16'), 'missing_value': -32767, '_FillValue': -32767, 'scale_factor': 0.0011997040993123216, 'add_offset': 269.40377331689564}

[rabernat]

Ok thanks. I don't see anything obvious yet. Let's try one more.

print(ERA5_t2m.t2m.sel(latitude = lats[0], longitude = lons[0], method='nearest')
print(ERA5_t2m_AZ.t2m.sel(latitude = lats[0], longitude = lons[0], method='nearest')

[rabernat]

Ideally the MRE would not involve any files at all but only code! Something like this:

• Create a dummy Xarray dataset from scratch with the encoding that causes the issue.
• Save to netcdf
• Re-open, and show that it is different from the origin dataset

However, developing this is a lot of work, which you might not have time for. If you can reproduce the issue with a few small files, then you could just attach them directly to the github issue I believe.

Thanks for taking the time to help us improve xarray.

[etsmith14]

I created a single time series of temperature for Phoenix (1950 - 2021). Phx_Temperature.nc is the netcdf with no issues, Phx_Temperature_bad.nc is the netcdf with the issues. Hopefully, this data will suffice.

If the zip file doesn't work, I added the files to a box folder and made a publicly accessible link: https://epri.box.com/s/k0l5zfulcngim3ahluowiwlm33b5a6hd

Test.zip

[rabernat]

Thanks so much for taking the time to do that! 🙏

I realize it was not clear what I was asking of you. In order for this to be actionable by our developers, a bit more is required.

1. Go to https://github.com/pydata/xarray/issues
2. Open an issue following the template
3. Craft a minimal reproducible example of the bug using the attached files. Provide complete code that anyone can copy and paste to reproduce the bug. The goal is to make it as simple as possible for a potential contributor to explore and diagnose the issue. This will be a modified version of the code at the top of this post.

I know this is a lot of work on your part. You may decide you don't have time, which is also fine. But this is the process that will lead most directly to the bug getting fixed. This is how open source moves forward. 🙏

[etsmith14]

Okay, will do! Thanks for clarifying!

[veenstrajelmer]

Just for completeness: #7039