Why is loading data from a file so slow?

[krikru]

I have a netCDF file that is derived from the ERA-5 dataset, originally downloaded using the era5cli command line tool. When I use ncview to view the contents of the file, and I press the "fast forward" button to play back the data time point per time point, I can play back approximately 1 month per second, which corresponds to about 30 time points per second (this specific file contains daily data).

However, when I load the exact same file into an iris cube with iris.load_cube, then load each of the first ten time slices into a NumPy array (cube[idx:idx+1].data where I loop idx from 0 to 9), the loop takes 6.43 seconds, resulting in an effective frame rate of 1.56, which is far from the estimated 30 frames per second obtained with ncview.

How comes the data is loaded so slowly with iris? It should be possible to load it much quicker, judging from how fast ncview can play back the content, so how comes that it doesn't load that quickly?

(For anyone who is wondering, the netCDF file is chuncked, with each time point being one chunk (which is also what I'm loading separately into my script), has the format netCDF-4, and follows the CF-1.7 convention (whatever that is), according to ncdump. I think it is also compressed, but I can't see that in the ncdump output, and I can't find the script I used to produce the file at the moment.)

[trexfeathers]

Hi @krikru, thanks for getting in touch. We don't have an immediate explanation - we would need to do some investigation. What specifically are you downloading via era5cli? We can try and replicate your observations.

[trexfeathers]

Thanks for the detail @krikru. I didn't realise I also needed credentials to use these downloading API's, so I can't precisely replicate, but fellow developers have still pointed out a possible explanation...

If you print your xvar variable, you will see:

[... values with dtype=float32]

... or similar. This indicates that the data has ALREADY BEEN LOADED into memory - making all subsequent operations much faster. If you change your open_dataset() call to: open_dataset(file_path, chunks=-1, ...), you will see:

dask.array<open_dataset-...

... which shows that the Dataset still has lazy (Dask) data. Only then will the test be equivalent between Iris, Dask and Xarray.

[krikru]

Hi @trexfeathers and sorry for not getting back to you earlier.

When I print xvar, I indeed get the text [7405818768 values with dtype=float32].

But I don't think the time slice have already been loaded, since I'm loading them with xarray first and then loading them with iris and lastly with iris+dask (even though I mentioned them in the opposite order here). I am loading the netCDF file to a cube with iris first and then calling the core_data method of the loaded cube, but I don't think any of that loads anything more than metadata, because loading the entire cube would take several minutes, while doing this takes about 0.25 seconds (which is less than it takes to load the ten time slices with iris or with iris+dask).

If I re-run the tests immediately afterwards, by using the same xvar, dask_data and cube, the numbers are very similar, and if I run the tests without loading the time slices with iris and iris+dask at all, I'm also getting very a similar number for xarray.

[krikru]

@trexfeathers Yes, you need credentials to use the downloading API, but maybe you can reproduce my results approximately with other netCDF files too. I did generate the netCDF file that I'm running the test on with iris, after all. 🙂

[trexfeathers]

Hi @krikru, thanks for your patience.

I've finally had time to try replicating your observations, like you suggested.

I won't think any more about whether Xarray is loading lazily or not: in my script your observations hold true with both lazy AND eager loading.

The implementation of lazy loading is complex, so it is inevitable that Xarray and Iris have made different choices about how to do this. @pp-mo and I have spent an hour trying different things but have not isolated the cause of the performance difference, however I'm not especially surprised:

• Iris' implementation is simpler / less sophisticated
• Xarray has a larger contributer community and improves at a faster rate
• Some of Xarray's most prominent members are specifically interested in very large lazy operations
  Iris'/SciTools' main recent focuses have been CF-parsing (including UGRID) and regridding

Obviously this lack of performance is disappointing. We can try and restructure our approach to lazy data if a speed up is important to you and other users. However, as part of the Iris ❤ Xarray philosophy: we encourage people to use Xarray where it is the more appropriate tool, and we are also working hard on the ncdata project (still in development) to improve the experience of switching between Iris & Xarray to make use of their relative strengths. Using Xarray where appropriate is hopefully a way of you getting the performance you need immediately, without having to wait for future Iris development.

Expand for my code

from time import perf_counter

import iris
import numpy as np
import xarray


class CatchTime:
    # https://stackoverflow.com/a/69156219
    def __init__(self, message_prefix: str):
        self.readout = message_prefix

    def __enter__(self):
        self.time = perf_counter()
        return self

    def __exit__(self, type_, value, traceback):
        self.time = perf_counter() - self.time
        self.readout += f"Time: {self.time:.3f} seconds"
        print(self.readout)


file_path = iris.sample_data_path("A1B_north_america.nc")

with CatchTime("Iris loading: "):
    my_cube = iris.load_cube(file_path)
assert my_cube.has_lazy_data()

with CatchTime("Xarray loading: "):
    my_ds = xarray.open_dataset(file_path, chunks=-1)
my_var = my_ds.data_vars["air_temperature"]
assert not isinstance(my_var.variable._data, np.ndarray)

print("---")

with CatchTime("Iris realisation: "):
    for i in range(10):
        assert isinstance(
            my_cube[i:i+1].data,
            np.ndarray,
        )

with CatchTime("Xarray realisation: "):
    for i in range(10):
        assert isinstance(
            my_var.isel(time=slice(i, i+1)).values,
            np.ndarray,
        )

Iris loading: Time: 1.161 seconds
Xarray loading: Time: 0.141 seconds
---
Iris realisation: Time: 0.089 seconds
Xarray realisation: Time: 0.019 seconds

[trexfeathers]

Ping @pp-mo