Dimension order impact on simple workflow (dataset -> sel -> to_dataframe)

[theomasson]

Hi everybody, we encountered an issue with the use of dataset, I hope it can be helpful, with my colleagues we think it's more relevant in bug than in the help section. Thank you in advance and sorry if it's not in the good place.

What happened:
We use netcdf4 to store large datasets on a google bucket. These dataset have at least 3 dimension (latitude, longitude, time) and sometimes more. The size of files is often several Go with 3 to 100 variables.
In order to comply with netcdf4 recommendation, all our recent files contain variables with the following dimension : [time, latitude, longitude]
From these netcdf4, we usually want to extract a point (in latitude/longitude) for all times and convert it to a dataframe to use it. So the workflow is pretty simple:

data = xr.open_dataset('my_data.nc')
point = data.sel(latitude=40.163, longitude=9.569, method="nearest")
df = point.to_dataframe()

When y made this operation on our historic of files, the time to do these three step was very low (3-5ms) on a part of the files, and very very slow on the other part (60-120s !!). The only thing which change between these files is the fact that the first have [latitude, longitude, time] as variable dimension when the other ones (the slower ones) have [time, latitude, longitude]

More precisely, when trying to reproduce the problem localy (without distant bucket) I found a coeff 3 in time between these two option (3 variable dataset, 1.7Go file), so we keep a ms time to request a point. But if I use a low-speed device or again the distant bucket, the time increase drastically. From my observation it seems that in the second case the to_dataframe part (or a .values also) load the entire dataset. So with distant files it download the entire file, which is coherent with colleague comparison of the time (for low speed connection my 50s can be converted to 10-15min easily).

Test where made with decode_times=True,decode_cf=True in the open_dataset, without improvement. The only thing which actually work is to transpose dimension after opening the dataset. (Edit, actually did not work as expected, should have miss something)

What you expected to happen:
Same request time with both dimension order. Not downloading file to extract one point.

Minimal Complete Verifiable Example:
The MCVE is not easy as test are more relevant with a distant storage.

import xarray as xr
import time
import numpy as np
import pandas as pd
import datetime


latitude = np.linspace(-60, 60, 600)
longitude = np.linspace(-180, 180, 1200)
times = pd.date_range(start=datetime.datetime(2021,1,1,0,0),periods=100,freq='1H')
#Time in third
var1 = np.random.random((latitude.shape[0],longitude.shape[0],times.shape[0]))
var2 = np.random.random((latitude.shape[0],longitude.shape[0],times.shape[0]))
var3 = np.random.random((latitude.shape[0],longitude.shape[0],times.shape[0]))
dataset1 = xr.Dataset(
            {
                'var1': (
                    ['latitude', 'longitude', 'time'],
                    var1,
                    {'units': 'W.m**-2', 'long_name': 'Global Horizontal Irradiance'},
                ),
                'var2': (
                    ['latitude', 'longitude', 'time'],
                    var2,
                    {'units': 'W.m**-2', 'long_name': 'Diffuse Horizontal Irradiance'},
                ),
                'var3': (
                    ['latitude', 'longitude', 'time'],
                    var3,
                    {'units': 'W.m**-2', 'long_name': 'Direct Normal Irradiance'},
                ),
            },
            coords={
                'latitude': (['latitude'], latitude, {'units': 'degrees_north', 'long_name': 'Latitude'}),
                'longitude': (['longitude'], longitude, {'units': 'degrees_east', 'long_name': 'Longitude'}),
                'time': times,
            },
        )
dataset1.to_netcdf('/bucket/time_third.nc')


#Time in first
var1 = np.random.random((times.shape[0],latitude.shape[0],longitude.shape[0]))
var2 = np.random.random((times.shape[0],latitude.shape[0],longitude.shape[0]))
var3 = np.random.random((times.shape[0],latitude.shape[0],longitude.shape[0]))

dataset2 = xr.Dataset(
            {
                'var1': (
                    ['time','latitude', 'longitude'],
                    var1,
                    {'units': 'W.m**-2', 'long_name': 'Global Horizontal Irradiance'},
                ),
                'var2': (
                    ['time','latitude', 'longitude'],
                    var2,
                    {'units': 'W.m**-2', 'long_name': 'Global Horizontal Irradiance'},
                ),
                'var3': (
                    ['time','latitude', 'longitude'],
                    var3,
                    {'units': 'W.m**-2', 'long_name': 'Global Horizontal Irradiance'},
                ),
            },
            coords={
                'latitude': (['latitude'], latitude, {'units': 'degrees_north', 'long_name': 'Latitude'}),
                'longitude': (['longitude'], longitude, {'units': 'degrees_east', 'long_name': 'Longitude'}),
                'time': times,
            },
        )
dataset2.to_netcdf('/bucket/time_first.nc')



### Load and extract
print('time in third (fast)')
tps = time.time()
path_new_format = '/bucket/time_third.nc'
data_new_format = xr.open_dataset(path_new_format)
site_new_format = data_new_format.sel(latitude=42.1732, longitude=9.429554, method='nearest')
print('open + sel ' + str(time.time()- tps))
tps = time.time()
dataframe_new_format = site_new_format.to_dataframe()
print(' to dataframe ' + str(time.time()- tps))

print('time in first (slow)')
tps = time.time()
path_new_format = '/bucket/time_first.nc'
data_new_format = xr.open_dataset(path_new_format,decode_times=True,decode_cf=True)
site_new_format = data_new_format.sel(latitude=42.1732, longitude=9.429554, method='nearest')
print('open + sel ' + str(time.time()- tps))
tps = time.time()
dataframe_new_format = site_new_format.to_dataframe()
print(' to dataframe ' + str(time.time()- tps))

Results for local file (confirmed with a mean over 100 iteration):
time in third (fast)
sel 0.0072247982025146484
to dataframe 0.0017933845520019531
time in first (slow)
sel 0.007394075393676758
to dataframe 0.0036745071411132812

Results for bucket file (my connection is around 30Mo/s for downloading):
time in third (fast)
sel 0.5228869915008545
to dataframe 0.45035338401794434
time in first (slow)
sel 0.5486464500427246
to dataframe 54.501604318618774

Environment:

Output of xr.show_versions()

INSTALLED VERSIONS

commit: None
python: 3.7.7 (default, Mar 26 2020, 15:48:22)
[GCC 7.3.0]
python-bits: 64
OS: Linux
OS-release: 5.4.0-77-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: fr_FR.UTF-8
LOCALE: fr_FR.UTF-8
libhdf5: 1.10.4
libnetcdf: 4.6.2

xarray: 0.18.0
pandas: 1.2.4
numpy: 1.19.2
scipy: 1.4.1
netCDF4: 1.5.1.2
pydap: None
h5netcdf: 0.11.0
h5py: 2.10.0
Nio: None
zarr: 2.8.1
cftime: 1.5.0
nc_time_axis: None
PseudoNetCDF: None
rasterio: 1.0.21
cfgrib: 0.9.8.2
iris: None
bottleneck: None
dask: 2021.06.2
distributed: 2021.06.2
matplotlib: 3.1.3
cartopy: None
seaborn: None
numbagg: None
pint: None
setuptools: 52.0.0.post20210125
pip: 21.1.2
conda: None
pytest: 6.2.4
IPython: 7.19.0
sphinx: None

[max-sixty]

Hi @theomasson

I'm not that experienced in this area but I'll try an respond briefly.

I get the same results as you for the local files. That does suggests the slowness is unrelated to the dimension order — is there anything else that's different between the files?

From my observation it seems that in the second case the to_dataframe part (or a .values also) load the entire dataset. So with distant files it download the entire file, which is coherent with colleague comparison of the time (for low speed connection my 50s can be converted to 10-15min easily).

IIUC, open_dataset with no chunks argument will load the whole file. Are you seeing something different?

[max-sixty]

Can you show the code and output of open_dataset in the code that has the problem?

A reproduction of the error is probably required though. Unless I'm misunderstanding something, this reproduction doesn't reproduce the problem you're experiencing, is that correct?

[theomasson]

The provided code allow you to reproduce the problem, but to really see it it's easier to use a distant storage (like a google cloud or amazon bucket). Even so if you save de files in local you can see in the provided result a coeff 7 between the dataset with time in third (faster) and in first (slower). Using the bucket just allow you to have time in second and not in millisecond to compare (but the gap increase to a 120 time slower for the dataset saved in a distant filesystem).

I put two operational file (950Mo each) on my drive you can download locally, you will find the same order (6-7 on my computer with ssd) for a loc over a coordinate then a conversion to a datafram. If you have a distant file system you can try with the same files (I cannot gave access to my company gcpproject obviously). Also be careful, if you try several time, the file can be in memory and you face completely different results (it's very tricky to test, all we know is that the first load, which is critical for us, can be very slow)
First one (slow) : https://drive.google.com/file/d/18_Xl6fYicLvgGkfCucqiwsvLZ6Bu6x-C/view?usp=sharing
Second one (fast) : https://drive.google.com/file/d/1mdp0OwiRbL5hDu8RxsxhEzLfxo8M_cS8/view?usp=sharing

import xarray as xr
data = xr.open_dataset(path_file)
site = data.sel(latitude=42.1732, longitude=9.429554, method='nearest')
dataframe = site.to_dataframe()

I know that locally the gap (1ms to 7ms) look very low, but we made several thousand of sel/to_dataframe over different files. If locally a 7 time coeff is important, imagine when we face 60s per point with distant files ...

[max-sixty]

I get the same results on the repro:

In [9]: %%timeit
   ...: data_new_format = xr.open_dataset('temp/time_third.nc')
   ...: data_new_format.sel(latitude=42.1732, longitude=9.429554, method='nearest')
   ...: site_new_format.to_dataframe()
3.96 ms ± 385 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [10]: %%timeit
    ...: data_new_format = xr.open_dataset('temp/time_first.nc')
    ...: data_new_format.sel(latitude=42.1732, longitude=9.429554, method='nearest')
    ...: site_new_format.to_dataframe()
3.74 ms ± 293 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Do you have a repro that shows a difference over many runs (i.e. use something like %%timeit, not a one-off), and doesn't involve downloading data?

[theomasson]

Sorry for the late answer, I made 1000 loop over the files that I have send, but also over the file that you can generate with the first code. I found some differences in the request time :

first file (slow) :0.009751518353981737
second file (fast)   :0.006421096535977995
time_first :0.008170099526993
time_third :0.006355700689036894

But all of this again is time with file on local filesystem. Here the file with time in first dimension are about 30% slower, when using distant storage it can be so worse. I don't see how I can made a access to distant file for the example, I will see if I can made a bucket and giving you credential to test the real bad situation

[max-sixty]

Please can you share copy & pastable code that generates that result?

I'm not the expert on this area, and I'm not sure how nc files are accessed — others will know more. Possibly the order of the dimensions impacts whether the relevant parts of the file are contiguous — and so would require different numbers of reads.

I'm not sure a reproducible example is sufficient to evaluate this, but it would definitely be necessary. Thanks.