Skip to content

slow performance with open_mfdataset #1385

Closed

Description

@rabernat

We have a dataset stored across multiple netCDF files. We are getting very slow performance with open_mfdataset, and I would like to improve this.

Each individual netCDF file looks like this:

%time ds_single = xr.open_dataset('float_trajectories.0000000000.nc')
ds_single
CPU times: user 14.9 ms, sys: 48.4 ms, total: 63.4 ms
Wall time: 60.8 ms

<xarray.Dataset>
Dimensions:  (npart: 8192000, time: 1)
Coordinates:
  * time     (time) datetime64[ns] 1993-01-01
  * npart    (npart) int32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
Data variables:
    z        (time, npart) float32 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 ...
    vort     (time, npart) float32 -9.71733e-10 -9.72858e-10 -9.73001e-10 ...
    u        (time, npart) float32 0.000545563 0.000544884 0.000544204 ...
    v        (time, npart) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
    x        (time, npart) float32 180.016 180.047 180.078 180.109 180.141 ...
    y        (time, npart) float32 -79.9844 -79.9844 -79.9844 -79.9844 ...

As shown above, a single data file opens in ~60 ms.

When I call open_mdsdataset on 49 files (each with a different time dimension but the same npart), here is what happens:

%time ds = xr.open_mfdataset('*.nc', )
ds
CPU times: user 1min 31s, sys: 25.4 s, total: 1min 57s
Wall time: 2min 4s

<xarray.Dataset>
Dimensions:  (npart: 8192000, time: 49)
Coordinates:
  * npart    (npart) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
  * time     (time) datetime64[ns] 1993-01-01 1993-01-02 1993-01-03 ...
Data variables:
    z        (time, npart) float64 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 ...
    vort     (time, npart) float64 -9.717e-10 -9.729e-10 -9.73e-10 -9.73e-10 ...
    u        (time, npart) float64 0.0005456 0.0005449 0.0005442 0.0005437 ...
    v        (time, npart) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
    x        (time, npart) float64 180.0 180.0 180.1 180.1 180.1 180.2 180.2 ...
    y        (time, npart) float64 -79.98 -79.98 -79.98 -79.98 -79.98 -79.98 ...

It takes over 2 minutes to open the dataset. Specifying concat_dim='time' does not improve performance.

Here is %prun of the open_mfdataset command.

         748994 function calls (724222 primitive calls) in 142.160 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       49   62.455    1.275   62.458    1.275 {method 'get_indexer' of 'pandas.index.IndexEngine' objects}
       49   47.207    0.963   47.209    0.963 base.py:1067(is_unique)
      196    7.198    0.037    7.267    0.037 {operator.getitem}
       49    4.632    0.095    4.687    0.096 netCDF4_.py:182(_open_netcdf4_group)
      240    3.189    0.013    3.426    0.014 numeric.py:2476(array_equal)
       98    1.937    0.020    1.937    0.020 {numpy.core.multiarray.arange}
4175/3146    1.867    0.000    9.296    0.003 {numpy.core.multiarray.array}
       49    1.525    0.031  119.144    2.432 alignment.py:251(reindex_variables)
       24    1.065    0.044    1.065    0.044 {method 'cumsum' of 'numpy.ndarray' objects}
       12    1.010    0.084    1.010    0.084 {method 'sort' of 'numpy.ndarray' objects}
5227/4035    0.660    0.000    1.688    0.000 collections.py:50(__init__)
       12    0.600    0.050    3.238    0.270 core.py:2761(insert)
12691/7497    0.473    0.000    0.875    0.000 indexing.py:363(shape)
   110728    0.425    0.000    0.663    0.000 {isinstance}
       12    0.413    0.034    0.413    0.034 {method 'flatten' of 'numpy.ndarray' objects}
       12    0.341    0.028    0.341    0.028 {numpy.core.multiarray.where}
        2    0.333    0.166    0.333    0.166 {pandas._join.outer_join_indexer_int64}
        1    0.331    0.331  142.164  142.164 <string>:1(<module>)

It looks like most of the time is being spent on reindex_variables. I understand why this happens...xarray needs to make sure the dimensions are the same in order to concatenate them together.

Is there any obvious way I could improve the load time? For example, can I give a hint to xarray that this reindex_variables step is not necessary, since I know that all the npart dimensions are the same in each file?

Possibly related to #1301 and #1340.

Activity

  1. rabernat commented on Apr 26, 2017

    @rabernat
    ContributorAuthor

    cc: @geosciz, who is helping with this project.

  2. shoyer commented on Apr 26, 2017

    @shoyer
    Member

    For example, can I give a hint to xarray that this reindex_variables step is not necessary

    Yes, adding an boolean argument prealigned which defaults to False to concat seems like a very reasonable optimization here.

    But more generally, I am a little surprised by how slow pandas.Index.get_indexer and pandas.Index.is_unique are. This suggests we should add a fast-path optimization to skip these steps in reindex_variables:

    if not index.is_unique:
    raise ValueError(
    'cannot reindex or align along dimension %r because the '
    'index has duplicate values' % name)
    indexer = get_indexer(index, target, method, tolerance)

    Basically, if index.equals(target), we should just set indexer = np.arange(target.size). Although, if we have duplicate values in the index, the operation should arguably fail for correctness.

  3. rabernat commented on Mar 2, 2018

    @rabernat
    ContributorAuthor

    An update on this long-standing issue.

    I have learned that open_mfdataset can be blazingly fast if decode_cf=False but extremely slow with decode_cf=True.

    As an example, I am loading a POP datataset on cheyenne. Anyone with access can try this example.

    base_dir = '/glade/scratch/rpa/'
    prefix = 'BRCP85C5CN_ne120_t12_pop62.c13b17.asdphys.001'
    code = 'pop.h.nday1.SST'
    glob_pattern = os.path.join(base_dir, prefix, '%s.%s.*.nc' % (prefix, code))
    
    def non_time_coords(ds):
        return [v for v in ds.data_vars
                if 'time' not in ds[v].dims]
    
    def drop_non_essential_vars_pop(ds):
        return ds.drop(non_time_coords(ds))   
    
    # this runs almost instantly
    ds = xr.open_mfdataset(glob_pattern, decode_times=False, chunks={'time': 1},
                           preprocess=drop_non_essential_vars_pop, decode_cf=False)

    And returns this

    <xarray.Dataset>
    Dimensions:     (d2: 2, nlat: 2400, nlon: 3600, time: 16401, z_t: 62, z_t_150m: 15, z_w: 62, z_w_bot: 62, z_w_top: 62)
    Coordinates:
      * z_w_top     (z_w_top) float32 0.0 1000.0 2000.0 3000.0 4000.0 5000.0 ...
      * z_t         (z_t) float32 500.0 1500.0 2500.0 3500.0 4500.0 5500.0 ...
      * z_w         (z_w) float32 0.0 1000.0 2000.0 3000.0 4000.0 5000.0 6000.0 ...
      * z_t_150m    (z_t_150m) float32 500.0 1500.0 2500.0 3500.0 4500.0 5500.0 ...
      * z_w_bot     (z_w_bot) float32 1000.0 2000.0 3000.0 4000.0 5000.0 6000.0 ...
      * time        (time) float64 7.322e+05 7.322e+05 7.322e+05 7.322e+05 ...
    Dimensions without coordinates: d2, nlat, nlon
    Data variables:
        time_bound  (time, d2) float64 dask.array<shape=(16401, 2), chunksize=(1, 2)>
        SST         (time, nlat, nlon) float32 dask.array<shape=(16401, 2400, 3600), chunksize=(1, 2400, 3600)>
    Attributes:
        nsteps_total:  480
        tavg_sum:      64800.0
        title:         BRCP85C5CN_ne120_t12_pop62.c13b17.asdphys.001
        start_time:    This dataset was created on 2016-03-14 at 05:32:30.3
        Conventions:   CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-curren...
        source:        CCSM POP2, the CCSM Ocean Component
        cell_methods:  cell_methods = time: mean ==> the variable values are aver...
        calendar:      All years have exactly  365 days.
        history:       none
        contents:      Diagnostic and Prognostic Variables
        revision:      $Id: tavg.F90 56176 2013-12-20 18:35:46Z mlevy@ucar.edu $
    

    This is roughly 45 years of daily data, one file per year.

    Instead, if I just change decode_cf=True (the default), it takes forever. I can monitor what is happening via the distributed dashboard. It looks like this:
    image

    There are more of these open_dataset tasks then there are number of files (45), so I can only presume there are 16401 individual tasks (one for each timestep), which each takes about 1 s in serial.

    This is a real failure of lazy decoding. Maybe it can be fixed by #1725, possibly related to #1372.

    cc Pangeo folks: @jhamman, @mrocklin

  4. shoyer commented on Mar 2, 2018

    @shoyer
    Member

    @rabernat How does performance compare if you call xarray.decode_cf() on the opened dataset? The adjustments I recently did to lazy decoding should only help once the data is already loaded into dask.

  5. rabernat commented on Mar 9, 2018

    @rabernat
    ContributorAuthor

    Calling ds = xr.decode_cf(ds, decode_times=False) on the dataset returns instantly. However, the variable data is wrapped in the adaptors, effectively destroying the chunks

    >>> ds.SST.variable._data
    LazilyIndexedArray(array=DaskIndexingAdapter(array=dask.array<_apply_mask, shape=(16401, 2400, 3600), dtype=float32, chunksize=(1, 2400, 3600)>), key=BasicIndexer((slice(None, None, None), slice(None, None, None), slice(None, None, None))))

    Calling getitem on this array triggers the whole dask array to be computed, which would takes forever and would completely blow out the notebook memory. This is because of #1372, which would be fixed by #1725.

    This has actually become a major showstopper for me. I need to work with this dataset in decoded form.

    Versions

    Details INSTALLED VERSIONS ------------------ commit: None python: 3.6.4.final.0 python-bits: 64 OS: Linux OS-release: 3.12.62-60.64.8-default machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: en_US.UTF-8

    xarray: 0.10.1
    pandas: 0.22.0
    numpy: 1.13.3
    scipy: 1.0.0
    netCDF4: 1.3.1
    h5netcdf: 0.5.0
    h5py: 2.7.1
    Nio: None
    zarr: 2.2.0a2.dev176
    bottleneck: 1.2.1
    cyordereddict: None
    dask: 0.17.1
    distributed: 1.21.3
    matplotlib: 2.1.2
    cartopy: 0.15.1
    seaborn: 0.8.1
    setuptools: 38.4.0
    pip: 9.0.1
    conda: None
    pytest: 3.3.2
    IPython: 6.2.1

  6. shoyer commented on Mar 9, 2018

    @shoyer
    Member

    OK, so it seems that we need a change to disable wrapping dask arrays with LazilyIndexedArray. Dask arrays are already lazy!

  7. shoyer commented on Nov 10, 2018

    @shoyer
    Member

    Was this fixed by #2047?

  8. chuaxr commented on Nov 15, 2018

    @chuaxr

    I can confirm that

    ds = xr.open_mfdataset(data_fnames,chunks={'lat':20,'time':50,'lon':24,'pfull':11},\
                          decode_cf=False)
    ds = xr.decode_cf(ds)
    

    is much faster (seconds vs minutes) than

    ds = xr.open_mfdataset(data_fnames,chunks={'lat':20,'time':50,'lon':24,'pfull':11})
    
    

    . For reference, data_fnames is a list of 5 files, each of which is ~75 GB.

  9. shoyer commented on Nov 15, 2018

    @shoyer
    Member

    @chuaxr I assume you're testing this with xarray 0.11?

    It would be good to do some profiling to figure out what is going wrong here.

  10. chuaxr commented on Nov 15, 2018

    @chuaxr

    Yes, I'm on 0.11.

    Nothing displays on the task stream/ progress bar when using open_mfdataset, although I can monitor progress when, say, computing the mean.

    The output from %time using decode_cf = False is

    CPU times: user 4.42 s, sys: 392 ms, total: 4.82 s
    Wall time: 4.74 s
    

    and for decode_cf = True:

    CPU times: user 11.6 s, sys: 1.61 s, total: 13.2 s
    Wall time: 3min 28s
    

    Using xr.set_options(file_cache_maxsize=1) doesn't make any noticeable difference.

    If I repeat the open_mfdataset for another 5 files (after opening the first 5), I occasionally get this warning:
    distributed.utils_perf - WARNING - full garbage collections took 24% CPU time recently (threshold: 10%)

    I only began using the dashboard recently; please let me know if there's something basic I'm missing.

  11. shoyer commented on Nov 16, 2018

    @shoyer
    Member

    @chuaxr What do you see when you use %prun when opening the dataset? This might point to the bottleneck.

    One way to fix this would be to move our call to decode_cf() in open_dataset() to after applying chunking, i.e., to switch up the order of operations on these lines:

    def maybe_decode_store(store, lock=False):
    ds = conventions.decode_cf(
    store, mask_and_scale=mask_and_scale, decode_times=decode_times,
    concat_characters=concat_characters, decode_coords=decode_coords,
    drop_variables=drop_variables)
    _protect_dataset_variables_inplace(ds, cache)
    if chunks is not None:
    from dask.base import tokenize
    # if passed an actual file path, augment the token with
    # the file modification time
    if (isinstance(filename_or_obj, basestring) and
    not is_remote_uri(filename_or_obj)):
    mtime = os.path.getmtime(filename_or_obj)
    else:
    mtime = None
    token = tokenize(filename_or_obj, mtime, group, decode_cf,
    mask_and_scale, decode_times, concat_characters,
    decode_coords, engine, chunks, drop_variables)
    name_prefix = 'open_dataset-%s' % token
    ds2 = ds.chunk(chunks, name_prefix=name_prefix, token=token)
    ds2._file_obj = ds._file_obj
    else:
    ds2 = ds
    return ds2

    In practice, is the difference between using xarray's internal lazy array classes for decoding and dask for decoding. I would expect to see small differences in performance between these approaches (especially when actually computing data), but for constructing the computation graph I would expect them to have similar performance. It is puzzling that dask is orders of magnitude faster -- that suggests that something else is going wrong in the normal code path for decode_cf(). It would certainly be good to understand this before trying to apply any fixes.

  12. chuaxr commented on Nov 16, 2018

    @chuaxr

    Sorry, I think the speedup had to do with accessing a file that had previously been loaded rather than due to decode_cf. Here's the output of prun using two different files of approximately the same size (~75 GB), run from a notebook without using distributed (which doesn't lead to any speedup):

    Output of
    %prun ds = xr.open_mfdataset('/work/xrc/AM4_skc/atmos_level.1999010100-2000123123.sphum.nc',chunks={'lat':20,'time':50,'lon':12,'pfull':11})

    
              780980 function calls (780741 primitive calls) in 55.374 seconds
     
        Ordered by: internal time
     
        ncalls  tottime  percall  cumtime  percall filename:lineno(function)
             7   54.448    7.778   54.448    7.778 {built-in method _operator.getitem}
        764838    0.473    0.000    0.473    0.000 core.py:169(<genexpr>)
             3    0.285    0.095    0.758    0.253 core.py:169(<listcomp>)
             2    0.041    0.020    0.041    0.020 {cftime._cftime.num2date}
             3    0.040    0.013    0.821    0.274 core.py:173(getem)
             1    0.027    0.027   55.374   55.374 <string>:1(<module>)
    
    
    

    Output of
    %prun ds = xr.open_mfdataset('/work/xrc/AM4_skc/atmos_level.2001010100-2002123123.temp.nc',chunks={'lat':20,'time':50,'lon':12,'pfull':11},
    decode_cf=False)

     
              772212 function calls (772026 primitive calls) in 56.000 seconds
     
        Ordered by: internal time
     
        ncalls  tottime  percall  cumtime  percall filename:lineno(function)
             5   55.213   11.043   55.214   11.043 {built-in method _operator.getitem}
        764838    0.486    0.000    0.486    0.000 core.py:169(<genexpr>)
             3    0.185    0.062    0.671    0.224 core.py:169(<listcomp>)
             3    0.041    0.014    0.735    0.245 core.py:173(getem)
             1    0.027    0.027   56.001   56.001 <string>:1(<module>)
    

    /work isn't a remote archive, so it surprises me that this should happen.

  13. shoyer commented on Nov 16, 2018

    @shoyer
    Member
  14. 28 remaining items

  15. shoyer commented on Aug 28, 2023

    @shoyer
    Member

    @fraserwg do you know what the performance bottleneck is in your case? (i.e., if you interrupt the comptuation, what is being computed?)

  16. dcherian commented on Aug 31, 2023

    @dcherian
    Contributor

    It is very common for different netCDF files in a "dataset" (a folder) to be encoded differently so we can't set decode_cf=False by default.

    there's probably something else going on under the hood that's causing the slowness of open_mfdataset at present.

    There's

    1. slowness in reading coordinate information from every file. We have parallel to help a bit here.
    2. slowness in combining each file to form a single Xarray dataset. By default, we do lots of consistency checking by reading data. Xarray allows you to control this, data_vars='minimal', coords='minimal', compat='override' is a common choice.

    What your describing sounds like a failure of lazy decoding or acftime slowdown (example)which should be fixed. If you can provide a reproducible example, that would help.

  17. fraserwg commented on Sep 7, 2023

    @fraserwg

    @shoyer I will double check what the bottle neck is and report back. @dcherian interestingly, the parallel option doesn't seem to do anything when decode_cf=True. From looking at the dask dashboard it seems to load each file sequentially, with each opening being carried out by a different worker but not concurrently. I will see what I can do minimal example wise!

  18. ashjbarnes commented on Feb 21, 2024

    @ashjbarnes

    I'm another person who stumbled across this thread, and found that decode_cf = False fix to work really well.

    I appreciate that we can't set this by default, but maybe this could be put into the docstring of open_mfdataset directly? It appears to be passed as a kwarg, so is hard to find despite it being such a helpful fix for so many people in this thread!

    It also seems like the decode_cf step is done in serial. My dask cluster has plenty of workers but when decode_cf is set to True it only processes one of my (many) files at a time. Switching to decode_cf = False and the task stream shows my entire cluster being utilised. Perhaps this is part of the issue?

    xarray version: 2024.2.0

  19. dcherian commented on Feb 22, 2024

    @dcherian
    Contributor

    @ashjbarnes Are you able to share two small files that illustrate the issue?

  20. ashjbarnes commented on Feb 22, 2024

    @ashjbarnes

    @dcherian Thanks, the files are shared in this google drive folder.

    This is a spatial subsample of the same files. Perhaps it's a bad idea to store static data like the lat/lon coordinates in each file? The overall size of this is so tiny compared to the 4D data that I left the coordinates there for convenience but I'm not sure whether this has broader implications.

    In my tests, running:

    xr.open_mfdataset(PATH,decode_times = False,parallel = True,decode_cf = False)
    

    on ~3000 files of 300mb each had an order of magnitude speedup over the same command with decode_cf = True on xarray 2024.2.0. The real files are chunked on disk in time and the yb dimension

  21. dcherian commented on Feb 22, 2024

    @dcherian
    Contributor

    Wow, thank you!

    This is an amazing bug. The defaults say data_vars="all", coords="different" which means always concatenate all data_vars along the concat dimensions (here inferred to be "time") but only concatenate coords if they differ in the different files.

    When decode_cf=False , lat ,lon are data_vars and get concatenated without any checking or reading.
    When decode_cf=True, lat, lon are promoted to coords, then get checked for equality across all files. The two variables get read sequentially from all files. This is the slowdown you see.

    Once again, this is a consequence of bad defaults for concat and open_mfdataset.

    I would follow https://docs.xarray.dev/en/stable/user-guide/io.html#reading-multi-file-datasets and use data_vars="minimal", coords="minimal", compat="override" which will only concatenate those variables with the time dimension, and skip any checking for variables that don't have a time dimension (simply pick the variable from the first file).

  22. ashjbarnes commented on Feb 22, 2024

    @ashjbarnes

    Well done @dcherian great find! Changing the defaults does seem like a really good idea in this case

  23. hannahzafar commented on Mar 13, 2024

    @hannahzafar

    @dcherian I am yet another person stumbling on this problem. Unfortunately decode_cf = False seems to override decode_times=True (https://stackoverflow.com/questions/77243075/xarray-wont-decode-times-from-netcdf-file-with-decode-cf-false-even-if-decode) so you cannot use that fix if you want to maintain datetime objects.

  24. dcherian commented on Mar 13, 2024

    @dcherian
    Contributor

    A simpler fix is probably decode_coords=False

  25. rabernat commented on Mar 14, 2024

    @rabernat
    ContributorAuthor

    This issue is almost seven years old! It has been "fixed" many times since my original post, but people keep finding new ways to make it reappear. 😆

    It seems like having better diagnostics / logging of what is happening under the hood with open_mfdataset is what people really need. Maybe even some sort of utility to pre-scan the files and figure out if they are easily openable or not.

  26. jtomfarrar commented on Mar 14, 2024

    @jtomfarrar

    This issue is almost seven years old! It has been "fixed" many times since my original post, but people keep finding new ways to make it reappear. 😆

    It seems like having better diagnostics / logging of what is happening under the hood with open_mfdataset is what people really need. Maybe even some sort of utility to pre-scan the files and figure out if they are easily openable or not.

    Both of those seem like great ideas. Maybe there could be a verbose or logging mode to help users identify what is wrong with the files (e.g., where the time is being spent and whether that seems suspicious). It is probably true that people (like me) will keep finding new ways to generate problematic netcdf files. (I'm sure we can think of something even worse than 20 Hz data referenced to a time origin 75 years ago).

  27. added a commit that references this issue on Oct 29, 2024
  28. dcherian commented on Aug 1, 2025

    @dcherian
    Contributor

    We've just merged #10062 from @jsignell , so after the next release (v2025.08.0) you should be able to use xr.set_options(use_new_combine_kwarg_defaults=True) to get much more consistent behaviour that does not load data by default.

    If you're stopping by this issue, please try your workload with this setting.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions