OK, we might check if that depends on the data size or on the number of groups, or both.

It scales with data size it seems, but: even if you reduce data size to 1 element, after 50 iterations a single write goes up to 150ms already (whereas it's a few milliseconds in an empty file). These 150ms is the pure 'file traversal' etc part; the rest (of the 2 seconds) is the part where it seemingly reads stuff - which scales with data. Ideally it should just stay at <10ms all the time.

I wonder if it would help to use the same underlying h5py.File or h5netcdf.File when appending.

I don't think it's about what's happening in the current Python's process, which instances are being cached or not, it's about the general logic.

For instance, in the example above, if you run it once (e.g. set the range to 50); and then run it but comment out the block that clears the file, and set the range to 50-100. The very first dataset written the second time will be already very slow, slower than the last dataset written the first time - which means it's not about reusing the same File instance.

Without digging into implementational details, my logic as a library user would be this:

• If I write one dataset to file1 and another dataset to file2 using to_netcdf(), to different groups
• And then I simply combine the two hdf5 files using some external tools (again, datasets stored in different groups)
• I will be able to read them both perfectly well using open_dataset() or load_dataset()
• This implies that the datasets can be written just fine independently without knowing about each other
• Why then those writing functions (flush in particular) traverse and read the entire file every time?

Just checked it out.

| Number of datasets in file | netCDF4 (ms/write) | h5netcdf (ms/write) | | --- | --- | --- | | 1 | 4 | 11 | | 250 | 142| 1933 |

Here's the minimal example, try running this:

```python import time import xarray as xr import numpy as np import h5py

arr = xr.DataArray(np.random.RandomState(0).randint(-100, 100, (50_000, 3)), dims=['x', 'y']) ds = xr.Dataset({'arr': arr})

filename = 'test.h5' save = lambda group: ds.to_netcdf(filename, engine='h5netcdf', mode='a', group=str(group))

with h5py.File(filename, 'w') as _: pass

for i in range(250): t0 = time.time() save(i) print(time.time() - t0) ```

Is it possible to use .to_netcdf() without h5netcdf.File touching any of the pre-existing data or attempting to read it or traverse it? This will inevitably cause quadratic slowdowns as you write multiple datasets to the file - and that's what seems to be happening.

Or at least, don't traverse anything above the current root group that the dataset is being written into.

@kmuehlbauer Just installed h5netcdf=0.10.0, here's the timings when there's 200 groups in file - store.close() takes 92.4% of time again:

1078 1 1.0 1.0 0.0 try: 1079 # TODO: allow this work (setting up the file for writing array data) 1080 # to be parallelized with dask 1081 2 221642.0 110821.0 4.2 dump_to_store( 1082 1 2.0 2.0 0.0 dataset, store, writer, encoding=encoding, unlimited_dims=unlimited_dims 1083 ) 1084 1 3.0 3.0 0.0 if autoclose: 1085 store.close() 1086 1087 1 1.0 1.0 0.0 if multifile: 1088 return writer, store 1089 1090 1 6.0 6.0 0.0 writes = writer.sync(compute=compute) 1091 1092 1 1.0 1.0 0.0 if path_or_file is None: 1093 store.sync() 1094 return target.getvalue() 1095 finally: 1096 1 2.0 2.0 0.0 if not multifile and compute: 1097 1 4857912.0 4857912.0 92.6 store.close()

And here's _lookup_dimensions(): (note that it only takes half of the time, there's tons of other time spent in File.flush() which I don't understand):

``` Timer unit: 1e-06 s

Total time: 2.44857 s File: .../python3.8/site-packages/h5netcdf/core.py Function: _lookup_dimensions at line 92

Line # Hits Time Per Hit % Time Line Contents

92                                               def _lookup_dimensions(self):
93       400      65513.0    163.8      2.7          attrs = self._h5ds.attrs
94       400       6175.0     15.4      0.3          if "_Netcdf4Coordinates" in attrs:
95                                                       order_dim = _reverse_dict(self._parent._dim_order)
96                                                       return tuple(
97                                                           order_dim[coord_id] for coord_id in attrs["_Netcdf4Coordinates"]
98                                                       )
99

100 400 44938.0 112.3 1.8 child_name = self.name.split("/")[-1] 101 400 5006.0 12.5 0.2 if child_name in self._parent.dimensions: 102 return (child_name,) 103
104 400 350.0 0.9 0.0 dims = [] 105 400 781.0 2.0 0.0 phony_dims = defaultdict(int) 106 1400 166093.0 118.6 6.8 for axis, dim in enumerate(self._h5ds.dims): 107 # get current dimension 108 1000 119507.0 119.5 4.9 dimsize = self.shape[axis] 109 1000 2459.0 2.5 0.1 phony_dims[dimsize] += 1 110 1000 34345.0 34.3 1.4 if len(dim): 111 1000 2001071.0 2001.1 81.7 name = _name_from_dimension(dim) 112 else: 113 # if unlabeled dimensions are found 114 if self._root._phony_dims_mode is None: 115 raise ValueError( 116 "variable %r has no dimension scale " 117 "associated with axis %s. \n" 118 "Use phony_dims=%r for sorted naming or " 119 "phony_dims=%r for per access naming." 120 % (self.name, axis, "sort", "access") 121 ) 122 else: 123 # get dimension name 124 name = self._parent._phony_dims[(dimsize, phony_dims[dimsize] - 1)] 125 1000 1820.0 1.8 0.1 dims.append(name) 126 400 512.0 1.3 0.0 return tuple(dims) ```

There's some absolutely obscure things here, e.g. h5netcdf.core.BaseVariable._lookup_dimensions:

For 0 datasets:

``` Timer unit: 1e-06 s

Total time: 0.005034 s File: .../python3.8/site-packages/h5netcdf/core.py Function: _lookup_dimensions at line 86

Line # Hits Time Per Hit % Time Line Contents

86                                               def _lookup_dimensions(self):
87         2        633.0    316.5     12.6          attrs = self._h5ds.attrs
88         2         53.0     26.5      1.1          if '_Netcdf4Coordinates' in attrs:
89                                                       order_dim = _reverse_dict(self._parent._dim_order)
90                                                       return tuple(order_dim[coord_id]
91                                                                    for coord_id in attrs['_Netcdf4Coordinates'])
92                                           
93         2        471.0    235.5      9.4          child_name = self.name.split('/')[-1]
94         2         51.0     25.5      1.0          if child_name in self._parent.dimensions:
95                                                       return (child_name,)
96                                           
97         2          4.0      2.0      0.1          dims = []
98         7       1671.0    238.7     33.2          for axis, dim in enumerate(self._h5ds.dims):
99                                                       # TODO: read dimension labels even if there is no associated

100 # scale? it's not netCDF4 spec, but it is unambiguous... 101 # Also: the netCDF lib can read HDF5 datasets with unlabeled 102 # dimensions. 103 5 355.0 71.0 7.1 if len(dim) == 0: 104 raise ValueError('variable %r has no dimension scale ' 105 'associated with axis %s' 106 % (self.name, axis)) 107 5 1772.0 354.4 35.2 name = _name_from_dimension(dim) 108 5 18.0 3.6 0.4 dims.append(name) 109 2 6.0 3.0 0.1 return tuple(dims) ```

For 200 datasets:

``` Timer unit: 1e-06 s

Total time: 2.34179 s File: .../python3.8/site-packages/h5netcdf/core.py Function: _lookup_dimensions at line 86

Line # Hits Time Per Hit % Time Line Contents

86                                               def _lookup_dimensions(self):
87       400      66185.0    165.5      2.8          attrs = self._h5ds.attrs
88       400       6106.0     15.3      0.3          if '_Netcdf4Coordinates' in attrs:
89                                                       order_dim = _reverse_dict(self._parent._dim_order)
90                                                       return tuple(order_dim[coord_id]
91                                                                    for coord_id in attrs['_Netcdf4Coordinates'])
92                                           
93       400      45176.0    112.9      1.9          child_name = self.name.split('/')[-1]
94       400       5006.0     12.5      0.2          if child_name in self._parent.dimensions:
95                                                       return (child_name,)
96                                           
97       400        317.0      0.8      0.0          dims = []
98      1400     168708.0    120.5      7.2          for axis, dim in enumerate(self._h5ds.dims):
99                                                       # TODO: read dimension labels even if there is no associated

100 # scale? it's not netCDF4 spec, but it is unambiguous... 101 # Also: the netCDF lib can read HDF5 datasets with unlabeled 102 # dimensions. 103 1000 35653.0 35.7 1.5 if len(dim) == 0: 104 raise ValueError('variable %r has no dimension scale ' 105 'associated with axis %s' 106 % (self.name, axis)) 107 1000 2012597.0 2012.6 85.9 name = _name_from_dimension(dim) 108 1000 1640.0 1.6 0.1 dims.append(name) 109 400 400.0 1.0 0.0 return tuple(dims) ```

• All of the time in store.close() is, in its turn spent in CachingFileManager.close()
• That time is spent in h5netcdf.File.close()
• All of which is spent in h5netcdf.File.flush()

h5netcdf.File.flush() when there's 0 datasets in file:

``` 0.21619391441345215

Timer unit: 1e-06 s

Total time: 0.006862 s File: .../python3.8/site-packages/h5netcdf/core.py Function: flush at line 689

Line # Hits Time Per Hit % Time Line Contents

689 def flush(self): 690 1 4.0 4.0 0.1 if 'r' not in self._mode: 691 1 111.0 111.0 1.6 self._set_unassigned_dimension_ids() 692 1 3521.0 3521.0 51.3 self._create_dim_scales() 693 1 3224.0 3224.0 47.0 self._attach_dim_scales() 694 1 2.0 2.0 0.0 if not self._preexisting_file and self._write_ncproperties: 695 self.attrs._h5attrs['_NCProperties'] = _NC_PROPERTIES ```

h5netcdf.File.flush() when there's 200 datasets in file (758 times slower):

``` Timer unit: 1e-06 s

Total time: 4.55295 s File: .../python3.8/site-packages/h5netcdf/core.py Function: flush at line 689

Line # Hits Time Per Hit % Time Line Contents

689 def flush(self): 690 1 3.0 3.0 0.0 if 'r' not in self._mode: 691 1 1148237.0 1148237.0 25.2 self._set_unassigned_dimension_ids() 692 1 462926.0 462926.0 10.2 self._create_dim_scales() 693 1 2941779.0 2941779.0 64.6 self._attach_dim_scales() 694 1 2.0 2.0 0.0 if not self._preexisting_file and self._write_ncproperties: 695 self.attrs._h5attrs['_NCProperties'] = _NC_PROPERTIES ```

@shoyer This problem persisted all of this time, but since I faced it again, I did a bit of digging. (it's strange noone else noticed it so far as it's pretty bad)

I've line-profiled this snippet for various number of datasets already written to file (xarray.backends.api.to_netcdf):

https://github.com/pydata/xarray/blob/8452120e52862df564a6e629d1ab5a7d392853b0/xarray/backends/api.py#L1075-L1094

| Number of datasets in file | dump_to_store() | store_open() | store.close() | | --- | --- | --- | --- | | 0 | 88% | 1% | 10% | | 50 | 18% | 2% | 80% | | 200 | 4% | 2% | 94% |

The above can be measured simply in a notebook via %lprun -f xarray.backends.api.to_netcdf test_func(). The writing was done in mode='a', with blosc:zstd compression. All datasets are written into different groups (i.e. by passing group=...).