Might be worth an issue over at numpy with the example from the test.

numpy/numpy#23886

The issue is already resolved over at numpy which is really great! It was also marked as backport. @headtr1ck How are these issues resolved currently or how do we track removing the ignore?

Might be worth an issue over at numpy with the example from the test.

https://github.com/numpy/numpy/issues/23886

@vrishk Sorry for the delay here and thanks for bringing this to attention. We now have at least two requests which might move this forward (moving ensure_dtype_not_object into the backends). But this would need some discussion first, how to do this.

@mktippett Thanks for raising this. The issue should be cleared after #7888 is merged.

@tomwhite I've added tests to check the backend code for vlen string dtype metadadata. Also had to add specific check for the h5py vlen string metadata. I think we've covered everything for the proposed change to allow empty vlen strings dtype metadata.

I'm looking at the mypy error and do not have the slightest clue what and where to change. Any help appreciated.

Thanks @ghiggi for your comment.

The problem is we have at least two contradicting user requests here, see #7328 and #7862.

I'm sure there is a solution to accommodate both sides.

Thanks for trying. I can't think of any downsides for the netcdf4-fix, as it just adds the needed metadata to the object-dtype. But you never know, so it would be good to get another set of eyes on it.

So it looks like the changes here with the fix in my branch will get your issue resolved @tomwhite, right?

I'm a bit worried, that this might break other users workflows, if they depend on the current conversion to floating point for some reason. Also other backends might rely on this feature. Especially because this has been there since the early days when xarray was known as xray.

@dcherian What would be the way to go here?

There is also a somehow contradicting issue in #7868.

My main question here is, why is dask not trying to retrieve the object types from dtype.metadata? Or does it and fail for some reason?.

@tomwhite I've put a commit with changes to zarr/netcdf4-backends which should preserve the dtype metadata here: https://github.com/kmuehlbauer/xarray/tree/preserve-vlen-string-dtype.

I'm not really sure if that is the right location, but as it was already present that location at netcdf4-backend I think it will do.

@tomwhite Special casing on netcdf4 backend should be possible, too.

But it might need fixing at zarr backend, too:

python ds = xr.Dataset({"a": np.array([], dtype=xr.coding.strings.create_vlen_dtype(str))}) print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ds.to_zarr("a.zarr") print("\n### Loading ###") with xr.open_dataset("a.zarr", engine="zarr") as ds: print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ```python dtype: object metadata: {'element_type': <class 'str'>}

Loading

dtype: object metadata: None ```

Could you verify the above example, please? I'm relatively new to zarr :grimacing:

@ghiggi Glad it works, but we still have to check if that is the correct location for the fix, as it's not CF specific.

Thanks @tomwhite for the PR. I've only quickly checked the approach, which looks reasonable. But those changes have implications on several locations of the backend code, which we would have to sort out.

Considering this example:

```python import numpy as np import xarray as xr print(f"creating dataset with empty string array") print("-----------------------------------------") dtype = xr.coding.strings.create_vlen_dtype(str) ds = xr.Dataset({"a": np.array([], dtype=dtype)}) print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ds.to_netcdf("a.nc", engine="netcdf4")

print("\nncdump") print("-------") !ncdump a.nc

engines = ["netcdf4", "h5netcdf"] for engine in engines: with xr.open_dataset("a.nc", engine=engine) as ds: print(f"\nloading with {engine}") print("-------------------") print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ```

```python creating dataset with empty string array

dtype: object metadata: {'element_type': <class 'str'>}

ncdump

netcdf a { dimensions: a = UNLIMITED ; // (0 currently) variables: string a(a) ; data: }

loading with netcdf4

dtype: object metadata: None

loading with h5netcdf

dtype: object metadata: {'vlen': <class 'str'>} ```

Engine netcdf4 does not roundtrip here, losing the dtype metadata information. There is special casing for h5netcdf backend, though.

The source is actually located in open_store_variable of netcdf4 backend, when the underlying data is converted to Variable (which does some object dtype twiddling).

Unfortunately I do not have an immediate solution here.

@tomwhite Sorry for the delay here. I'll respond shortly on your PR #7862, but we might have to reiterate here later

@ghiggi Thanks for getting this back into action. I got dragged away from the one string object issue in #7654. I'll split this out and add a PR.

@ghiggi I'd appreciate if you could test your workflows against #7869. Your example and the one over in #7652 are working AFAICT.

The example below is only based on Variable and the cf encode/decode variable functions.

```python import xarray as xr import numpy as np

create DataArray

times = [np.datetime64("2000-01-01", "ns"), np.datetime64("NaT")] da = xr.DataArray(times, dims=["time"], name="foo") da.encoding["dtype"] = np.float64 da.encoding["_FillValue"] = 20.0

extract Variable

source_var = da.variable print("---------- source_var ------------------") print(source_var) print(source_var.encoding)

encode Variable

encoded_var = xr.conventions.encode_cf_variable(source_var) print("\n---------- encoded_var ------------------") print(encoded_var)

decode Variable

decoded_var = xr.conventions.decode_cf_variable("foo", encoded_var) print("\n---------- decoded_var ------------------") print(decoded_var.load()) ```

```python /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:618: RuntimeWarning: invalid value encountered in cast int_num = np.asarray(num, dtype=np.int64) /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:254: RuntimeWarning: invalid value encountered in cast flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype( /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:254: RuntimeWarning: invalid value encountered in cast flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(

---------- source_var ------------------ <xarray.Variable (time: 2)> array(['2000-01-01T00:00:00.000000000', 'NaT'], dtype='datetime64[ns]') {'dtype': <class 'numpy.float64'>, '_FillValue': 20.0} dtype num float64

---------- encoded_var ------------------ <xarray.Variable (time: 2)> array([ 0., 20.]) Attributes: units: days since 2000-01-01 00:00:00 calendar: proleptic_gregorian _FillValue: 20.0

---------- decoded_var ------------------ <xarray.Variable (time: 2)> array(['2000-01-01T00:00:00.000000000', 'NaT'], dtype='datetime64[ns]') {'_FillValue': 20.0, 'units': 'days since 2000-01-01 00:00:00', 'calendar': 'proleptic_gregorian', 'dtype': dtype('float64')} ```

@spencerkclark With current master I get the following RuntimeWarning running your code example:

• on encoding (calling to_netcdf()):

python /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:618: RuntimeWarning: invalid value encountered in cast int_num = np.asarray(num, dtype=np.int64)

• on decoding (calling open_dataset()):

python /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:254: RuntimeWarning: invalid value encountered in cast flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype( /home/kai/miniconda/envs/xarray_311/lib/python3.11/site-packages/xarray/coding/times.py:254: RuntimeWarning: invalid value encountered in cast flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(

The latter was discussed in #7098 (casting float64 to int64), the former was aimed to be resolved with this PR.

I'll try to create a test case using Variable and the respective encoding/decoding functions without involving IO (per your suggestion @spencerkclark).

Thanks @spencerkclark for taking the time. NaN has been written to disk (as you assumed). Let's have another try next week.

@maxhollmann I'm sorry, I'm still finding my way into Xarray. I've taken a closer look at #2377, especially https://github.com/pydata/xarray/issues/2377#issuecomment-415074188.

There @shoyer suggested to just use:

python data = duck_array_ops.where_method(data, ~mask, fill_value)

instead of python data[mask] = fill_value

I've checked and it works nicely with your test. That way we would get away without the flags test and the special handling will take place in duck_array_ops. Would be great if someone can double check.

combine_first uses fillna under the hood -> #3570

This is resolved in recent netcdf-c/netcdf4-python and works with recent Xarray.

@maxhollmann We might get at least some more views on this. There have been discussions on handling masked arrays and we should make sure this is exactly the solution we want to have.

@dcherian This changes as_compatible_data. Could you please have another look here? I'm a bit unclear about the implications.

@mx-moth Yes, this casting should be fixed.

I'm adding a bit of context here, as this might need to be solved in combination with #7098 and #7827. #7098 removes undefined casting for decoding. In #7827 there are efforts to do this for encoding, too.

As cast_to_int_if_safe is called for encoding as well as decoding I'm not sure if all cases have been catched by these two PR.

One issue on decoding is that at least for datetime64 based times the calculated time_deltas are currently converted to float64 in the presence of NaT (although NaT can perfectly be expressed as int64). It would be great if you could try your PR on top of #7827 (which includes #7098) to see if that fixes the errors in this PR.

@alimanfoo The slow code stems from my changes in #7400. Obviously the performance drop did not manifest in the tests/ benchmarks.

In #7824 @Illviljan is tackling concat performance.

@dcherian You were right from the beginning, changing order for decoding and handling _FillValue in CFDatetimeCoder seems to be one working solution with minimal code changes.

If the CI is happy I'll add tests to cover for the nanosecond issues in #7817.

Yet another idea would be to add and Engines heading on https://docs.xarray.dev/en/stable/ecosystem.html where engines/backends and there respective packages can be listed. The error could include a link to that page.

Maybe it would also help to rephrase the error, something along the lines

"Engine rasterio is not available. Please install the needed package. Engines [xxx, yyy, zzz] are available."

Thanks for the heads-up, @spencerkclark. No worries, I need to apply some changes anyway as it turns out.

All tests have passed. Rebased now on latest main. The issue described in #7817 is resolved. Ready for first reviews.

I've reset the order of coders to the initial behaviour. Instead the times are special cased in the CFMaskCoder. Locally it works, but I'll only trust the CI.

Setting status back to draft for now, still evaluating solutions for the CF encoding/decoding.

I've converted to draft for now, as I'm still evaluating solutions for the CF encoding/decoding.

@spencerkclark I'd appreciate if you could have a look here. All but one test pass, but I can't immediately see what that test is doing. Looks like mismatched dtypes on the attributes. If you have any suggestions how to possibly improve, please let me know. I've not added tests here, yet.

@dcherian Yes, I've setup a prototype in #7827. But the overall solution doesn't look that nice. The handling of fill_value has still to be done in CFMaskCoder.

Also #7098 is needed for this.

@gauteh No worries, glad it works now!

@paul0207 Thanks for providing the datafiles. I can't reproduce on my machine. Please provide more information, the output of xr.show_versions() would help and a complete traceback of the error you are experiencing.

A complete list of installed Python Packages would be nice (eg. by pip list), too.

Another couple of questions to get some more insight:

• Does this happen only with these special files, or do you experience this every time?
• Does the problem persists when specifying engine="netcdf4" or engine="h5netcdf" in the call to open_mfdataset?
• Does this also happen, if you open the files one-by-one (with xr.open_dataset) and combine the Datasets with xr.concat?

@gauteh Yes, please provide as much information as possible. It is also of interest, how you installed the package and what Python environment you are using (eg. system python, conda, venv etc.)

@gauteh You would probably have to delete this line:

https://github.com/gauteh/hidefix/blob/main/python/hidefix/xarray.py#L192

As @headtr1ck already explained, it is all handled via plugin system to be able to handle duplicate engine names on discovery by the python metadata.

cc @spencerkclark @DocOtak I've tried to at least find one example which incarnates as bug. Nevertheless the transformation from int to float in CFMaskCoder should be avoided.

We might think about special casing time data in CFMaskCoder, or handle masking of time data in CFDatetimeCoder/CFTimedeltaCoder.

@christine-e-smit Great this works on you side with the proposed patch in #7098.

Nevertheless, we've identified three more issues here in the debugging process which can now be handled one by one. So again, thanks for your contribution here.

As in #7098, citing @dcherian:

I think the real solution here is to explicitly handle NaNs during the decoding step. We do want these to be NaT in the output.

There are three more issues revealed here when using datetime64:

• if _FillValue is set in encoding, it has to be of same type/resolution as the times in the array
• If _FillValue is provided, we need to provide dtype and units to which fit our data, eg. if the _FillValue is referenced to unix-epoch the unit's should be equivalent
• when encoding in the presence of NaT the data array is converted to floating point with NaN, which is problematic for the subsequent conversion to int64

This is indeed an issue with scale_factor and add_offset as @d70-t has already mentioned.

That is not a problem per se, but those attributes are obviously different for different files. When concatenating only the first files's attributes survive. That might already be the source of the above problem, as it might slightly change values.

An even bigger problem is, when the dynamic range of the decoded data (min/max) doesn't overlap. Then the data might be folded from the lower border to the upper border or vica versa.

I've put an example into #5739. The suggestion for now is as @keewis comment to drop encoding in such cases and use floating point values for writing. You might use the available compression options for floating point data.

Xref: #5739

@christine-e-smit I've created an fresh environment with only xarray and zarr and it still works on my machine. I've then followed the Darwin idea and digged up #6191 (I've got those casting warnings from exactly the line you were referring to). Comment https://github.com/pydata/xarray/issues/6191#issuecomment-1209567966 should explain what happens here.

tl;dr citing @DocOtak

The short explanation is that the time conversion functions do an astype(np.int64) or equivalent cast on arrays that contain nans. This is undefined behavior and very soon, doing this will start to emit RuntimeWarnings.

There is also an open PR #7098.

Thanks @christine-e-smit for sticking with me to find the root-cause here by providing detailed information and code examples. :+1:

@christine-e-smit One more idea, you might delete the zarr folder before re-creating (if you are not doing that already). I've removed the complete folder before any new write (by putting eg. !rm -rf xarray_and_units.zarr at the beginning of the notebook-cell).

It would also be great if you could run the code from https://github.com/pydata/xarray/issues/7790#issuecomment-1529894939 and post the output here, just for the sake of comparison (please delete the zarr-folder before if it exists). Thanks!

@christine-e-smit I've plugged your code into a fresh notebook, here is my output:

```python

xarray created with NaT fill value

<xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02

xarray created read with NaT fill value

<xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02 {} {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -9223372036854775808, 'units': 'nanoseconds since 1970-01-01', 'calendar': 'proleptic_gregorian', 'dtype': dtype('int64')} ```

The output seems OK on my side. I've no idea why the data isn't correctly decoded as NaT on your side. I've checked that my environment is comparable to yours. The only difference remaining is you are on Darwin arm64 whereas I'm on Linux.

``` INSTALLED VERSIONS

commit: None python: 3.11.2 | packaged by conda-forge | (main, Mar 31 2023, 17:51:05) [GCC 11.3.0] python-bits: 64 OS: Linux OS-release: 5.4.0-144-generic machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: de_DE.UTF-8 LOCALE: ('de_DE', 'UTF-8') libhdf5: 1.14.0 libnetcdf: None

xarray: 2023.4.2 pandas: 2.0.1 numpy: 1.24.3 scipy: 1.10.1 netCDF4: None pydap: None h5netcdf: 1.1.0 h5py: 3.8.0 Nio: None zarr: 2.14.2 cftime: None nc_time_axis: None PseudoNetCDF: None iris: None bottleneck: None dask: 2023.3.2 distributed: 2023.3.2 matplotlib: None cartopy: None seaborn: None numbagg: None fsspec: 2023.3.0 cupy: None pint: None sparse: None flox: None numpy_groupies: None setuptools: 67.6.1 pip: 23.0.1 conda: None pytest: 7.2.2 mypy: 0.982 IPython: 8.12.0 sphinx: None ```

Unfortunately, I think you may have also gotten some wires crossed? You set the time fill value to 1900-01-01, but then use NaT in the actual array?

Yes, I use NaT because I want to check if the encoder does correctly translate NaT to the provided _FillValue on write.

So from your last example I'm assuming you would like to have the int64 representation of NaT as _FillValue, right? I'll try to adapt this, and see what I get

So, after some debugging I think I've found two issues here with the current code.

First, we need to give the fillvalue with a fitting resolution. Second, we have an issue with inferring the units from the data (if not given).

Here is some workaround code which (finally, :crossed_fingers:) should at least write and read correct data (added comments below):

```python

Create a numpy array of type np.datetime64 with one fill value and one date

FIRST ISSUE WITH _FillValue

we need to provide ns resolution here too, otherwise we get wrong fillvalues (day-reference)

time_fill_value = np.datetime64("1900-01-01 00:00:00.00000000", "ns") time = np.array([np.datetime64("NaT", "ns"), '2023-01-02 00:00:00.00000000'], dtype='M8[ns]')

Create a dataset with this one array

xr_time_array = xr.DataArray(data=time,dims=['time'],name='time') xr_ds = xr.Dataset(dict(time=xr_time_array))

print("******") print("Created with fill value 1900-01-01") print(xr_ds["time"])

Save the dataset to zarr

location_new_fill = "from_xarray_new_fill.zarr"

SECOND ISSUE with inferring units from data

We need to specify "dtype" and "units" which fit our data

Note: as we provide a _FillValue with a reference to unix-epoch

we need to provide a fitting units too

encoding = { "time":{"_FillValue":time_fill_value, "dtype":np.int64, "units":"nanoseconds since 1970-01-01"} } xr_ds.to_zarr(location_new_fill, mode="w", encoding=encoding)

xr_read = xr.open_zarr(location_new_fill) print("******") print("Read back out of the zarr store with xarray") print(xr_read["time"]) print(xr_read["time"].attrs) print(xr_read["time"].encoding)

z_new_fill = zarr.open('from_xarray_new_fill.zarr','r', ) print("******") print("Read back out of the zarr store with zarr")

print(z_new_fill["time"]) print(z_new_fill["time"].attrs) print(z_new_fill["time"][:]) ```

```python

Created with fill value 1900-01-01 <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02

Read back out of the zarr store with xarray <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02 {} {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -2208988800000000000, 'units': 'nanoseconds since 1970-01-01', 'calendar': 'proleptic_gregorian', 'dtype': dtype('int64')}

Read back out of the zarr store with zarr <zarr.core.Array '/time' (2,) int64 read-only> <zarr.attrs.Attributes object at 0x7f086ab8e710> [-2208988800000000000 1672617600000000000] ```

@christine-e-smit Please let me know, if the above workaround gives you correct results in your workflow. If so, then we can think about how to automatically align fillvalue-resolution with data-resolution and what needs to be done to correctly deduce the units.

```python xr_ds.to_zarr(location_new_fill,encoding=encoding)

xr_read = xr.open_zarr(location) print("******") print("Read back out of the zarr store with xarray") print(xr_read["time"]) print(xr_read["time"].encoding) ```

@christine-e-smit Is this just a remnant of copy&paste? The above code writes to location_new_fill, but reads from location.

Here is my code and output for comparison (using latest zarr/xarray):

```python

Create a numpy array of type np.datetime64 with one fill value and one date

time_fill_value = np.datetime64("1900-01-01") time = np.array([np.datetime64("NaT"), '2023-01-02'], dtype='M8[ns]')

Create a dataset with this one array

xr_time_array = xr.DataArray(data=time,dims=['time'],name='time') xr_ds = xr.Dataset(dict(time=xr_time_array))

print("******") print("Created with fill value 1900-01-01") print(xr_ds["time"])

Save the dataset to zarr

location_new_fill = "from_xarray_new_fill.zarr" encoding = { "time":{"_FillValue":time_fill_value,"dtype":np.int64} } xr_ds.to_zarr(location_new_fill, encoding=encoding)

xr_read = xr.open_zarr(location_new_fill) print("******") print("Read back out of the zarr store with xarray") print(xr_read["time"]) print(xr_read["time"].encoding) ```

```python

Created with fill value 1900-01-01 <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02

Read back out of the zarr store with xarray <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02 {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -25567, 'units': 'days since 2023-01-02 00:00:00', 'calendar': 'proleptic_gregorian', 'dtype': dtype('int64')} ```

This doesn't look correct either. At least the decoded _FillValue or the units are wrong. So -25567 is 1900-01-01 when referenced to of unix-epoch (Question: Is zarr time based on unix epoch?). When read back via zarr only this would decode into:

python <xarray.DataArray 'time' (time: 2)> array(['1953-01-02T00:00:00.000000000', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]')

I totally agree with @christine-e-smit, this is all very confusing. As said at the beginning, I have little knowledge of zarr. I'm currently digging into cf encoding/decoding which made me jump on here.

AFAICT, it looks like already the encoding has a problem, at least the data on disk is already not what we expect. It seems that somehow the xarray cf_encoding/decoding is not well aligned with the zarr writing/reading of datetimes.

@sbiner Sorry for the massive delay here. It doesn't have changed much since creation of your issue. Xarray doesn't take the netcdf default fill values into account (there are reasons, which @shoyer has explained in https://github.com/pydata/xarray/pull/5680#issuecomment-895455163 and https://github.com/pydata/xarray/pull/5680#issuecomment-895508489).

On write it just uses NaN as _FillValue (in case no specific encoding is given).

Xref: #2374, #7723, #5680

The code is there since #867 by @shoyer which was committed almost 7 years ago.

I've no idea what's the purpose for packing tuples into 0d arrays but as there are also tests for it in the above PR I'm assuming there is one real reason. Maybe @shoyer can chime in here to shed some light?

@wangshuaicumt Did you get along with this issue? If this is still unresolved it would be great if you could provide the data or a MCVE.

@AlxndrLhr I suppose your original issue is resolved. Please reopen or create a new issue if you still have problems with this.

It looks like this is no issue any more with recent versions of the stack. At least I can't reproduce this. @mjwillson Please reopen, if you still encounter problems while plotting.

@leicunxing-rs Sorry for the delay here. Your issue might be connected with concatenation/merge of several files containing packed data with different scale_factor/add_offset. See issue #5739 for more details (there they also merge different ERA5 datasets, hence the idea).

@dougrichardson Sorry for the delay. If you are still interested in the source of this issue here is what I found:

The root cause is different scale_factor and add_offset in the source files.

When merging only the .encoding of the first dataset survives. This leads to wrongly encoded file for the may-dates. But why is this so?

The issue is with the packed dtype ("int16") and the particular values of scale_factor/add_offset.

For feb the dynamic range is (228.96394336525748, 309.9690856933594) K whereas for may it is (205.7644192729947, 311.7797088623047) K.

Now we can clearly see that all our values which are above 309.969 K will be folded to the lower end (>229 K).

To circumvent that you have at least two options:

• change scale_factor and add_offset values in the variables .encoding before writing to appropriate values which cover your whole dynamic range
• drop scale_factor/add_offset (and other CF related attributes) from .encoding to write floating point values

It might be nice to have checks for that in the encoding steps, to prevent writing erroneous values. So this is not really a bug, but might be less impactful when encoding is dropped on operations (see discussion in #6323).

@floriankrb Sorry for the long delay. If you are still interested in the source of the issue, here is what I found:

By default Xarray will promote any data variable which shares it's name with a dimension to a coordinate. That accounts for ['number', 'time', 'step', 'heightAboveGround', 'latitude', 'longitude']. valid_time is a two dimensional coordinate (by CF standard) and is a coordinate here because t2m data variable has a corresponding coordinates-attribute containing valid_time. In the decoding-step valid_time gets added to the .coords. The attribute is removed from t2m's attrs and kept in t2m.encoding. So far so good.

By renaming number to n that coordinates attribute (in encoding) does not change as well. So when the data is written, t2m will still hold number in it's coordinates-attribute (on disk).

The issue manifests on subsequent read as now the decoding-step tries to align the found coordinates with the available data variables. As number is not available, no coordinate from that string will be taken into account as coordinate (note the all on line 444):

https://github.com/pydata/xarray/blob/0f4e99d036b0d6d76a3271e6191eacbc9922662f/xarray/conventions.py#L439-L447

This can easily be observed by looking into t2m.attrs where the coordinates remains instead of being preserved in .encoding.

So the source of all problems here is that the renaming number -> n was missed for coordinates-attribute of t2m's .encoding.

Can't reproduce with recent xarray/matplotlib/cartopy. Looks like this has been resolved.

python import xarray as xr import cartopy.crs as ccrs ds = xr.tutorial.load_dataset('air_temperature') ds = ds.sel(lon = slice(250, 300)) air = ds['air'] transform = ccrs.PlateCarree() projection = ccrs.Mercator(air.lon.values.mean(), air.lat.values.min(), air.lat.values.max()) p = air.isel(time=[0,1]).plot(transform = transform, aspect = ds.dims['lon']/ds.dims['lat'], col = 'time', col_wrap = 1, subplot_kws = {'projection': projection}) for ax in p.axs.flat: ax.set_extent((air.lon.values.min(), air.lon.values.max(), air.lat.values.min(), air.lat.values.max()), crs = transform) ax.set_aspect('equal', 'box')

Please reopen, if this is still an issue.

Thanks @dcherian for filling in the details.

I've digged up some more related issues: #2265, #3942, #4045

IIUC, #4684 did a great job to iron out much of these issues, but as it looks like only in the case when no NaT is within the time array (cc @spencerkclark). @christine-e-smit If you have no NaT in your time array then you can just omit encoding completely and Xarray will use int64 per default and your data should be fine on disk.

In the presence of NaT it looks like one workaround to circumvent that issue for the time being is to add the dtype in addition to _FillValue when writing out to zarr :

python encoding = { "time":{"_FillValue": time_fill_value, "dtype": np.int64} xr_ds.to_zarr(location, encoding=encoding) }

One note to this: Xarray is deducing the units from the current time data. So for the above example it will result in 'days since 2023-01-02 00:00:00' where days would now be the resolution in the file. If you want the resolution to be nanoseconds on disk units would need to be added to the encoding.

python encoding = { "time":{"_FillValue": time_fill_value, "dtype": np.int64, 'units': 'nanoseconds since 2023-01-02'} } xr_ds.to_zarr(location, encoding=encoding)

@christine-e-smit It would be great if you could confirm that from your side (some sanity check needed on my side).

@christine-e-smit I see, thanks for the details. AFAICT from the code it looks like zarr is special-cased in some ways compared to other backends. I'd really rely on some zarr-expert shedding light here and over at #7776.

@zoj613 Thanks for raising this.

The root-cause is that the tuple is returned from as_compatible_data as single element array:

python import xarray as xr print(xr.core.variable.as_compatible_data((2, 3, 4))) python array((2, 3, 4), dtype=object) This then breaks with the error you are seeing. I'm not quite sure if this is a bug in the code, a bug in the doc or no bug at all. But as a tuple is easily wrapped by np.array there should be a reason why Xarray is currently not able to digest tuples.

As we can see from the above output, in netCDF4-python scaling is adapting the dtype to unsigned, not masking. This is also reflected in the docs unidata.github.io/netcdf4-python/#Variable.

Do we know why this is so?

TL;DR: NETCDF3 detail to allow (signal) unsigned integer, still used in recent formats

• more discussion details on this over at https://github.com/Unidata/netcdf4-python/issues/656
• at NetCDF Users Guide on packed data:

A conventional way to indicate whether a byte, short, or int variable is meant to be interpreted as unsigned, even for the netCDF-3 classic model that has no external unsigned integer type, is by providing the special variable attribute _Unsigned with value "true". However, most existing data for which packed values are intended to be interpreted as unsigned are stored without this attribute, so readers must be aware of packing assumptions in this case. In the enhanced netCDF-4 data model, packed integers may be declared to be of the appropriate unsigned type.

My suggestion would be to nudge the user by issuing warnings and link to new to be added documentation on the topic. This could be in line with the cf-coding conformance checks which have been discussed yesterday in the dev-meeting.

Xref: discussion #7776, which got no attention up to now.

@christine-e-smit

So, I'm no expert for zarr, but it turns out that your NaT was converted to -9.223372036854776e+18 in the encoding step. I'm assuming that zarr is converting NaT as the format doesn't allow to use NaT directly, so it chooses a (default) value.

The _FillValue is not lost, but it will be preserved in the .encoding-dict of the underlying Variable:

python xr_read = xr.open_zarr(location) print("******************") print("No fill value") print(xr_read["time"]) print(xr_read["time"].encoding) ```python

No fill value <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02 {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -9.223372036854776e+18, 'units': 'days since 2023-01-02 00:00:00', 'calendar': 'proleptic_gregorian', 'dtype': dtype('float64')} ```

You might also check this without decoding (decode_cd=False):

python with xr.open_zarr(location, decode_cf=False) as xr_read: print("******************") print("No fill value") print(xr_read["time"]) print(xr_read["time"].encoding) ```python

No fill value <xarray.DataArray 'time' (time: 2)> array([-9.223372e+18, 0.000000e+00]) Coordinates: * time (time) float64 -9.223e+18 0.0 Attributes: calendar: proleptic_gregorian units: days since 2023-01-02 00:00:00 _FillValue: -9.223372036854776e+18 {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, 'dtype': dtype('float64')} ```

Maybe a zarr-expert can chime in here, what's the best practice for time-fill_values.

@maxhollmann I've checked and memory served well, the following issue might be related: #2377. It looks like your use-case is at least connected to @gerritholl's. It would be great if you could add your original use case (as MCVE, if possible) to get more details.

A special case (masked integer arrays) is discussed in #3955. As this might give additional information, it might not exactly fit your problem.

@maxhollmann I'll have a look into this, I think I've seen something like this some time ago.

Maybe you can add the tests to the PR or as comment? This might get more attention and will really help to debug.

I've marked this by accident, sorry @maxhollmann. Let us know when you feel this is ready

This is how netCDF4-python handles this data with different parameters:

python import netCDF4 as nc with nc.Dataset("http://dap.ceda.ac.uk/thredds/dodsC/neodc/esacci/snow/data/scfv/MODIS/v2.0/2010/01/20100101-ESACCI-L3C_SNOW-SCFV-MODIS_TERRA-fv2.0.nc") as ds_dap: v = ds_dap["scfv"] print(v) print("\n- default") print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- maskandscale False") ds_dap.set_auto_maskandscale(False) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- mask/scale False") ds_dap.set_auto_mask(False) ds_dap.set_auto_scale(False) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- mask True / scale False") ds_dap.set_auto_mask(True) ds_dap.set_auto_scale(False) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- mask False / scale True") ds_dap.set_auto_mask(False) ds_dap.set_auto_scale(True) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- mask True / scale True") ds_dap.set_auto_mask(True) ds_dap.set_auto_scale(True) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") print("\n- maskandscale True") ds_dap.set_auto_mask(False) ds_dap.set_auto_scale(False) ds_dap.set_auto_maskandscale(True) v = ds_dap["scfv"] print(f"variable dtype: {v.dtype}") print(f"first 2 elements: {v[0, 0, :2].dtype} {v[0, 0, :2]}") print(f"last 2 elements: {v[0, 0, -2:].dtype} {v[0, 0, -2:]}") ```python <class 'netCDF4._netCDF4.Variable'> int8 scfv(time, lat, lon) _Unsigned: true _FillValue: -1 standard_name: snow_area_fraction_viewable_from_above long_name: Snow Cover Fraction Viewable units: percent valid_range: [ 0 -2] actual_range: [ 0 100] flag_values: [-51 -50 -46 -41 -4 -3 -2] flag_meanings: Cloud Polar_Night_or_Night Water Permanent_Snow_and_Ice Classification_failed Input_Data_Error No_Satellite_Acquisition missing_value: -1 ancillary_variables: scfv_unc grid_mapping: spatial_ref _ChunkSizes: [ 1 1385 2770] unlimited dimensions: time current shape = (1, 18000, 36000) filling off

• default variable dtype: int8 first 2 elements: uint8 [215 215] last 2 elements: uint8 [215 215]

• maskandscale False variable dtype: int8 first 2 elements: int8 [-41 -41] last 2 elements: int8 [-41 -41]

• mask/scale False variable dtype: int8 first 2 elements: int8 [-41 -41] last 2 elements: int8 [-41 -41]

• mask True / scale False variable dtype: int8 first 2 elements: int8 [-- --] last 2 elements: int8 [-- --]

• mask False / scale True variable dtype: int8 first 2 elements: uint8 [215 215] last 2 elements: uint8 [215 215]

• mask True / scale True variable dtype: int8 first 2 elements: uint8 [215 215] last 2 elements: uint8 [215 215]

• maskandscale True variable dtype: int8 first 2 elements: uint8 [215 215] last 2 elements: uint8 [215 215] ```

First, the dataset was created with filling off (read more about that in the netcdf file format specs https://docs.unidata.ucar.edu/netcdf-c/current/file_format_specifications.html). This should not be a problem for the analysis, but it tells us that all data points should have been written to somehow.

As we can see from the above output, in netCDF4-python scaling is adapting the dtype to unsigned, not masking. This is also reflected in the docs https://unidata.github.io/netcdf4-python/#Variable.

If Xarray is trying to align with netCDF4-python it should separate mask and scale as netCDF4-python is doing. It does this already by using different coders but it doesn't separate it API-wise.

We would need a similar approach here for Xarray with additional kwargs scale and mask in addition to mask_and_scale. We cannot just move the UnsignedCoder out of mask_and_scale and apply it unconditionally.

@dcherian The main issue here is that we have two different CF things which are applied, Unsigned and _FillValue/missing_value.

For netcdf4-python the values would just be masked and the dtype would be preserved. For xarray it will be cast to float32 because of the _FillValue/missing_value.

I agree, moving the Unsigned Coder out of mask_and_scale should help in that particular case.

@dcherian Yes, that would work.

We would want to check the different attributes and apply the coders only as needed. That might need some refactoring. I'm already wrapping my head around this for several weeks now.

Then you are somewhat deadlocked. mask_and_scale=False will also deactivate the Unsigned decoding.

You might be able to achieve what want by using decode_cf=False (completely deactivate cf decoding). Then you would have to remove _FillValue attribute as well as missing_value attribute from the variables. Finally, you can run xr.decode_cf(ds) to correctly decode your data.

I'll add a code example tomorrow if no one beats me to it.

@Articoking

As both variables have a _FillValue attached xarray converts these values to NaN effectively casting to float32 in this case.

You might inspect the .encoding-property of the respective variables to get information of the source dtype.

You can deactivate the automatic conversion by adding kwarg mask_and_scale=False.

There is more information in the docs https://docs.xarray.dev/en/stable/user-guide/io.html

OK it seems this is ready for a first round of reviews.

A bit of added context:

Currently there is no dedicated function for checking for CF standard conformance. The idea is to read as much as possible also non-standard conforming data files, but restrict writing non-standard conforming files.

The implemented function ensure_scale_offset_conformance takes a strict keyword argument, which is True when encoding and False when decoding. If strict=True it will raise errors if there is a mismatch with the standard and when strict=False it will issue warnings.

I've only had to adapt a few tests which where not conforming to standard on encoding to align with that. I've observed some warnings in the test suite which we might to have a look into.

One idea would be to fix erroneous scale_factor/add_offset with our best fitting estimate. This is already done for list-type scale_factor/add_offset.

I will follow-up with checks for CFMaskCoder.

It's also possible to use the custom BackendEntrypoint-class directly in the call to xr.open_dataset with the engine keyword.

I think the equivalent incantation would be (note the different order of arguments in xr.where):

python da = xr.DataArray(np.arange(10)) print(xr.where(da < 5, da, 0).values) print(da.where(da < 5, 0).values) [0 1 2 3 4 0 0 0 0 0] [0 1 2 3 4 0 0 0 0 0]

@ChristmasZCY Please have a look at the documentation about string encoding

https://docs.xarray.dev/en/stable/user-guide/io.html#string-encoding

Good chance that this gives you the needed information.

Rebased on top main after merge of #7719. This is ready for review. It's a one-liner actually :grin:

@JoerivanEngelen Thanks for taking the time. Much appreciated.

This looks like it is ready to go. This will surely help further refactoring encode_cf_variable/decode_cf_variable. At least while working on it I spotted several locations where inconsistencies can be ironed out. A neat mostly flaw-free encoding/decoding is needed especially with regard to #6323.

I'm still convinced this could be fixed for floating point data.

Generally its worse if we obey some default fill values but not others, because it becomes quite confusing to a user.

I think this depends from which side you look at it :-) My point here is, we do not have to submissively obey to default fill values, but just use them when decoding. This only need to happen if no _FillValue is attached to the variable. By doing this we ensure that these missing values are mapped to np.nan (as it is expected by users).

In further course we can just apply the xarray standard np.nan when writing out. We need to document that in that case exact roundtrip isn't possible (it also isn't currently possible, in this example).

Consider this example:

```python dtype = "f4" with nc.Dataset("test-fillvalues-01.nc", mode="w") as ds: x = ds.createDimension("x", 10) test_fillval_fillon = ds.createVariable("test_fillval_fillon", dtype, ("x",), fill_value=nc.default_fillvals[dtype]) test_fillval_fillon[:5] = np.array([0.0, nc.default_fillvals[dtype], np.nan, 1.0, 8.0], dtype=dtype) test_nofillval_fillon = ds.createVariable("test_nofillval_fillon", dtype, ("x",), fill_value=None) test_nofillval_fillon[:5] = np.array([0.0, nc.default_fillvals[dtype], np.nan, 1.0, 8.0], dtype=dtype)

with nc.Dataset("test-fillvalues-01.nc") as ds: print("\n read with netCDF4-python") print("---------------------------") print(ds["test_fillval_fillon"]) print(ds["test_fillval_fillon"][:]) print(ds["test_nofillval_fillon"]) print(ds["test_nofillval_fillon"][:])

with xr.open_dataset("test-fillvalues-01.nc").load() as ds: print("\n read with xarray") print("---------------------------") print(ds["test_fillval_fillon"]) print(ds["test_fillval_fillon"][:]) print(ds["test_nofillval_fillon"]) print(ds["test_nofillval_fillon"][:]) python read with netCDF4-python

<class 'netCDF4._netCDF4.Variable'> float32 test_fillval_fillon(x) _FillValue: 9.96921e+36 unlimited dimensions: current shape = (10,) filling on [0.0 -- nan 1.0 8.0 -- -- -- -- --] <class 'netCDF4._netCDF4.Variable'> float32 test_nofillval_fillon(x) unlimited dimensions: current shape = (10,) filling on, default _FillValue of 9.969209968386869e+36 used [0.0 -- nan 1.0 8.0 -- -- -- -- --]

read with xarray-python

<xarray.DataArray 'test_fillval_fillon' (x: 10)> array([ 0., nan, nan, 1., 8., nan, nan, nan, nan, nan], dtype=float32) Dimensions without coordinates: x <xarray.DataArray 'test_nofillval_fillon' (x: 10)> array([0.00000e+00, 9.96921e+36, nan, 1.00000e+00, 8.00000e+00, 9.96921e+36, 9.96921e+36, 9.96921e+36, 9.96921e+36, 9.96921e+36], dtype=float32) Dimensions without coordinates: x ```

The only difference between these two variables is that on the first the _FillValue is declared, and on the other the default _FillValue is used. So if xarray obeys (by CF standard) the first it should also obey the second.

This might just work, if these cases the default fillvalue is used for decoding to np.nan, and declared that np.nan will be the new _FillValue. Does that make sense?

Now, this is interesting! It looks like those FillValue issues are following me. What did change that this now materializes here, all of a sudden.

Update: Small change - big issue. Checked for fv_exists instead of not fv_exists :grimacing:

The recommendation is to use _FillValue if there is only one value describing missing/fillvalue.

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#missing-data

It's also written that missing_value is

This attribute is not treated in any special way by the library or conforming generic applications, but is often useful documentation and may be used by specific applications.

https://docs.unidata.ucar.edu/netcdf-c/current/attribute_conventions.html

Not sure, if xarray is a conforming generic application or a specific application.

@dcherian Great, a duplicate. :-( Sorry I must have overlooked that one.

It's somewhat counter-intuitive to get differing results when using netcdf4-python and xarray. Would be a good idea to document this behaviour.

It looks like it might at least be resolved for floating point source data:

Let's take the above simple example. We have np.nan written to the file, but the netcdf representation on disk uses a default (undeclared by attribute) _FillValue for unwritten parts.

For the netcdf4-python user the np.nan will not be masked, but the unfilled parts will be masked.

For xarray the default fillvalue won't be masked, appearing as valid data, which it is not. On subsequent writes np.nan will be introduced as the new fillvalue (by attribute), effectively changing the meaning of the default fillvalues.

Wouldn't it make sense then, to transform these default fill values to np.nan on read too, instead of giving the a seemingly meaningful value? Maybe yet another keyword switch, use_default_fillvalues?

There should be at least a warning on read, in these situations, that there are undefined values in the dataset which were never written and which will not be masked.

If the dataset contains unwritten parts, and a default fillvalue is used, in turn meaning the data creator did this by purpose (by not setting a _FillValue) it can mean several things:

• The creators data does actually not have missing values which need declaring, but it means, that his data will get masked for default fillvalue entries (maybe they doesn't know about this, but that might be unlikely).
• The creator doesn't care at all, with same conclusion as above.
• The creator purposefully uses default fillvalue as missing value, since they use this as a means of saving disk space. But this could also be done, by just defining that as _FillValue attribute at creation time, if I`m not mistaken.

I'm still convinced this could be fixed for floating point data.

Please check #7720 if that fixes the conversion problems. Thanks.

If it helps to minimize interoperability issues I'm all in for the change. One thing I would maybe do is wait for the next version. With the current PR we would end up with two different build numbers with differing behaviour, which might confuse folks.

But I'd rely on @ocefpaf's expertise.

As explained I've created two PR (#7719 and #7720) for the "easy" changes from this PR. Would be great, if those could go in fast. Thanks!

@dcherian Just a heads-up: I find this PR getting more and more involved at different parts of the machinery and hard to follow for reviewers. I'll split this up and start with the more or less undisputed changes.

Still hunting for corner cases and issues inside encode_cf_variable/decode_cf_variable.

It looks like I already see some light again. Not sure, if this is the last iteration, but the testsuite is still running green with added and enhanced tests, which is not that bad.

Unfortunately https://github.com/pydata/xarray/issues/2304 is still an issue for now. I'll clarify that later with an added test.

While trying to create a test which specifically tests _choose_float_dtype I've found some issues with checking for availability of scale_factor/add_offset. Now testing for None.

This is now ready for another round of reviews, @dcherian, @Illviljan and @mankoff.

As @mankoff already pointed out, xarray is very generous to try to encode/decode non CF conforming data. This makes things a bit complicated as some issues only surface in rare corner cases.

I've tried to be as explicit in _choose_float_dtype, also added comments/tests where needed.

I'm finding the typing a bit hard. It seems that mypy can't derive the correct types from return types in certain cases.

If at first you don't succeed... It looks like we have something working here.

Some more typing and maybe some more tests covering the cases with scale_factor/add_offset/_FillValue non-conforming CF and we should be good to go. Or do I miss something?

The latest changes brake #1840 again. We have two contradicting forces here, which need to be aligned.

To fix this, I think logic in _choose_float_dtype should be updated to look at encoding['dtype'] (if available) instead of dtype, in order to understand how the data was originally stored.

This is aimed at in #7654

@Illviljan I'm not able to figure out the typing if I want to use Data-types as functions to convert python numbers to array scalars. If you have any suggestion how to fix this, please let me know.

@dcherian @Illviljan Thanks for the first round of review. I've rebased everything on latest main. Now the code moving from conventions.py to coding.variable.py is correct. I've also removed the functions which have been converted to VariableCoders and adapted the tests.

To sum up this PR, it does:

• convert functions to VariableCoders along @shoyer's TODO: https://github.com/pydata/xarray/blob/1c81162755457b3f4dc1f551f0321c75ec9daf6c/xarray/conventions.py#L298-L302 https://github.com/pydata/xarray/blob/1c81162755457b3f4dc1f551f0321c75ec9daf6c/xarray/conventions.py#L393-L405
• preserve boolean dtype within encoding: https://github.com/pydata/xarray/issues/7652#issuecomment-1476956975
• deterrmine cf packed dtype from scale_factor/add_offset

7691, #2304

, the PR seems to solve my specific issue without changing the encoding

Great, thanks for testing.

@euronion There is a potential fix for your issue in #7654. It would be great, if you could have a closer look and test against that PR.

@dcherian @basnijholt

After the dev-meeting I've taken a step back and first implemented the coders as mentioned in @shoyer's ToDo.

I've fixed the one bool->int issue and it now derives the dtype for ScaleOffset coding from scale_factor add_offset.

I've improved some test with regard to the scale/offset issue.

I'll concentrate on the string fillvalue issues in a follow up PR.

MCVE:

python fname = "test-7691.nc" import netCDF4 as nc with nc.Dataset(fname, "w") as ds0: ds0.createDimension("t", 5) ds0.createVariable("x", "int16", ("t",), fill_value=-32767) v = ds0.variables["x"] v.set_auto_maskandscale(False) v.add_offset = 278.297319296597 v.scale_factor = 1.16753614203674e-05 v[:] = np.array([-32768, -32767, -32766, 32767, 0]) with nc.Dataset(fname) as ds1: x1 = ds1["x"][:] print("netCDF4-python:", x1.dtype, x1) with xr.open_dataset(fname) as ds2: x2 = ds2["x"].values ds2.to_netcdf("test-7691-01.nc") print("xarray first read:", x2.dtype, x2) with xr.open_dataset("test-7691-01.nc") as ds3: x3 = ds3["x"].values print("xarray roundtrip:", x3.dtype, x3)

python netCDF4-python: float64 [277.9147410535744 -- 277.9147644042972 278.67988586425815 278.297319296597] xarray first read: float32 [277.91476 nan 277.91476 278.6799 278.29733] xarray roundtrip: float32 [ nan nan nan 278.6799 278.29733] I've confirmed that correctly promoting to float64 in CFMaskCoder solves this issue.

As this doesn't surface that often it might just happen here by accident. If the _FillValue/missing_value would be -32768 then the issue would not manifest.

So for NetCDF the default fillvalue for NC_SHORT (int16) is -32767. That means the promotion to float32 instead the needed float64 is the problem here (floating point precision).