What happened?

The following code yields a segfault on my machine (and many other machines with a similar environment)

``` import xarray filename = 'tiny.nc.txt' engine = "netcdf4"

dataset = xarray.open_dataset(filename, engine=engine)

i = 0 for i in range(60): xarray.open_dataset(filename, engine=engine) ```

tiny.nc.txt mrc.nc.txt

What did you expect to happen?

Not to segfault.

Minimal Complete Verifiable Example

1. Generate some netcdf4 with my application.
2. Trim the netcdf4 file down (load it, and drop all the vars I can while still reproducing this bug)
3. Try to read it.

```Python import xarray from tqdm import tqdm filename = 'mrc.nc.txt' engine = "h5netcdf" dataset = xarray.open_dataset(filename, engine=engine)

for i in tqdm(range(60), desc=f"filename={filename}, enine={engine}"): xarray.open_dataset(filename, engine=engine)

engine = "netcdf4"

dataset = xarray.open_dataset(filename, engine=engine) for i in tqdm(range(60), desc=f"filename={filename}, enine={engine}"): xarray.open_dataset(filename, engine=engine)

filename = 'tiny.nc.txt'

engine = "h5netcdf" dataset = xarray.open_dataset(filename, engine=engine) for i in tqdm(range(60), desc=f"filename={filename}, enine={engine}"): xarray.open_dataset(filename, engine=engine)

engine = "netcdf4"

dataset = xarray.open_dataset(filename, engine=engine) for i in tqdm(range(60), desc=f"filename={filename}, enine={engine}"): xarray.open_dataset(filename, engine=engine) ```

hand crafting the file from start to finish seems to not segfault: ``` import xarray import numpy as np engine = 'netcdf4'

dataset = xarray.Dataset()

coords = {} coords['image_x'] = np.arange(1, dtype='int') dataset = dataset.assign_coords(coords)

dataset['image'] = xarray.DataArray( np.zeros((1,), dtype='uint8'), dims=('image_x',) )

%%

dataset.to_netcdf('mrc.nc.txt')

%%

dataset = xarray.open_dataset('mrc.nc.txt', engine=engine)

for i in range(10): xarray.open_dataset('mrc.nc.txt', engine=engine)

```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.
• [X] Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

Python i=0 passes i=1 mostly segfaults, but sometimes it can take more than 1 iteration

Anything else we need to know?

At first I thought it was deep in hdf5, but I am less convinced now

xref: https://github.com/HDFGroup/hdf5/issues/3649

Environment

Is your feature request related to a problem?

I'm creating duckarrays for various file backed datastructures for mine that are naturally "chunked". i.e. different parts of the array may appear in completely different files.

Using these "chunks" and the "strides" algorithms can better decide on how to iterate in a convenient manner.

For example, an MP4 file's chunks may be defined as being delimited by I frames, while images stored in a TIFF may be delimited by a page.

So for me, chunks are not so useful for parallel computing, but more for computing locally and choosing the appropriate way to iterate through a large arrays (TB of uncompressed data).

Describe the solution you'd like

I think a default Chunk manager could simply implement compute as np.asarray as a default instance, and be a catchall to all other instances.

Advanced users could then go in an reimplement their own chunkmanager, but I was unable to use my duckarrays that incldued a chunk property because they weren't associated with any chunk manager.

Something as simple as:

```patch diff --git a/xarray/core/parallelcompat.py b/xarray/core/parallelcompat.py index c009ef48..bf500abb 100644 --- a/xarray/core/parallelcompat.py +++ b/xarray/core/parallelcompat.py @@ -681,3 +681,26 @@ class ChunkManagerEntrypoint(ABC, Generic[T_ChunkedArray]): cubed.store """ raise NotImplementedError() + + +class DefaultChunkManager(ChunkMangerEntrypoint): + def init(self) -> None: + self.array_cls = None + + def is_chunked_array(self, data: Any) -> bool: + return is_duck_array(data) and hasattr(data, "chunks") + + def chunks(self, data: T_ChunkedArray) -> T_NormalizedChunks: + return data.chunks + + def compute(self, data: T_ChunkedArray | Any, kwargs) -> tuple[np.ndarray, ...]: + raise tuple(np.asarray(d) for d in data) + + def normalize_chunks(self, args, kwargs): + raise NotImplementedError() + + def from_array(self, *args, kwargs): + raise NotImplementedError() + + def apply_gufunc(self, args, *kwargs): + raise NotImplementedError()

```

Describe alternatives you've considered

I created my own chunk manager, with my own chunk manager entry point.

Kinda tedious...

Additional context

It seems that this is related to: https://github.com/pydata/xarray/pull/7019

xref: https://github.com/pydata/xarray/issues/8733

What happened?

I'm not too sure this is a bug report, but I figured I would share some of the investigation I've done on the topic of writing large datasets to netcdf.

For clarity, the usecase I'm considering is writing large in-memory array to persistant storage on Linux.

• Array size 4-100GB
• File Format Netcdf
• Dask: No.
• Hardware: Some very modern SSD that can write more than 1GB/s (Sabrent Rocket 4 Plus for example).
• Operating System: Linux
• Xarray version: 0.21.1 (or around there)

The symptoms are two fold: 1. The write speed is slow. About 1GB/s, much less than the 2-3 GB/s you can get with other means. 2. The Linux disk cache just keeps filling up.

Its quite hard to get good performance from systems, so I"m going to put a few more constraints on the type of data we are are writing: 1. The underlying numpy array must be alight to the linux Page boundary of 4096 bytes. 2. The underlying numpy array must have been pre-faulted and not swapped. (Do not use np.zeros, it doesn't fault the memory)

I feel like these two options are rather easy to get to as I'll show in my example.

What did you expect to happen?

I want to be able to write at 3.2GB/s with my shiny new SSD.

I want to leave my RAM unused when I'm archiving to disk.

Minimal Complete Verifiable Example

```Python import numpy as np import xarray as xr

def empty_aligned(shape, dtype=np.float64, align=4096): if not isinstance(shape, tuple): shape = (shape,)

dtype = np.dtype(dtype)
size = dtype.itemsize
# Compute the final size of the array
for s in shape:
    size *= s

a = np.empty(size + (align - 1), dtype=np.uint8)
data_align = a.ctypes.data % align
offset = 0 if data_align == 0 else (align - data_align)
arr = a[offset:offset + size].view(dtype)
# Don't use reshape since reshape might copy the data.
# This is the suggested way to assign a new shape with guarantee
# That the data won't be copied.
arr.shape = shape
return arr

dataset = xr.DataArray( empty_aligned((4, 1024, 1024, 1024), dtype='uint8'), name='mydata').to_dataset()

Fault and write data to this dataset

dataset['mydata'].data[...] = 1

%time dataset.to_netcdf("test", engine='h5netcdf') %time dataset.to_netcdf("test", engine='netcdf4') ```

Relevant log output

Both output about 3.5s equivalent to just about 1GB/s.

To get to about 3 ish GB/s (taking about 1.27s to write a 4GB array). One needs to do a few things:

1. You must align the underlying data to disk.
2. h5netcdf (h5py) backend https://github.com/h5py/h5py/pull/2040
3. netcdf4: https://github.com/Unidata/netcdf-c/pull/2206
4. You must use a driver that bypasses the operating system cache
5. https://github.com/h5py/h5py/pull/2041

For the h5netcdf backend you would have to add the following kwargs to h5netcdf constructor kwargs = { "invalid_netcdf": invalid_netcdf, "phony_dims": phony_dims, "decode_vlen_strings": decode_vlen_strings, 'alignment_threshold': alignment_threshold, 'alignment_interval': alignment_interval, }

Anything else we need to know?

The main challenge is that while writing aligned data this way is REALLY fast, writing small chunks and unaligned data becomes REALLY slow.

Personally, I think that someone might be able to write a new HDF5 driver that does better optimization, I feel like this can help people loading large datasets which seems to be a large part of the community of xarray users.

Environment

``` INSTALLED VERSIONS

commit: None python: 3.9.9 (main, Dec 29 2021, 07:47:36) [GCC 9.4.0] python-bits: 64 OS: Linux OS-release: 5.13.0-30-generic machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: ('en_US', 'UTF-8') libhdf5: 1.12.1 libnetcdf: 4.8.1

xarray: 0.21.1 pandas: 1.4.0 numpy: 1.22.2 scipy: 1.8.0 netCDF4: 1.5.8 pydap: None h5netcdf: 0.13.1 h5py: 3.6.0.post1 Nio: None zarr: None cftime: 1.5.2 nc_time_axis: None PseudoNetCDF: None rasterio: None cfgrib: None iris: None bottleneck: None dask: 2022.01.1 distributed: None matplotlib: 3.5.1 cartopy: None seaborn: None numbagg: None fsspec: 2022.01.0 cupy: None pint: None sparse: None setuptools: 60.8.1 pip: 22.0.3 conda: None pytest: None IPython: 8.0.1 sphinx: None ```

h5py includes some additions of mine that allow you to use the DIRECT driver and I am using a version of HDF5 that is built with the DIRECT driver.

Is your feature request related to a problem?

calling np.empty_like seems to be instantiating the whole array.

```python from xarray.tests import InaccessibleArray import xarray as xr import numpy as np

array = InaccessibleArray(np.zeros((3, 3), dtype="uint8")) da = xr.DataArray(array, dims=["x", "y"])

np.empty_like(da) ```

python Traceback (most recent call last): File "/home/mark/t.py", line 8, in <module> np.empty_like(da) File "/home/mark/mambaforge/envs/dev/lib/python3.9/site-packages/xarray/core/common.py", line 165, in __array__ return np.asarray(self.values, dtype=dtype) File "/home/mark/mambaforge/envs/dev/lib/python3.9/site-packages/xarray/core/dataarray.py", line 732, in values return self.variable.values File "/home/mark/mambaforge/envs/dev/lib/python3.9/site-packages/xarray/core/variable.py", line 614, in values return _as_array_or_item(self._data) File "/home/mark/mambaforge/envs/dev/lib/python3.9/site-packages/xarray/core/variable.py", line 314, in _as_array_or_item data = np.asarray(data) File "/home/mark/mambaforge/envs/dev/lib/python3.9/site-packages/xarray/tests/__init__.py", line 151, in __array__ raise UnexpectedDataAccess("Tried accessing data") xarray.tests.UnexpectedDataAccess: Tried accessing data

Describe the solution you'd like

I'm not too sure. This is why I raised this as a "feature" and not a bug.

On one hand, it is pretty hard to "get" the underlying class.

Is it a:

• numpy array
• a lazy thing that looks like a numpy array?
• a dask array when it is dask?

I think that there are also some nuances between:

1. Loading an nc file from a file (where things might be handled by dask even though you don't want them to be)
2. Creating your xarray from in memory.

Describe alternatives you've considered

for now, i'm trying to avoid empty_like or zeros_like.

In general, we haven't seen much benefit from dask and cuda still needs careful memory management.

Additional context

No response

What happened?

When calling reset_coords on a dataset that is loaded from disk, the coordinates are not removed from the encoding of the variable.

This means, that at save time they will be resaved as coordinates... annoying. (and erroneous)

What did you expect to happen?

No response

Minimal Complete Verifiable Example

```Python import xarray as xr

dataset = xr.Dataset( data_vars={'images': (('y', 'x'), np.zeros((10, 2)))}, coords={'zar': 1} )

dataset.to_netcdf('foo.nc', mode='w')

%%

foo_loaded = xr.open_dataset('foo.nc')

foo_loaded_reset = foo_loaded.reset_coords()

%%

assert 'zar' in foo_loaded.coords assert 'zar' not in foo_loaded_reset.coords assert 'zar' in foo_loaded_reset.data_vars foo_loaded_reset.to_netcdf('bar.nc', mode='w')

%% Now load the dataset

bar_loaded = xr.open_dataset('bar.nc') assert 'zar' not in bar_loaded.coords, 'zar is erroneously a coordinate'

%%

This is the problem

assert 'zar' not in foo_loaded_reset.images.encoding['coordinates'].split(' '), "zar should not be in here" ```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

No response

Anything else we need to know?

for _, variable in obj._variables.items(): coords_in_encoding = set(variable.encoding.get('coordinates', ' ').split(' ')) variable.encoding['coordinates'] = ' '.join(coords_in_encoding - set(names)) suggested fix in dataset.py, reset_coords

https://github.com/pydata/xarray/blob/513ee34f16cc8f9250a72952e33bf9b4c95d33d1/xarray/core/dataset.py#L1734

Environment

Is your feature request related to a problem?

I'm trying to extend the features of xarray for a new backend I'm developing internally. The main use case that we are trying to open a multi 100's of GB dataset, slice out a smaller dataset (10s of GB) and write it.

However, when we try to use functions like prepare_variable, the way they are currently written, they implicitely instantiate the whole data, (potentially 10s of GB) which incurs a huge "time cost" at a surprising (to me) point in the code.

https://github.com/pydata/xarray/blob/6e77f5e8942206b3e0ab08c3621ade1499d8235b/xarray/backends/h5netcdf_.py#L338

Describe the solution you'd like

Would it be possible to just remove the second return value from prepare_variable? It isn't particuarly "useful" and easy to obtain from the inputs to the function.

Describe alternatives you've considered

I'm proably going to create a new method, with a not so well chosen name like prepare_variable_no_data that does the above, but only for my backend. My code path that needs this only uses our custom backend.

Additional context

I think this would be useful, in general for other users that need more out of memory computation. I've found that you really have to "buy into" dask, all the way to the end, if you want to see any benefits. As such, if somebody used a dask array, this would create a serial choke point in:

https://github.com/pydata/xarray/blob/6e77f5e8942206b3e0ab08c3621ade1499d8235b/xarray/backends/common.py#L308

What is your issue?

In https://github.com/pydata/xarray/pull/7221 I showed that a major contributor the slowdown in inserting a new element was the cost associated with an internal only debugging assert statement.

The benchmarks results 7221 and 7222 are pretty useful to look at.

Thank you for encouraging the creation of a "benchmark" so that we can monitor the performance of element insertion.

Unfortunately, that was the only "free" lunch I got.

A few other minor improvements can be obtained with: https://github.com/pydata/xarray/pull/7222

However, it seems to me that the fundamental reason this is "slow" is because element insertion is not so much "insertion" as it is: * Dataset Merge * Dataset Replacement of the internal methods.

This is really solidified in the https://github.com/pydata/xarray/blob/main/xarray/core/dataset.py#L4918

In my benchmarks, I found that in the limit of large datasets, list comprehensions of 1000 elements or more were often used to "search" for variables that were "indexed" https://github.com/pydata/xarray/blob/ca57e5cd984e626487636628b1d34dca85cc2e7c/xarray/core/merge.py#L267

I think a few speedsups can be obtained by avoiding these kinds of "searches" and list comprehensions. However, I think that the dataset would have to provide this kind of information to the merge_core routine, instead of the merge_core routine recreating it all the time.

Ultimately, I think you trade off "memory footprint" (due to the potential increase of datastructures you keep around) of a dataset, and "speed".

Anyway, I just wanted to share where I got.

Is your feature request related to a problem?

I briefly tried to see if any issue was brought up but couldn't.

I'm hoping to be able to pass additional keyword arguments to the engine when using to_netcdf. https://xarray.pydata.org/en/stable/generated/xarray.open_dataset.html

However, it doesn't seem to easy to do so.

Similar to how open_dataset has an additional **kwargs parameter, would it be reasonable to add a similar parameter, maybe engine_kwargs to the to_netcdf to allow users to pass additional parameters to the engine?

Describe the solution you'd like

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

dataset = xr.DataArray( data=np.zeros(3), name="hello" ).to_dataset()

dataset.to_netcdf("my_file.nc", engine="h5netcdf", engine_kwargs={"decode_vlen_strings=True"}) ```

Describe alternatives you've considered

One could forward the additional keyword arguments with **kwargs. I just feel like this makes things less "explicit".

Additional context

h5netcdf emits a warning that is hard to disable without passing a keyword argument to the constructor. https://github.com/h5netcdf/h5netcdf/issues/132

Also, for performance reasons, it might be very good to tune things like the storage data alignment.

See upstream: https://github.com/zarr-developers/zarr-python/issues/551

It seems that using a ZipStore creates 1 byte objects for Unicode string attributes.

For example, saving the same Dataset with a DirectoryStore and a Zip Store creates an attribute for a unicode array with 20 bytes in size in the first, and 1 byte in size in the second.

In fact, ubuntu file roller isn't even allowing me to extract the files.

I have a feeling it is due to the note in the zarr documentation

Note that Zip files do not provide any way to remove or replace existing entries.

https://zarr.readthedocs.io/en/stable/api/storage.html#zarr.storage.ZipStore

MCVE Code Sample

ZipStore python import xarray as xr import zarr x = xr.Dataset() x['hello'] = 'world' x with zarr.ZipStore('test_store.zip', mode='w') as store: x.to_zarr(store) with zarr.ZipStore('test_store.zip', mode='r') as store: x_read = xr.open_zarr(store).compute()

Working Directory Store example python import xarray as xr import zarr x = xr.Dataset() x['hello'] = 'world' x store = zarr.DirectoryStore('test_store2.zarr') x.to_zarr(store) x_read = xr.open_zarr(store) x_read.compute() assert x_read.hello == x.hello

Expected Output

The string metadata should work.

Output of xr.show_versions()