Operations on DataTree objects should not create duplicate coordinates on sub-trees

[shoyer]

What is your issue?

This is very obviously an issue using the repr from #9470:

In [20]: tree = DataTree(Dataset(coords={'x': [1, 2, 3]}))

In [21]: tree['/foo'] = DataTree(Dataset({'bar': ('x', [4, 5, 6])}))

In [22]: tree
Out[22]:
<xarray.DataTree>
Group: /
│   Dimensions:  (x: 3)
│   Coordinates:
│     * x        (x) int64 24B 1 2 3
└── Group: /foo
        Dimensions:  (x: 3)
        Data variables:
            bar      (x) int64 24B 4 5 6

In [23]: tree * 2
Out[23]:
<xarray.DataTree>
Group: /
│   Dimensions:  (x: 3)
│   Coordinates:
│     * x        (x) int64 24B 1 2 3
└── Group: /foo
        Dimensions:  (x: 3)
        Coordinates:
          * x        (x) int64 24B 1 2 3
        Data variables:
            bar      (x) int64 24B 8 10 12

Instead, the result for tree * 2 should look like:

<xarray.DataTree>
Group: /
│   Dimensions:  (x: 3)
│   Coordinates:
│     * x        (x) int64 24B 1 2 3
└── Group: /foo
        Dimensions:  (x: 3)
        Data variables:
            bar      (x) int64 24B 8 10 12

We probably also need to also revisit other Dataset methods that are ported to DataTree via map_over_subtree (xref #9472). Some of these (e.g., arithmetic, aggregations) can likely easily be corrected simply by mapping over nodes with inherited=False. Others (e.g., indexing) will need more careful consideration.

[shoyer]

As a starting point, it might make sense to disable most automatically ported DataTree methods, at least those that can't be implemented via map_over_subtree with inherited=False.

[TomNicholas]

I'm kind of surprised that the existing map_over_subtree tests didn't catch this.

[shoyer]

I'm kind of surprised that the existing map_over_subtree tests didn't catch this.

This could be because assert_identical (prior to #9473) doesn't currently check this.

[TomNicholas]

In the meeting just now we decided that this is a deficiency in the data model, not just in the implementation of map_over_subtree. We should change the model so that identical coordinates in children are automatically removed, so that the above would be automatically de-duplicated at construction time.

This is important to do as otherwise almost every user will run into this problem, especially using this pattern:

dt['lat'] = ...  # coord present on root

ds = dt[path].dataset
result_ds = some_operation(ds)
dt[path].dataset = result_ds  # coord now present on both root and descendant!

Currently this would duplicate the 'lat' coordinate between dt and dt[path], which would lead to very confusing behaviour, e.g.

del dt[path]['lat']
'lat' in dt[path]  # would return true!

(because only the duplicated coordinate on the child would have been deleted, not the original on the root, which would still be inherited).

IO

Disallowing duplicated coordinates in general would subtly affect round-tripping behaviour. It would no longer be possible to write out a netCDF/Zarr file with duplicated coordinates using DataTree, because DataTree could not represent them in the first place.

Round-tripping is then broken for any file with duplicated coordinates, as reading the file would de-duplicate the coordinates, and the information wouldn't be propagated in order to re-duplicate them before writing.

Mostly this wouldn't matter often - there is not much reason to define duplicated coordinates in netCDF/Zarr in the first place.

It might lead to subtle trickiness for files containing CF metadata that refers to other variables in different groups using absolute paths, but that is not something that xarray promises to handle anyway, as interpreting metadata contents is not part of xarray's data model (beyond what happens during decoding).

However again open_groups could still open such files without removing duplicates, and one could imagine writing out group-by-group to a single output file, and that should roundtrip perfectly.

All of this is a consequence of the fact that whilst very similar, the data models of netCDF, Zarr, and DataTree are not identical.

cc @owenlittlejohns @flamingbear @eni-awowale

---

Not sure exactly what the internal implementation of de-deduplication would be (perhaps some kind of NonDuplicatedChainMap?), but @shoyer said he is on it 🙂 A small change to map_over_subtree might also still be required to close this issue.

[shoyer]

It should be straightforward enough to handle duplicated index coordinates. Non-identical coordinates with an index are already disallowed between parent and child nodes, because the parent and child nodes would fail the alignment check in _pre_attach.

However, there is one tricky case: duplicated coordinates without an associated index. Here's an example of how these currently manifest themselves:

In [33]: tree = DataTree(Dataset(coords={'x': 0}), children={'child': DataTree(Dataset({'y': 1}))})

In [34]: tree
Out[34]:
<xarray.DataTree>
Group: /
│   Dimensions:  ()
│   Coordinates:
│       x        int64 8B 0
└── Group: /child
        Dimensions:  ()
        Data variables:
            y        int64 8B 1

In [35]: tree.map_over_subtree(lambda ds: ds)
Out[35]:
<xarray.DataTree>
Group: /
│   Dimensions:  ()
│   Coordinates:
│       x        int64 8B 0
└── Group: /child
        Dimensions:  ()
        Coordinates:
            x        int64 8B 0
        Data variables:
            y        int64 8B 1

Notice the duplicated coordinates /x and /child/x.

For consistency with the rest of Xarray, we should try to handle this like duplicate coordinates that result from other Xarray operations like arithmetic. As noted in #9481, the current strategy requires evaluating coordinate values to see if they are equal or not in order to decide where they appear. This is attractive for DataTree as well, but unfortunately is not compatible with lazy evaluation.

Instead, I propose that we switch to the equivalent of compat='override, similar to the proposal for Dataset/DataArray coordinates proposed in #9481. On DataTree, if an operation results in a node with a (non-index) coordinate with the same name as a parent coordinate, the child coordinate should be silently dropped.

An exception to this rule are operations that explicitly assign a child coordinate, e.g., DataTree.coords.__setitem__. Here we should raise an error (noting the existence of the parent coordinate) rather than being a silent no-op.

[eni-awowale]

Thanks for this great write up! I think in the long run this will probably help to avoid confusion and is probably the best way forward. But from more of a DAAC archiving, service provider and data validation perspective (lol I wear a lot of hats) this makes me a little nervous. At GES DISC we use DataTree to do data validation checks against server hosted services and cloud hosted services. The part that makes me nervous is that using DataTree to open netCDF4 files would be modifying the orginal file by removing the duplicate coordinate variables at each node. I don't think this will actually break anything in the backend for us but there is some complexity especialy with metadata that was mentioned.

Would folks be opposed to some kind of flag that's like inherit=True or de_duplicate=True? I think we could set that to being the default behavior and give the user an option to turn it off.

[TomNicholas]

Would folks be opposed to some kind of flag that's like inherit=True or de_duplicate=True? I think we could set that to being the default behavior and give the user an option to turn it off.

The proposed change here would make such a flag impossible - your intended result of de_duplicate=False would now be forbidden by the DataTree data model.

What we could perhaps do though is add a flag to constructors / openers that would raise on encountering duplicated coordinates instead of silently de-duplicating them. e.g.

dt = open_datatree(unvalidated.nc, duplicate_coords='raise')

where the error message tells you which coordinates are duplicated, and refers you to open_groups. The default could still be 'ignore' or maybe 'warn'. There's precedent for treating potential errors this way, see for example

xarray/xarray/namedarray/core.py

Line 1001 in aeaa082

missing_dims: ErrorOptionsWithWarn = "raise",

[shoyer]

What we could perhaps do though is add a flag to constructors / openers that would raise on encountering duplicated coordinates instead of silently de-duplicating them. e.g.

Yes, this is a great idea!

I would actually suggest such a stricter mode for open_datatree by default. Users who write data from Xarray will not encounter the issue, because Xarray will not be able to represent such datasets. So it's only something that will effect loading data from an external data provider, and it seems better to force users to make an intentional choice in such cases.

The cases where we should be more lenient by default (automatically dropping conflicts) are situations where users are constructing a new DataTree from a collection of Dataset objects, which were likely created using inherited=True.

[TomNicholas]

If we add the flag kwarg to DataTree.from_dict() then it can just be passed down from open_datatree immediately, and we can also explicitly write DataTree.from_dict(..., duplicate_coords='ignore') inside map_over_subtree to make our intent clear internally 🙂

[TomNicholas]

So summarising another meeting's worth of discussion on this... (including special guest @castelao)

The de-duplication idea has some issues.

1. map_over_subtree vs DataTree.__setitem__

As implied above, any de-duplication should ideally occur in the data model itself, not as a special feature of map_over_subtree. That means that the de-duplication also has to be applied for DataTree.__setitem__.

2. De-duplication vs Overriding

Comparing two coordinate variables to decide if they are de-duplicated could just be done by comparing names, as in #9510 (comment). But removing anything of the same name (and doing so in __setitem__) basically completely forbids "overriding" inherited coordinates with a new coordinate of the same name but different values. This would be a valid data model for DataTree, but an even more restrictive one than what we currently have (the ChainMap implementation we currently have allows for overriding any coordinate).

3. De-duplicating lazy variables by id

If you don't compare by names you have to compare something else. You could compare the id of the underlying Variable, i.e. id(Variable), but this will not recognize copies as being duplicates, even if their values are the same. e.g.

def mean(ds: Dataset) -> Dataset:
    return ds.mean()  # this creates a shallow copy of all the variables!

dt.map_over_subtree(mean)  # this would therefore still end up with duplicated coordinates

4. Only inherit index variables?

You can instead identify duplicated coordinates by comparing values directly, but this implies loading the variable into memory. Mostly our coordinate variables will be backed by in-memory indexes, so we can compare those and all would be fine. So a data model of "only inherit index-backed coordinates" works quite nicely, in that you can then always cheaply do comparison of inherited coordinates to check for de-duplication, and you can distinguish overridden from duplicated coordinates.

5. Accessing non-indexed inherited variables?

The fly in that ointment is that it's possible to have coordinate variables that are not backed by indexes, and as these can still be multi-dimensional they can still lazily point to large amounts of data.

If the data model is now "only inherit coordinate variables backed by indexes", then it limits the usefulness of inheritance.

If we try to do inheritance of non-indexed coordinate variables, we can't use the solution from step (4).

[kmuehlbauer]

My only concern is what @eni-awowale already mentioned, removing duplicated coordinates from child-groups will break the original file data model. Overriding coordinates in child groups is also a CF conventions feature, it should at least be handled correctly via some switch/flag.

[shoyer]

I made a few slides going through my proposed fix (only inheriting index variables):
https://docs.google.com/presentation/d/1zBjEsihBhK_U972jxHwaAZBbzS1-hd3aDLnO9uu2Ob4/edit?usp=sharing

[TomNicholas]

I realised that the .sel example in #9501 is interesting because you can navigate around the unintuitive behaviour by indexing in such a way that preserves the IndexVariable:

Click to see example

In [6]: ds = xr.tutorial.open_dataset("air_temperature").drop_attrs()

In [7]: ds_daily = ds.resample(time="D").mean("time")

In [8]: ds_weekly = ds.resample(time="W").mean("time")

In [9]: ds_monthly = ds.resample(time="ME").mean("time")

In [12]: dt = xr.DataTree.from_dict(
  ...:     {
  ...:         "/": ds.drop_dims("time"),
  ...:         "daily": ds_daily.drop_vars(["lat", "lon"]),
  ...:         "weekly": ds_weekly.drop_vars(["lat", "lon"]),
  ...:         "monthly": ds_monthly.drop_vars(["lat", "lon"]),
  ...:     }
  ...: )

In [13]: dt
Out[13]: 
<xarray.DataTree>
Group: /
│   Dimensions:  (lat: 25, lon: 53)
│   Coordinates:
│     * lat      (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0
│     * lon      (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0
├── Group: /daily
│       Dimensions:  (time: 730, lat: 25, lon: 53)
│       Coordinates:
│         * time     (time) datetime64[ns] 6kB 2013-01-01 2013-01-02 ... 2014-12-31
│       Data variables:
│           air      (time, lat, lon) float64 8MB 241.9 242.3 242.7 ... 295.9 295.5
├── Group: /weekly
│       Dimensions:  (time: 105, lat: 25, lon: 53)
│       Coordinates:
│         * time     (time) datetime64[ns] 840B 2013-01-06 2013-01-13 ... 2015-01-04
│       Data variables:
│           air      (time, lat, lon) float64 1MB 245.3 245.2 245.0 ... 296.6 296.2
└── Group: /monthly
      Dimensions:  (time: 24, lat: 25, lon: 53)
      Coordinates:
        * time     (time) datetime64[ns] 192B 2013-01-31 2013-02-28 ... 2014-12-31
      Data variables:
          air      (time, lat, lon) float64 254kB 244.5 244.7 244.7 ... 297.7 297.7

In [14]: dt.sel(lat=75, lon=300)
Out[14]: 
<xarray.DataTree>
Group: /
│   Dimensions:  ()
│   Coordinates:
│       lat      float32 4B 75.0
│       lon      float32 4B 300.0
├── Group: /daily
│       Dimensions:  (time: 730)
│       Coordinates:
│           lat      float32 4B 75.0
│           lon      float32 4B 300.0
│         * time     (time) datetime64[ns] 6kB 2013-01-01 2013-01-02 ... 2014-12-31
│       Data variables:
│           air      (time) float64 6kB 242.7 245.6 244.9 249.8 ... 254.8 255.6 256.8
├── Group: /weekly
│       Dimensions:  (time: 105)
│       Coordinates:
│           lat      float32 4B 75.0
│           lon      float32 4B 300.0
│         * time     (time) datetime64[ns] 840B 2013-01-06 2013-01-13 ... 2015-01-04
│       Data variables:
│           air      (time) float64 840B 247.2 251.7 256.2 261.4 ... 249.8 248.2 255.7
└── Group: /monthly
      Dimensions:  (time: 24)
      Coordinates:
          lat      float32 4B 75.0
          lon      float32 4B 300.0
        * time     (time) datetime64[ns] 192B 2013-01-31 2013-02-28 ... 2014-12-31
      Data variables:
          air      (time) float64 192B 254.0 252.8 256.9 258.7 ... 265.1 261.8 251.7


In [15]: dt.sel(lat=[75], lon=[300])
Out[15]: 
<xarray.DataTree>
Group: /
│   Dimensions:  (lat: 1, lon: 1)
│   Coordinates:
│     * lat      (lat) float32 4B 75.0
│     * lon      (lon) float32 4B 300.0
├── Group: /daily
│       Dimensions:  (time: 730, lat: 1, lon: 1)
│       Coordinates:
│         * time     (time) datetime64[ns] 6kB 2013-01-01 2013-01-02 ... 2014-12-31
│       Data variables:
│           air      (time, lat, lon) float64 6kB 242.7 245.6 244.9 ... 255.6 256.8
├── Group: /weekly
│       Dimensions:  (time: 105, lat: 1, lon: 1)
│       Coordinates:
│         * time     (time) datetime64[ns] 840B 2013-01-06 2013-01-13 ... 2015-01-04
│       Data variables:
│           air      (time, lat, lon) float64 840B 247.2 251.7 256.2 ... 248.2 255.7
└── Group: /monthly
      Dimensions:  (time: 24, lat: 1, lon: 1)
      Coordinates:
        * time     (time) datetime64[ns] 192B 2013-01-31 2013-02-28 ... 2014-12-31
      Data variables:
          air      (time, lat, lon) float64 192B 254.0 252.8 256.9 ... 261.8 251.7

The first call dt.sel(lat=75, lon=300) behaves counter-intuitively in that it duplicates the lat and lon coords. However dt.sel(lat=[75], lon=[300]) keeps the coordinate inheritance in the way you would expect. This only works because dt.sel(lat=[75], lon=[300]) does not actually get rid of the Index underlying the lat/lon coordinate variables, so the result is still subject to the inheritance rules.

In some ways this is weird, but it also kind of makes sense: dimensions are inherited, and the dt.sel(lat=75, lon=300) call removes the dimensions (by making lat/lon into scalars). dt.sel(lat=[75], lon=[300]) keeps the dimensions by indexing to create 1D arrays instead of scalars.

I wondered what people thought about encouraging these kind of index-preserving operations as an escape hatch, especially @castelao.

EDIT: This is really the same thing as what @kmuehlbauer and @castelao were saying last meeting - that knowing that they could add length-1 dimensions to scalars to have them be inherited was a least a way to have "scalars" be inherited...

[shoyer]

For sel() and isel() in particular, I think we should have special case logic to keep the indexed coordinates at the same level where they were defined originally.

[TomNicholas]

I'm not sure I totally follow. You mean so that even if the result of indexing is a non-indexed coordinate (e.g. a scalar), special-case logic in isel/sel de-deuplicates that result coordinate automatically.

[shoyer]

I'm not sure I totally follow. You mean so that even if the result of indexing is a non-indexed coordinate (e.g. a scalar), special-case logic in isel/sel de-deuplicates that result coordinate automatically.

Yes, exactly. Coordinates defined at a higher level could be automatically be excluded from the indexing result.

This is similar to how Dataset.sel() is not exactly equivalent to mapping over DataArray.sel() -- dimensions that do not apply to a particular array as dropped from the indexers. (Incidentally, we should add this feature for DataTree, too!)

[TomNicholas]

(Incidentally, we should add this feature for DataTree, too!)

I think that's already tracked in #8949. This would be a nice thing to get in before release, because it's definitely going to annoy people who try to use datatrees with inheritance for analysis.

[TomNicholas]

This should be closed by #9555