Indexing issues with lazy loading

[jurjendejong]

Hi all,

First thanks for this great package! Started using it as an alternative Python backend recently, and impressed by how much it simplifies my work.

Recently I came across the issue that results might differ when lazy loading.
The .diff(dim='plan') does not produce correct results.
The problem does not appear with other axis (only 'plan'), and as soon as you trigger anything to actually load the data before applying the .diff() it resolves the issue.

Similar issues occur also with .shift(plan=1).

(sorry about the formating)

`
slf_open = xr.open_dataset(selafin_3D_file, engine="selafin")
slf_load = xr.load_dataset(selafin_3D_file, engine="selafin")

#%%

Get the delta_z between all the Z layers gives different results between .open_dataset() and .load_dataset()

z_levels_open = slf_open.Z.isel(time=0).drop_vars('time')
z_levels_load = slf_load.Z.isel(time=0).drop_vars('time')

dz_open = z_levels_open.diff(dim='plan')
dz_load = z_levels_load.diff(dim='plan')

xr.testing.assert_equal(dz_open, dz_load)

--> AssertionError

#%%

But if I query the Z-data itself all seems

xr.testing.assert_equal(z_levels_open, z_levels_load)

--> No AssertionError

#%%

If I test my original problem again, it has now disappeared

dz_open = z_levels_open.diff(dim='plan')
dz_load = z_levels_load.diff(dim='plan')

xr.testing.assert_equal(dz_open, dz_load)

--> No Assertion Error

#%%

Also if I request a .load() the issue has disappeared.

z_levels_open = slf_open.Z.isel(time=0).drop_vars('time').load()
z_levels_load = slf_load.Z.isel(time=0).drop_vars('time')

dz_open = z_levels_open.diff(dim='plan')
dz_load = z_levels_load.diff(dim='plan')

xr.testing.assert_equal(dz_open, dz_load)

--> No Assertion error`

[tomsail]

Thanks for this detailed traceback @jurjendejong
It must be that the data structure returned in Lazy Loading is different from the normal one:

xarray-selafin/xarray_selafin/xarray_backend.py

Line 300 in 6c99f28

data = indexing.LazilyIndexedArray(lazy_array)

I will look into this and adapt the tests accordingly

[tomsail]

@jurjendejong can you check the version you have for this error.
I want to make sure it's the last one (0.1.9) since #43 fixed a very similar problem

[tomsail]

I could indeed reproduce your case with the latest version.
However we talk about really small differences (of 10e-7 precision):

import hvplot.xarray
import xarray as xr
selafin_3D_file = "data/r3d_tidal_flats.slf"

slf_open = xr.open_dataset(selafin_3D_file, engine="selafin")
slf_load = xr.load_dataset(selafin_3D_file, engine="selafin")

# Get the delta_z between all the Z layers gives different results between .open_dataset() and .load_dataset()
z_levels_open = slf_open.Z.isel(time=0).drop_vars('time')
z_levels_load = slf_load.Z.isel(time=0).drop_vars('time')

dz_open = z_levels_open.diff(dim='plan')
dz_load = z_levels_load.diff(dim='plan')

dz_open.hvplot.scatter(x="x", y="y", c="Z")
dz_load.hvplot.scatter(x="x", y="y", c="Z")

Which should be intrinsically related to the Lazy Loading of the DataArray.

Instead of using xr.testing.assert_equal() you should use:

xr.testing.assert_allclose(dz_open, dz_load, rtol=10e-6)

Which fails for relative precision rtol=10e-7 and passes for rtol=10e-6

I guess I can close the issue since I don't think the problem is related to xarray-selafin

[jurjendejong]

Hej! Sorry for my terribly slow response. I had not expect a response this quick and apparently have no notifications :).

My issue was still valid with the latest version. I am not sure why your test did not reproduce it. Maybe the test model had uniform depth in all cells at timestep 0.

I have now created a pull request #55 that seems to fix the issue. It turned out to be a bit of a mess with the order of dimensions, and this gave errors when you started slicing over nodes/planes. Please run tests to validate that I did not break anything.

[tomsail]

Hi Jurjen, your implementation didn't pass CI. can you provide a selafin and minimal example of the bug?
I'll reopen the issue to trach this

[jurjendejong]

Good thing you have proper testing.
I cannot share my model, but will see if I can reproduce it with one of the test cases.

[jurjendejong]

I managed to reproduce it using this testmodel. I have included example and code.
I also believe that my latest commit might have resolved most errors.

import xarray as xr

inputfile = r'r3d_bump.slf'

slf_open = xr.open_dataset(inputfile, engine="selafin")
slf_load = xr.load_dataset(inputfile, engine="selafin")

z_levels_open = slf_open.Z.isel(time=0).drop_vars('time')
dz_open = z_levels_open.diff(dim='plan')

z_levels_load = slf_load.Z.isel(time=0).drop_vars('time')
dz_load = z_levels_load.diff(dim='plan')

xr.testing.assert_allclose(dz_open, dz_load,  rtol=1e-3)

r3d_bump.zip

[tomsail]

Thanks Jurjen for your commits, I try to be reactive to approve pushes so that they can be tested.
Note that a personal fork on your git would have allowed you to workaround this workflow approval, consider doing so in the future so that it more streamlined for all of us.

I'll have a look at the example you provided too.
I think it is relevant to add this case to the tests too:

• add the r3d_bump.slf in tests/data
• add it the files to test in tests/data/io_test.py

[jurjendejong]

Thanks for the input, much appreciated. I do believe that I work through a personal fork, but I might done something wrong...

Do you expect me to update the test cases? Not sure if the netcdf-read-write method that you use for many tests is as applicable in this case.

[tomsail]

I do believe that I work through a personal fork

Indeed sorry I missed it! You can thus activate CI tests in the Actions tab on your fork.

Do you expect me to update the test cases? Not sure if the netcdf-read-write method that you use for many tests is as applicable in this case.

I have a bit of overhead now, but from October on, I should have more time to dedicate to the packages. I'll review your work and the other PRs at that point.

[tomsail]

fixed with #56