issue_comments

Given that this error only is caused when Python is shutting down, which is exactly a case in which we do not need to clean up open file objects, maybe we can remove the __del__ instead?

Something like: ```python import atexit

@atexit.register def _remove_del_method(): # We don't need to close unclosed files at program exit, # and may not be able to do, because Python is cleaning up # imports. del CachingFileManager.del ```

(I have not tested this!)

I agree that this seems very hard to test!

Have you verfied that this fixes things at least on your machine?

If we can find cases where we know concurrent writes are unsafe, we can definitely start raising errors. Better to be safe than to suffer from silent data corruption!

For anyone following along, I released a small package for reading TensorStore data into Xarray: https://github.com/google/xarray-tensorstore

Can we delete the "Flexible indexes" meeting? It doesn't happen anymore.

Allowing for explicitly passing a function matching the einsum interface is certainly more flexible than a boolean or enum argument, so @TomNicholas's suggestion of einsum_func=np.einsum is the version I would suggest.

The overhead from optimizing contraction paths is probably very small relative to the overhead of Xarray in general, so I would support setting optimize=True by default in Xarray, and/or using opt-einsum automatically if it is installed. JAX always use opt-einsum (opt-einsum is actually a hard dependency) and I have never heard any complaints.

In the hypothetical invocation open_dataset(..., return_encoding=True), do you envision the returned encoding as being a separate returned object, or would it still be an attribute on the Dataset object?

My expectation was that this would be a separate object, e.g., dataset, encoding = xarray.open_dataset(..., return_encoding=True), where encoding is a dict providing the encoding on each variable, and which could be passed as the encoding argument into to_netcdf(). That said, I can see how keeping encoding as variable attributes could also be convenient.

"disable all encoding propagation by discarding encoding attributes once a Dataset has been modified" would be an intermediate step, on the route to removing encoding from Xarray's data model entirely entirely.

(As a side note, I would probably spell this as open_dataset_with_encoding rather than having a function with a variable return signature.)

@dschwoerer are you sure that you are actually calculating the same thing in both cases? What exactly do the values of slc[d] look like? I would test thing on smaller inputs to verify. My guess is that you are inadvertently calculating something different, recalling that Xarray's broadcasting rules differ slightly from NumPy's.

I agree, automatic dimension only ever really made sense for interactive usecases, where a user could see and fix the default names.

It's a little late to change the default now to raising an error instead, but maybe we could add a warning?

Is LazilyIndexedArray really a public API? I don't see it on the API docs page.

Personally I would not want to guarantee external stability/availability for this API in its current state.

I like open_zarr(...) because it's less typing than open_dataset(..., engine='zarr'). The automatic backend detection logic doesn't currently work for Zarr, and in every case it adds overhead, which could be significant in the case of remote storage backends like Zarr.

So personally I would rather go the other direction and add open_netcdf().

The inconsistency in the chunks argument is non-ideal, but that could be handled by a separate deprecation process.

I agree, datatree is an important data structure for Xarray. My preferred way to do this would be follow @rabernat's suggestion and to fork the code the existing repo into the Xarray main codebase.

My main concern is that we should carefully evaluate the datatree API to make sure we won't want to change it soon. Once we bring it into Xarray, there will be a higher expectation that the interface will remain stable.

If you care about memory usage, you should explicitly close files after you use them, e.g., by calling ds.close() or by using a context manager. Does that work for you?

I think it's OK to still require bottleneck for ffill and bfill:

1. There are no numerical concerns: these functions simply repeat numbers forward (or backwards).
2. There is no good alternative to using a loop, and writing the loop in NumPy would be probitively slow.

I'm actually trying to merge IndexedCoordinates with Coordinates but I'm stuck: the latter is abstract and I don't really see how I could refactor it together with DatasetCoordinates and DataArrayCoordinates

Coordinates is abstract because in the (current) Xarray data model, it doesn't actually store any data -- coordinates are stored in private attributes of the original Dataset (._variables and ._coord_names) or DataArray (._coords). So Coordinates needs to serve as a proxy for the data.

In the long term, I think we should refactor Dataset/DataArray to actually store data (coordinate variables, indexes and dimension sizes) on Coordinates, but that's a bigger refactor.

For now, it's worth noting that the current Coordinates class isn't actually exposed in Xarray's public API, just the DatasetCoordinates and DataArrayCoordinates classes (and not even their constructors). So the intermdiate step I would try is: 1. Rename the current Coordinates baseclass to AbstractCoordinates. 2. Rename your IndexedCoordinates class to Coordinates. Expose it in the public API. Make sure it can handle DatasetCoordinates and DataArrayCoordinates in the constructor. 3. Maybe: use some Python magic to make DatasetCorodinates/DataArrayCoordinates subclasses of the new Coordinates. Or maybe make them actual subclassses, overriding many of the methods (including the constructor).

Long term, do you think it would make sense to merge together Indexes, Coordinates and IndexedCoordinates? They are sort of all containers for the same thing.

what do you think about the approach proposed here? I'd like to check that with you before going further with tests, docs, etc.

Generally this looks great to me!

• How to avoid building any default index? It seems silly to add or use the indexes argument just for that purpose? We could address that later.

My suggestion would be:

• coords passed as a dict: create default indexes
• coords passed as IndexedCoordinates: do not create defaults

Alternatively to an IndexedCoordinates subclass I was wondering if we could reuse the Coordinates base class?

Yes, this makes more sense to me!

What if the Indexes class was a facade based on IndexedCoordinates instead of the other way around?

Yes, I also agree! This makes more sense.

I also like the idea of creating specific Grouper objects for different types of selection, e.g., UniqueGrouper (the default), BinGrouper, TimeResampleGrouper, etc.

IMO, it's not correctly implementing the rule as you phrased it. You said "still present", which isn't the case here since the coordinate wasn't present before.

Another way of describing the current behavior would be that xarray keeps around "every coordinate which could possibly still be valid," which is determined based upon dimension names.

The main challenge is that "Coordinate variables should not have their coordinates changed" doesn't really make sense in Xarray's data model. Only Dataset or DataArray objects have coordinates, which apply to the the entire Dataset/DataArray.

Let me give an example of why we might want to keep scalar coordinates around. Consider a Dataset where lat and lon need to be represented as 2D arrays, along x and y dimensions. If we index out a single lat/lon point, i.e., ds.isel(x=0, y=0) it would have scalar coordinates "x", "y", "lat" and "lon." If we now convert any of these to a DataArray, arguably all the coordinates are still valid.

The "robust" part is really just a modification to how the limits for color scales are chosen, i.e., ignoring the bottom and top 2% of the dtaa from the color scale. So it sounds like what you're hoping for is separate per-column or per-row color scaling?

This was an intentional design choice, back in the early days of Xarray.

The rule Xarray uses for choosing which coordinates to associate with a DataArray created from a Dataset or DataArray is "every coordinate whose dimensions are still present on the new DataArray." This includes scalar coordinates, which are always kept around (because their dimensions are always included).

What rule would you suggest instead? I agree that the behavior in this case "feels" wrong, but keep in mind that once time because a scalar coordinate, Xarray doesn't have any way of knowing that it used to have its own dimension.

The case where Bottleneck really makes a difference was for moving window statistics, where it uses a smarter algorithm than our current NumPy implementation, which creating a moving window view.

Otherwise, I agree, it probably isn't worth the trouble.

That said -- we could also switch to smarter NumPy based algorithms to implement most moving window calculations, e.g,. using np.nancumsum for moving window means.

Difference between empty and non-empty arrays comes from different logic used for empty arrays in Variable._getitem_with_mask.

The problem itself can be solved by fixing dtypes.maybe_promote to return fill_value=np.float32('nan') instead of fill_value=np.nan on dtype('float32') input.

Thanks for the excellent report!

I agree, this sounds like a good fix to me.

I think something like the following would work: Replace the return line of maybe_promote python return np.dtype(dtype), fill_value with python dtype = np.dtype(dtype) fill_value = dtype.type(fill_value) return dtype, fill_value

Use cases would be in any web service that would like to provide the final data values back to a user in JSON.

For what it's worth, I think your users will have a poor experience with encoded JSON data for very large arrays. It will be slow to compress and transfer this data.

In the long term, you would probably do better to transmit the data in some binary form (e.g., by calling tobytes() on the underlying np.ndarray objects, or by using Xarray's to_netcdf).

Thanks for report and the PR!

This really needs a "minimal complete verifiable" example (e.g., by creating and loading a Zarr array with random data) so others can verify your reported the performance gains: https://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports https://stackoverflow.com/help/minimal-reproducible-example

To be honest, this fix looks a little funny to me, because NumPy's own implementation of tolist() is so similar. I would love to understand what is going on.

If you can reproduce the issue only using NumPy, it could also make more sense to file this as a upstream bug report to NumPy. The NumPy maintainers are in a better position to debug tricky memory allocation issues involving NumPy.

Explicitly providing indexes is an advanced user feature.

Agreed. However, xr.Dataset(coords={"x": pandas_midx}) is something that presumably a lot of users rely on (it is used extensively in Xarray's tests) and that we should really deprecate IMO. If we don't provide a convenient alternative, I expect many of those users will complain.

I agree -- we should support this for backwards compatibility (even if we deprecate it).

it's easier to explicitly manipulate indexes in the form of a dict

While generally I also prefer handling plain dict objects over custom dict-like objects, here I don't see much reasons of manipulating Xarray index objects independently of their coordinate variables. Indexes allows keeping them tied together, and it is already returned by .xindexes.

EDIT -- For more context: initially an Indexes object was almost equivalent to a Frozen(obj._indexes). In #5692 I tried hard and struggled to keep dealing with separate dicts of indexes and indexed variables, but in the end it made things much easier to encapsulate the variables in Indexes, which is also used internally in different places.

OK, this totally makes sense.

I don't love that it is possible to express invalid states in Xarray's data model. This motivated the creation of assert_internal_invariants and currently mostly is a concern for Xarray's own developers, but when we exposes the indexes argument, it will be easier for users to make the same sort of errors.

I wonder if we should consider the broader refactor of merging the Indexes and Coordinates objects, and expose the constructor as a public API. For clarity, I'll call it CoordinatesAndIndexes for now, but it could likely reuse the public name of Coordinates.

This would have a number of benefits:

1. It's impossible to provide inconsistent coords and indexes, because there is no separate indexes argument.
2. Likewise, it is impossible to create inconsistent coordinates and indexes on an existing Xarray object.
3. All the logic for verifying consistent coords and indexes can go in one place, shared between Dataset/DataArray. (Yes, it would be annoying to refactor Dataset to merge in variables from CoordinatesAndIndexes rather than the current separate Dataset._variables)
4. The public API also becomes clearer: if users want default indexes, they can pass a dict of variables into coords. If they want to copy indexes from another object, they can pass in a CoordinatesAndIndexes object (either from another Xarray object or constructed directly).

I no longer remember why I added these checks, but I certainly did not expect to see this sort of performance penalty!

I'm thinking of only accepting one or more instances of Indexes as indexes argument in the Dataset and DataArray constructors

I would lean against this, only because it's easier to explicitly manipulate indexes in the form of a dict than an xarray.Indexes object.

Explicitly providing indexes is an advanced user feature. I think it's OK to require users to do a bit more work in this case and to not necessarily do consistency checks (beyond verifying that the coordinate variables exist).

This is correct -- CFDatetimeCoder.encode is not lazy, even if the inputs are Dask arrays.

We would welcome contributions to fix this. This would entail making the encode look similar to the decode method (using lazy_elemwise_func).

We would also need a fall-back method for determining appropriate time units without looking at the array values. Something like seconds since 1900-01-01T00:00:00 would probably be a reasonable choice.

Cubed should define a concatenate function, so that should be OK

I think we could fix this by marking CachingFileManager with typing.final

Anyone want to review here? I think this should be ready to go in.

~~Actually maybe we don't need to keep files open after pickling... let me give this one more try.~~

Nevermind, this didn't work -- it still results in failing tests with dask-distributed on Windows.

OK, after a bit more futzing tests are passing and I think this is actually ready to go in. I ended up leaving in the reference counting after all -- I couldn't figure out another way to keep files open after a pickle round-trip.

Yes yes -- the sooner we can get rid of MultiIndex special cases the better!

I added @cjauvin's integration test, and verified that the fix works for the scipy and h5netcdf backends.

Unfortunately, it doesn't work yet for the netCDF4 backend. I don't think we can solve this in Xarray without fixes netCDF4-Python or the netCDF-C library: https://github.com/Unidata/netcdf4-python/issues/1195

As a concrete example, suppose we have two datasets: 1. Hourly predictions for 10 days 2. Daily observations for a month.

```python import numpy as np import pandas as pd import xarray

predictions = xarray.DataArray( np.random.RandomState(0).randn(24*10), {'time': pd.date_range('2022-01-01', '2022-01-11', freq='1h', closed='left')}, ) observations = xarray.DataArray( np.random.RandomState(1).randn(31), {'time': pd.date_range('2022-01-01', '2022-01-31', freq='24h')}, ) ```

Today, if you compare these datasets, they automatically align: ```

predictions - observations <xarray.DataArray (time: 10)> array([ 0.13970698, 2.88151104, -1.0857261 , 2.21236931, -0.85490761, 2.67796423, 0.63833301, 1.94923669, -0.35832191, 0.23234996]) Coordinates: * time (time) datetime64[ns] 2022-01-01 2022-01-02 ... 2022-01-10 ```

With this proposed change, you would get an error, e.g., something like: ```

predictions - observations ValueError: xarray objects are not aligned along dimension 'time':
array(['2022-01-01T00:00:00.000000000', '2022-01-02T00:00:00.000000000', '2022-01-03T00:00:00.000000000', '2022-01-04T00:00:00.000000000', '2022-01-05T00:00:00.000000000', '2022-01-06T00:00:00.000000000', '2022-01-07T00:00:00.000000000', '2022-01-08T00:00:00.000000000', '2022-01-09T00:00:00.000000000', '2022-01-10T00:00:00.000000000', '2022-01-11T00:00:00.000000000', '2022-01-12T00:00:00.000000000', '2022-01-13T00:00:00.000000000', '2022-01-14T00:00:00.000000000', '2022-01-15T00:00:00.000000000', '2022-01-16T00:00:00.000000000', '2022-01-17T00:00:00.000000000', '2022-01-18T00:00:00.000000000', '2022-01-19T00:00:00.000000000', '2022-01-20T00:00:00.000000000', '2022-01-21T00:00:00.000000000', '2022-01-22T00:00:00.000000000', '2022-01-23T00:00:00.000000000', '2022-01-24T00:00:00.000000000', '2022-01-25T00:00:00.000000000', '2022-01-26T00:00:00.000000000', '2022-01-27T00:00:00.000000000', '2022-01-28T00:00:00.000000000', '2022-01-29T00:00:00.000000000', '2022-01-30T00:00:00.000000000', '2022-01-31T00:00:00.000000000'], dtype='datetime64[ns]') vs array(['2022-01-01T00:00:00.000000000', '2022-01-01T01:00:00.000000000', '2022-01-01T02:00:00.000000000', ..., '2022-01-10T21:00:00.000000000', '2022-01-10T22:00:00.000000000', '2022-01-10T23:00:00.000000000'], dtype='datetime64[ns]') ```

Instead, you would need to manually align these objects, e.g., with xarray.align, reindex_like() or interp_like(), e.g., ```

predictions, observations = xarray.align(predictions, observations) or observations = observations.reindex_like(predictions) or predictions = predictions.interp_like(observations) ```

To (partially) simulate the effect of this change on a codebase today, you could write xarray.set_options(arithmetic_join='exact') -- but presmably it would also make sense to change Xarray's other alignment code (e.g., in concat and merge).

IMO we could first align (hah) these choices to be the same:

the exact mode of automatic alignment (outer vs inner vs left join) depends on the specific operation.

The problem is that user expectations are actually rather different for different options:

• With data movement operations like xarray.merge, you expect to keep around all existing data -- so you want an outer join.
• With inplace operations that modify an existing Dataset, e.g., by adding new variables, you don't expect the existing coordinates to change -- so you want a left join.
• With computate based operations (like arithmatic), you don't have an expectation that all existing data is unmodified, so keeping around a bunch of NaN values felt very wasteful -- hence the inner join.

What do you think of making the default FloatIndex use a reasonable (hard to define!) rtol for comparisons?

This would definitely be a step forward! However, it's a tricky nut to crack. We would both need a heuristic for defining rtol (some fraction of coordinate spacing?) and a method for deciding what the resulting coordinates should be (use values from the first object?).

Even then, automatic alignment is often problematic, e.g., imagine cases where a coordinate is defined in separate units.

I like option (4). If a multi-coordinate index needs to care about order, it can implement that logic itself.

You might look into different multiprocessing modes: https://docs.python.org/3/library/multiprocessing.html#contexts-and-start-methods

It may also be that the NetCDF or HDF5 libraries were simply not written in a way that can support multi-processing. This would not surprise me.

BTW is there any advantage or difference in terms of cpu and memory consumption in opening the file only one or let it open by every process? I'm asking because I thought opening in every process was just plain stupid but it seems to perform exactly the same, so maybe I'm just creating a problem where there is none

I agree, maybe this isn't worth the trouble. I have not seen it done successfully before.

Will that work in the same way if I still use process_map, which uses concurrent.futures under the hood?

Yes it should, as long as you're using multi-processing under the covers.

If you do multi-threading, then you would want to use threading.Lock(). But I believe we already apply a thread lock by default.

Can you try explicitly passing in a multiprocessing lock into the open_dataset() constructor? Something like: python from multiprocessing import Lock ds = xarray.open_dataset(file, lock=Lock())

(We automatically select appropriate locks if using Dask, but I'm not sure how we would do that more generally...)

I am tempted to suggest that the right way to handle Awkward array is to treat "var" dimensions similar to NumPy's structured dtypes, with shape only handling non-variable dimensions. The uniform dimensions are the only ones for which Xarray's API is going to work properly out of the box, and Awkward array properly already has the right tools for working with ragged dimensions.

Either way, I would definitely encourage figuring out some actual use-cases before building this out :)

As I understand it, the main purpose here is to remove Xarray lazy indexing class.

Maybe call this get_duck_array(), just to be a little more descriptive?

the unpacked data should match the type of these attributes, which must both be of type float or both be of type double. An additional restriction in this case is that the variable containing the packed data must be of type byte, short or int. It is not advised to unpack an int into a float as there is a potential precision loss.

I find this is ambiguous. is float above referring to float16 or float32? Is double referring to float64?

Yes, I'm pretty sure "float" means single precision (np.float32), given that "double" certainly means double precision (no.float64).

If so, then they do recommend float64, as requested by the OP, because the test data is short and the scale_factor is float64 (a.k.a double?)

Yes, I believe so.

The broader discussion here is about CF compliance. I find the spec ambiguous and xarray non-compliant. So many tests rely on the existing behavior, that I am unsure how best to proceed to improve compliance. I worry it may be a major refactor, and possibly break things relying on the existing behavior. I'd like to discuss architecture. Should this be in a new issue, if this closes with PR #6851? Should there be a new keyword for cf_strict or something?

I think we can treat this a bug fix and just go forward with it. Yes, some people are going to be surprised, but I don't think it's distruptive enough that we need to go to a major effort to preserve backwards compatibility. It should already be straightforward to work around by setting decode_cf=False when opening a file and then explicitly calling xarray.decode_cf().

I agree, I think only setting a few indexes at a time would be normal. If we eventually need convenience methods for setting multiple indexes we can add those later.

This sounds great to me!

I don't think we need support for setting multiple indexes at once in a single method call. You can call set_xindex multiple times for that if needed.

I just toggled the "Require a pull request before merging" option

Sounds good to me. The challenge will be defining a parallel computing API that works across all these projects, with their slightly different models.

Reopening because my second example print(stacked.assign_coords(z=[1, 2, 3, 4])) is still broken with the same error message. It would be ideal to fix this before the release.

Maybe a separate project in xarray-contrib would make sense?

I would be reluctant to add this into Xarray proper if we need a new external dependency for reading XML files.

On Wed, Jul 6, 2022 at 2:37 PM David Huard @.***> wrote:

I've got a first draft that parses an NcML document and spits out an xarray.Dataset. It does not cover all the NcML syntax, but the essential elements are there.

It uses xsdata https://xsdata.readthedocs.io/en/latest/ to parse the XML, using a datamodel automatically generated from the NcML 2-2 schema. I've scrapped test files from the netcdf-java https://github.com/Unidata/netcdf-java repo to create a test suite.

Wondering what's the best place to host the code, tests and test data so others may give it a spin ?

— Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/issues/2697#issuecomment-1176775280, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJJFVW32WV5YKZZP7KFVBTVSX4BZANCNFSM4GRUVDBQ . You are receiving this because you were mentioned.Message ID: @.***>

The simplest option would probably be a custom Zarr store that raises an error if you try to look at array data. This could be implemented as a subclass of an existing Zarr store (e.g., the in memory store) that raises an error in __getitem__ is the filename of requests does not start with ..

We have some tests with InaccessibleVariableDataStore for this sort of thing, but I don't know immediately how to hook that into the Zarr backend.

Dataset.rename does both variables and dimensions. That seems useful in many cases. I think it also makes more sense than Dataset.drop does, given that variables and dimensions often use the same names -- whereas drop mixed up variable names and index values.

NumPy mostly uses axis instead of axes, which we could copy.

It would be helpful to understand if there are also other uses of dim/dims that are inconsistent. Which is the most common pattern?

Should we call it rename_vars or rename_coords?

The later might make more sense, but then it wouldn't mirror Dataset.

This might fit more naturally into interp() as a new method like "nearest-valid" rather than in sel().

The difference is that sel() only looks at indexes (and not the data) to select out a single value, whereas interp() can combine adjacent values in arbitrary ways.

error: "ndarray[Any, dtype[Any]]" has no attribute "rename"

Yes, it's worth discussing. I don't know if there will be a satisfying resolution, though.

Looking at this mur-sst dataset in particular, it stores time in chunks of size 5. That means fetching the 6443 time values requires 1288 separate HTTP requests -- no wonder it's so slow! If the time axis were instead stored in a single chunk of 51 KB, Xarray would only need 3 small size HTTP requests to load the lat, lon and time indexes, which would probably complete in a fraction of a second.

That said, I agree that this would be nice to have in general.

but another option (post explicit index refactor) might be an option for opening a dataset without creating indexes for 1D coordinates along dimensions.

It might indeed be worth considering this case too in #6392. Maybe indexes=None (default) to create default indexes for 1D coordinates and indexes={} (empty dictionary) to explicitly skip creating indexes?

+1 this syntax makes sense to me!

One issue with relying only on Array and Group as currently implemented in Zarr-Python is that we can create array nodes outside of any group subfolder. e.g. one can currently create an Array directly at path 'array1' and this would put the chunks under 'data/root/array1/', and metadata at 'meta/root/array1.array.json'. However, the root itself is not a Group. A group is basically a subfolder under root (e.g.' open_group with path = group1 creates '/meta/root/group1/' folder and 'meta/root/group1.group.json' metadata). There is no mechanism in the spec to open root directly as a Group!

is there an issue on the Zarr side where this is currently being discussed?

I opened up https://github.com/zarr-developers/zarr-python/issues/1039

Early versions of Xarray used to have lazy loading of data for indexes, but we removed this for the sake of simplicity. In principle we could restore lazy indexes, but another option (post explicit index refactor) might be an option for opening a dataset without creating indexes for 1D coordinates along dimensions.

Another way to solve this sort of challenges might be to load index data in parallel when using Dask. Right now I believe the data corresponding to indexes is always loaded eagerly, without using Dask.

All that said -- Do you have a specific example where this has been problematic? In my experience it has been pretty reasonable to use xarray.Dataset objects for schema-like templates, even with index data needing to be loaded eagerly. Possibly another Zarr chunking scheme for your index data could be more efficient?

We could raise an error asking the user to switch to swap_dims.

This seems like a good idea

In the long term, we like to decouple indexes from coordinate, and make something like the following work: ds.set_coords(['lon']).rename(x='lon').set_index('lon')

👍 this looks great to me!

For whatever reason, Windows seems to be much stricter about requiring file handles to be explicitly closed. So my guess is that this could be solved by using open_dataset() as a context manager.

I'm a little skeptical that it makes sense to add special case logic into Xarray in an attempt to keep NumPy's "OWNDATA" flag up to date. There are lots of places where we create views of data from existing arrays inside Xarray operations.

There are definitely cases where Xarray's internal operations do memory copies followed by views, which would also result in datasets with misleading "OWNDATA" flags if you look only at resulting datasets, e.g., DataArray.interp() which definitely does internal memory copies: ```

y = xarray.DataArray([1, 2, 3], dims='x', coords={'x': [0, 1, 2]}) y.interp(x=0.5).data.flags C_CONTIGUOUS : True F_CONTIGUOUS : True OWNDATA : False WRITEABLE : True ALIGNED : True WRITEBACKIFCOPY : False UPDATEIFCOPY : False ```

Overall, I just don't think this is a reliable way to trace memory allocation with NumPy. Maybe you could do better by also tracing back to source arrays with .base?

Still relevant!

I think this would need to get updated on the NumPy side. Ideally NumPy ufuncs would be typed to check for __array_ufunc__. Something like: ```python from typing import Protocol, TypeVar

class HasArrayUFunc(Protocol): def array_ufunc(ufunc, method, inputs, *kwargs): pass

ArrayOrHasArrayUFunc = TypeVar("ArrayOrHasArrayUFunc", ndarray, HasArrayUFunc)

def exp(x: ArrayOrHasArrayUFunc) -> ArrayOrHasArrayUFunc: ... ```

@gagebeni please open a new discussion for your issue: https://github.com/pydata/xarray/discussions

This is still relevant

I will take a look soon!

On Thu, Apr 14, 2022 at 6:23 PM Maximilian Roos @.***> wrote:

Hi @cisaacstern https://github.com/cisaacstern — thanks a lot and welcome to xarray!

This looks very coherent, as far as the context I have. Any thoughts from others who know the area better?

— Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/pull/6476#issuecomment-1099769304, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJJFVVACXN36KM4XREIBCLVFDAKVANCNFSM5TE2KQ4Q . You are receiving this because you were mentioned.Message ID: @.***>

Thanks @malmans2 !

One issue with relying only on Array and Group as currently implemented in Zarr-Python is that we can create array nodes outside of any group subfolder. e.g. one can currently create an Array directly at path 'array1' and this would put the chunks under 'data/root/array1/', and metadata at 'meta/root/array1.array.json'. However, the root itself is not a Group. A group is basically a subfolder under root (e.g.' open_group with path = group1 creates '/meta/root/group1/' folder and 'meta/root/group1.group.json' metadata). There is no mechanism in the spec to open root directly as a Group!

is there an issue on the Zarr side where this is currently being discussed?

• The v3 spec requires a path be specified when calling open_group or open_consolidated. This PR currently just sets a default group name of 'xarray' if one is not specified via the group kwarg to ZarrStore.open_group. I think that is convenient, but one could instead be stricter and raise an error in this case.

Does Zarr v3 have a notion of a "root" group? That feels like a more sensible default to me, both for Xarray and Zarr-Python

• If a string corresponding to a filesystem path or URL is used for store, then it is not possible to infer which version of the zarr spec is desired. In this case, the user must specify zarr_version to choose the zarr protocol version. The default of zarr_version=None will infer the version from a zarr BaseStore subclass when possible, otherwise defaulting to zarr_version=2 for backwards compatibility.

This sounds fine for now, but I am concerned that it will slow the adoption of Zarr v3. Eventually, we would presumably want to change the default to version 3, but this is difficult to do if it entirely breaks backwards compatibility.

My preference would be for the default behavior to try opening Zarr v2, and fall back to opening in v3 mode, even if this requires attempting to open a file from the store. This is similar to how Xarray handles other Zarr versioning issues (e.g., for consolidated metadata). Perhaps Zarr-Python could raise an informative error that we could catch if the Zarr version is incorrect, or even handle this behavior itself?

Could you also add brief updates to mention NCZarr support in the docstring for open_zarr and the user guide? In particular this paragraph should be updated:

Xarray can’t open just any zarr dataset, because xarray requires special metadata (attributes) describing the dataset dimensions and coordinates. At this time, xarray can only open zarr datasets that have been written by xarray. For implementation details, see Zarr Encoding Specification.

As it is currently it is also not possible to write a zarr which follows the GDAL ZARR driver conventions. Writing the _CRS attribute also results in a TypeError:

Can you elaborate? What API are you using to do the write: python, netcdf-c, or what?

This error message comes from Xarray and can be triggered by calling to_zarr(): https://github.com/pydata/xarray/blob/facafac359c39c3e940391a3829869b4a3df5d70/xarray/backends/api.py#L162

I don't think netCDF-C needs to be involved at all, which is why I suggested opening a separate issue.

@wankoelias could you kindly open a new issue for writing GDAL ZARR?

This is the second follow-up item in https://github.com/pydata/xarray/issues/6293

I think we could definitely experiment with relaxing this constraint now, although ideally we would continue to check off auditing all of the methods in that long list first.

I think this is probably fine without a deprecation cycle. This is a very easy fix for users.

@DennisHeimbigner I think it would be great to standardize NCZarr as a super-set of the "Xarray-Zarr" standard! I think Xarray should indeed be able to read such files. If you want to read a sub-group, you can read the sub-group in a separate call to xarray.open_zarr().

@rabernat I would not be opposed to adding support inside Xarray for reading NCZarr data, specifically to understand NCZarr's encoding of dimension names when using Zarr-Python. This wouldn't give 100% compatibility with NCZarr, but it would be very close (maybe just with incorrect dtypes for attributes) with a minimal amount of work. I don't think it would be a big deal to look for .nczvar files.

OK, in it goes! Big thanks to @benbovy for seeing this through :)

OK, I think we’re good to go here?

The data type restriction here seems to date back to the original PR adding support for appending. I turned up this comment that seems to summarize the motivation for this check: https://github.com/pydata/xarray/pull/2706#issuecomment-502481584

I think the original issue was that appending a fixed-width string could be a problem if the fixed-width does not match the width of the existing string dtype stored in Zarr.

This obviously doesn't apply in this case, because you are adding an entirely new variable. So I guess the check could be removed in that case.

Another challenge with changing the meaning of slice is handling partial slices, e.g., what does slice(500, None) mean? With a monotonic decreasing index, that would select values below 500, but ignoring underlying coordinate order it would presumably mean selecting values above 500.

I think the separate new API (e.g., xarray.Between or .sel_between()) is probably a better idea.

I would like to merge this PR very soon so it can get testing before the next release. If anyone has any remaining concerns, please speak up!

Well, even if we keep squeeze as an option, I think squeeze=False would be much more consistent default behavior :)

One complication with using sel() with slice objects is that you can do selection over non-monotonic indexes, merely based on matching bounds: ```

data = xarray.DataArray([1, 2, 3, 4, 5], dims=['x'], coords=[[5, 1, 4, 3, 2]]) data <xarray.DataArray (x: 5)> array([1, 2, 3, 4, 5]) Coordinates: * x (x) int64 5 1 4 3 2 data.sel(x=slice(1, 3))) <xarray.DataArray (x: 3)> array([2, 3, 4]) Coordinates: * x (x) int64 1 4 3 ```

If we change the semantics of slice in sel() to do filtering rather than be concerned about order (which does seem much less useful), we should probably deprecate the handling of non-monotonic ascending or descending indexes.

Alternatively, we could either do the dedicated indexing object like xarray.Between(lower, upper) or have a dedicated method for selecting between values, e.g., perhaps data.sel_between(x=(1, 3)) or data.sel_bounds(x=(1, 3)).

This is probably worth fixing if possible in a straightforward way. I don't think anyone is well served by matching the behavior of Python list indexing here -- it's a strange edge that case that indexing a list like x[5:0] returns an empty list.

Breaking changes will continue to be very rare, and whenever possible will be preceeded by deprecation or future warnings for multiple months.

Adding a join argument sounds good to me. I do not remember why the default is an outer join.

One thing that came up in our discussion about this in the developer meeting today is that we could also pretty easily expose a "low level" API for IO using dictionaries of xarray.Variable objects. This intermediate representation could be useful for cleaning up data into a form suitable for conversion into Dataset objects.

On Wed, Feb 16, 2022 at 11:39 PM Alessandro Amici @.***> wrote:

@TomNicholas https://github.com/TomNicholas (cc @mraspaud https://github.com/mraspaud)

Do you have use cases which one of these designs could handle but the other couldn't?

The two main classes of on-disk formats that, I know of, which cannot be always represented in the "group is a Dataset" approach are:

• in netCDF following the CF conventions for groups https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#groups, it is legal for an array to refer to a dimension or a coordinate in a different group and so arrays in the same group may have dimensions with the same name, but different size / coordinate values,
• the current spec for the Next-generation file formats (NGFF) https://ngff.openmicroscopy.org for bio-imaging has all scales of the same 5D data in the same group.

I don't have an example at hand, but my impression is that satellite products that use HDF5 file format also place arrays with inconsistent dimensions / coordinates in the same group.

— Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/issues/4118#issuecomment-1042656377, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJJFVT27QD4RQDYZ2N4W7TU3SQ3BANCNFSM4NQEIKFQ . You are receiving this because you were mentioned.Message ID: @.***>

For what it's wroth, the recommended way to do this is to explicitly close the Dataset with ds.close() rather than using del ds.

Or with a context manager, e.g., python for num in range(100): with xr.open_dataset('data.{}.nc'.format(num)) as ds: # do some stuff, but NOT assigning any data in ds to new variables ...

The reason why this isn't allowed is because it's ambiguous what to do with the other variables that are not restricted to the region (['cell', 'face', 'layer', 'max_cell_node', 'max_face_nodes', 'node', 'siglay'] in this case).

I can imagine quite a few different ways this behavior could be implemented:

1. Ignore these variables entirely.
2. Ignore variables if they also already exist, but write new ones.
3. Write or overwrite both new and existing these variables.
4. Write new variables. Ignore existing variables only if they already exist with the same values, and if not, raise an error.

I believe your proposal here (removing these checks from _validate_region) would achieve (3), but I'm not sure that's the best option.

(4) seems like perhaps the most user-friendly option, but checking existing variables can add significant overhead. When experimenting adding region support Xarray-Beam, I found many cases where it was easy to inadvertently make large parallel pipelines much slower by downloaded existing variables.

The current solution is not to do any of these, and to force the user to make an explicit choice by dropping new variables, or write them in a separate call to to_zarr. I think it would also be OK to let a user explicitly opt-in to one of these behaviors, but I don't think guessing what the user wants would be ideal.

In the long term (cc @benbovy) I think we would ideally split IndexVariable into two classes:

1. FrozenVariable which is just an immutable Variable, and thus that can be safely used for coordinates that have indexes.
2. PandasIndexArray which wraps pandas.Index objects to satisfy the np.ndarray interface. This is the object which could allow duck_array_ops.isin to use the pandas.Index.isin method.

Oh, I guess the challenge is that apply_ufunc operates on arrays, not indexes. I'm not entirely sure how to deal with this easily....

Yes, I think replacing this with something like lambda x, y: x.isin(y) if isinstance(x, pd.Index) else np.isin(x, y) could work

Have you seen xarray.save_mfdataset?

In principle, it was designed for exactly this sort of thing.

Let me ponder the linked issue. This was not an intentional feature for compute=False, so I'd like to be sure we can be committed to supporting it before we document it :)

E.g., I think skipping this line would save some of the users in my original post a lot of time.

I don't think that line adds any measurable overhead. It's just telling dask to delay computation of a single function.

For sure this would be worth elaborating on in the Xarray docs! I wrote a little bit about this in the docs for Xarray-Beam: see "One recommended pattern" in https://xarray-beam.readthedocs.io/en/latest/read-write.html#writing-data-to-zarr