@headtr1ck I just discovered that it is not per se a difference between floats/da, but it has to do with the creation of the new dimension (plipoints in this case), I have updated the MCVE in https://github.com/pydata/xarray/issues/7701#issue-1647883619.

Thanks for your feedback, that is interesting and helpful. I have tested the older xarray version on a laptop with an older environment. I assumed the xarray version was the difference, but I guess there is something else that is causing it if you cannot reproduce it.

Environment where it does work as expected

Since .interp() is where it goes wrong, my next hunch was to check the scipy versions. The working env contains scipy 1.9.1, the failing environment contains scipy 1.10.1. When I downgrade to scipy 1.9.1 (or 1.9.3), the behaviour is as expected independent of the xarray version (.interp() gives same results with floats/da as input). When I use scipy 1.10.0 the behaviour is unexpected (different .interp() results), just like scipy 1.10.1. I see there are many scipy.interpolate changes done in scipy with the scipy 1.10.0 release, so that could explain it. For completeness, the interp_with_da variable changes with this scipy version change, not the interp_with_float variable, so it is still a bit xarray related.

I see no related+recent scipy issues yet. Do you have a suggestion on how to proceed?

@headtr1ck the proposed changes were pushed agian. I think most should be fixed now.

@rjp23 with the latest update to the PR (thanks to @jklymak), your example code produced identical figures without changing it.

Thanks to @jklymak, there was an update in PR I created (https://github.com/pydata/xarray/pull/7553). @ghiggi with the code from this PR, your code shows identical plots (except for the first one, but that should be the case). Hopefully the PR can be merged somewhere soon.

Thanks @jklymak, that makes sense I guess. I have committed the changes again. Could you check if it can be merged?

Hi @etsmith14. The suggestion I did loses accuracy and depending on the variable this is not acceptable. However, recomputing scale_factor and add_offset is possible: https://github.com/ArcticSnow/TopoPyScale/issues/60#issuecomment-1460022033 It is more complicated than dropping the dtype, but it does keep the filesize small.

@etsmith14: another workaround is removing the scale_factor instead of the dtype. This keeps the file size small. However, there are slight offsets between the source and destination datasets, which is to be expected since the original value for the msl variable was in the range of 0.1/0.11 and removing it defaults to 1. For your variable, the scale_factor might also be completely different. However, maybe the scale_factor (and add_offset) can be replaced by something that works for all ERA5 data instead of a value very specific to a single dataset/period.

@rjp23: could you close the issue if this indeed resolves your problem?

@ghiggi: If I understand it correctly, your issue/examplecode covers multiple issues. Since one subissue might be using norm instead levels, I would recommend trimming down your example code do only show the remaining vmax issue.

The related issues https://github.com/pydata/xarray/issues/4061 and https://github.com/Deltares/xugrid/issues/49 are fixed by supplying levels=levels instead of norm=norm to ds.plot(), as suggested by https://github.com/pydata/xarray/pull/7553#discussion_r1117264787. However, I cannot judge if this also fixes the example from this issue, since the plots still all look different. However, commenting #da1.data[da1.data>=norm.vmax] = norm.vmax - 1 solves most of the differences. @ghiggi: could you check if this suggestion solves your issue indeed?

As suggested by https://github.com/pydata/xarray/pull/7553#discussion_r1117264787, pass levels=bounds instead of norm=norm to data.plot(). Your example code results in a plot as expected. Would this solve your issue?

``` import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as colors import xarray

example from https://matplotlib.org/3.1.1/tutorials/colors/colormapnorms.html

for colormap normalisation

N = 100 X, Y = np.mgrid[-3:3:complex(0, N), -2:2:complex(0, N)] Z1 = np.exp(-X2 - Y2) Z2 = np.exp(-(X - 1)2 - (Y - 1)2) Z = (Z1 - Z2) * 2

fig, ax = plt.subplots(2, 1, figsize=(8, 8)) ax = ax.flatten()

bounds = np.linspace(-1, 1, 10) norm = colors.BoundaryNorm(boundaries=bounds, ncolors=256)

ax[0].pcolormesh(X, Y, Z, norm=norm, cmap='RdBu_r')

now add data into dataset and plot it using same normalisation

data = xarray.DataArray(Z, dims=('x', 'y'), coords={'x': X[:,0], 'y': Y[0,:]}) data.plot(ax=ax[1], x='x', y='y', levels=bounds, add_colorbar=False) ```

I just combined @Huite's suggestion with splitting the if-statement. This works for both solving the issue and keeping the testcases in test_plot.py green. if levels is not None: cmap, newnorm = _build_discrete_cmap(cmap, levels, extend, filled) norm = newnorm if norm is None else norm if isinstance(norm, mpl.colors.BoundaryNorm): cmap, norm = _build_discrete_cmap(cmap, levels, extend, filled) This could replace this code: https://github.com/pydata/xarray/blob/1de881e5f5d4d7697403f9ab7686849ab02953bd/xarray/plot/utils.py#L307

However, a bit up in the code there is a if isinstance(norm, mpl.colors.BoundaryNorm) statement, which I guess could be combined. https://github.com/pydata/xarray/blob/1de881e5f5d4d7697403f9ab7686849ab02953bd/xarray/plot/utils.py#L277 if isinstance(norm, mpl.colors.BoundaryNorm): levels = norm.boundaries

I think it is a potential solution nevertheless, but some help is appreciated with the last steps. Also since the case of @ghiggi seems not to be solved with this fix. It does solve https://github.com/pydata/xarray/issues/4061 though.

I have been thinking about a desireable solution, but I have a bit of trouble with it. Besides removing dtype from encoding (resulting in floats being written), one could also change the scale_factor to a higher value (e.g. 0.5). Writing this to int does take half the disksize than releasing the int restriction and writing it to float32. Whatever you do, the data is altered at least slightly.

Apparently, the data cannot be properly written to integers after reading it. This is a bit odd I would say, would that mean that the scaling+offset of ERA5 data is that thightly chosen that when applying it to another dataset/month, the data would fall out of the integer reach? Would be great if this would "just work". At the moment, apparently reading and writing ERA5 data with xarray results in incorrect netcdf files. I expected xarray would work off the shelf with these type of data, it feels like xarray is designed for doing exactly these type of things.

I have also encountered an issue with reading of ERA5 data with open_mfdatset, writing it to_netcdf() and reading it again (https://github.com/Deltares/dfm_tools/issues/239). I was actually looking for a place to land this, and found your issue.

My expectation is that this is because the ERA5 data is saved as ints, but all files have different offsets/scalingfactors. Upon opening it with open_mfdataset(), the data is converted to floats and to the offset/scalingfactor of the first file. This is fine, but the issue occurs I think (and what you also mention) since {'dtype': 'int16'} is in the encoding. The file is written as ints and this seems to mess up the data. (all a theory)

A workaround is to remove the dtype from the encoding for all variables in the file (or update to float32), but that seems cumbersome.

Is there any update on this issue? I have been running into the same problem recently and am happy to see that this issue was already recognized by others.

Thanks for pointing me there