Repository: pydata/xarray

Issue:
Add option to choose the source of global attributes in mfdataset.
#### Code Sample

```python
import numpy as np
import xarray as xr

# prepare fake data
time = np.arange(1000)
data = time**2

# write to multiple netcdf files
for i in range(10):
    filename = 'ds{:d}.nc'.format(i)
    ds = xr.Dataset({'data': (['time'],  data[100*i:100*i+100])},
                    coords={'time': time[100*i:100*i+100]},
                    attrs={'history': 'Created ' + filename + '.'})
    ds.to_netcdf(filename)

# open as mfdataset
with xr.open_mfdataset('ds?.nc') as ds:
    print ds.history
```
#### Problem description

Currently, global attributes of multi-file datasets are taken from the first file in the list.

https://github.com/pydata/xarray/blob/b87b684b36cf5adbe4dca208aed0c69c44fc44c4/xarray/backends/api.py#L649

I think this is a problem in the context of consecutive model runs where history is appended in each subsequent run. When opening the results as mfdataset, history is taken from the first run file all subsequent history is lost.

NetCDF4 has a new keyword argument to set the `master_file` in a `MFDataset` (Unidata/netcdf4-python#835). Would it be possible to add a similar option in xarray?

#### Output of ``xr.show_versions()``

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 2.7.15.final.0
python-bits: 64
OS: Linux
OS-release: 4.14.65-1-MANJARO
machine: x86_64
processor: 
byteorder: little
LC_ALL: None
LANG: fr_CH.utf8
LOCALE: None.None

xarray: 0.10.8
pandas: 0.23.4
numpy: 1.15.0
scipy: 1.1.0
netCDF4: 1.4.2
h5netcdf: None
h5py: None
Nio: None
zarr: None
bottleneck: None
cyordereddict: None
dask: 0.18.2
distributed: None
matplotlib: 2.2.3
cartopy: 0.15.1
seaborn: None
setuptools: 40.0.0
pip: 18.0
conda: None
pytest: None
IPython: 5.8.0
sphinx: 1.7.6
</details>


Discussion:
juseg: I think I'm done. Can someone look at it? The `master_file` keyword argument is borrowed from NetCDF (see #2382 and Unidata/netcdf4-python#835), although `xarray`'s mechanism is independant.

The default is `0`, which is consistent with current `xarray` behaviour. When concatenating model output in time it would be more logical to use `-1`, i.e. preserve history of the last file.
dcherian: Thanks @juseg . I've left a few comments.

I see that this is your first PR. Welcome to xarray! and thanks for contributing :clap: 
juseg: This will add a new kwarg in `open_mfdataset`. The current name is `attrs_file`, defaulting to 0 (current behaviour). Or? Suggestions welcome.
TomNicholas: Thanks for this @juseg. The only problem I see is that a scalar number to specify the file only makes sense if it's a 1D list, but `open_mfdataset` can also accept nested list-of-lists (with `combine='nested'`), or ignore the order of the input entirely (with `combine='manual'`). What happens if you pass a list-of-lists of datasets? 

On the other hand specifying the particular filepath or object makes sense in all cases, so perhaps the easiest way to avoid ambiguity would be to restrict to that option? (The default would just be left as-is.)
juseg: @TomNicholas Thanks for bringing the discussion live again! I'm not sure what happens in those cases, but I'm confident the default behaviour is unchanged, i.e. the attributes file is 0, whatever that 0 means (see my first commit).

If this is an issue I would suggest to discuss in a separate thread, as I think it is independent from my changes. On the other hand I am eager to keep the file number option because (1) `attrs_file=-1` is the behaviour that I need (to ensure that history is always preserved) and (2) `attrs_file=0` is the current behaviour (again, whatever that means for nested lists).
TomNicholas: I'm not sure we should merge changes if we're unsure how they will behave in certain circumstances.

> On the other hand I am eager to keep the file number option because (1) attrs_file=-1 is the behaviour that I need

If we kept just the string specifier, you could still solve the problem of preserving the history:
```python
files_to_open = ['filepath1', 'filepath2']

ds = open_mfdataset(files_to_open, attrs_file=files_to_open[-1])
```
But then the option would always have clear and well-defined behaviour, even in more complex cases like `combine_by_coords`or `combine='nested'` with a >1D input file list.

juseg: @TomNicholas I've had a closer look at the code. Nested lists of file paths are processed by:

https://github.com/pydata/xarray/blob/f2b2f9f62ea0f1020262a7ff563bfe74258ffaa1/xarray/backends/api.py#L881-L882

Using the method defined in:

https://github.com/pydata/xarray/blob/f2b2f9f62ea0f1020262a7ff563bfe74258ffaa1/xarray/core/combine.py#L15-L46

In Python 3.7+ the list of `paths` is essentially flattened list e.g.:

```
>>> import xarray as xr
>>> paths = [list('abc'), list('def')]
>>> list(xr.core.combine._infer_concat_order_from_positions(paths).values())
['a', 'b', 'c', 'd', 'e', 'f']
```
Unfortunately the current code uses a dictionary which means that in Python 3.6- the order is not guaranteed preserved. This also implies that the current default `0` is not well defined in case of nested lists.

https://github.com/pydata/xarray/blob/f2b2f9f62ea0f1020262a7ff563bfe74258ffaa1/xarray/backends/api.py#L900
https://github.com/pydata/xarray/blob/f2b2f9f62ea0f1020262a7ff563bfe74258ffaa1/xarray/backends/api.py#L964

On the other hands the `ids` are essentially ND indexes that could perhaps be used...

```
>>> list(xr.core.combine._infer_concat_order_from_positions(paths).keys())
[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
```

Or should we just stick to file paths as you suggest? And leave the default as is (e.g. ambiguous for Python 3.6-)?
TomNicholas: Thanks @juseg .

> in Python 3.6- the order is not guaranteed preserved.

I think for python 3.6 and above the [order is preserved](https://stackoverflow.com/questions/39980323/are-dictionaries-ordered-in-python-3-6) isn't it?

> the current default 0 is not well defined in case of nested lists.

Yes, this is what I was thinking of.

> the ids are essentially ND indexes that could perhaps be used...

We *could* do this, and that's how we would solve it in general, but I don't really think it's worth the effort/complexity.

> Or should we just stick to file paths as you suggest? And leave the default as is

I think so - if we do this then users can still easily pick the attributes from the file of their choosing (solving the original issue), and if someone wants to be able to choose the `attrs_file` in another way later then we can worry about that then.
TomNicholas: I think there's nothing left to do here, thanks @juseg!

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 93df2b65..eacf8433 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -37,6 +37,9 @@ New Features
 - Added the ``count`` reduction method to both :py:class:`~core.rolling.DatasetCoarsen`
   and :py:class:`~core.rolling.DataArrayCoarsen` objects. (:pull:`3500`)
   By `Deepak Cherian <https://github.com/dcherian>`_
+- Add `attrs_file` option in :py:func:`~xarray.open_mfdataset` to choose the
+  source file for global attributes in a multi-file dataset (:issue:`2382`,
+  :pull:`3498`) by `Julien Seguinot <https://github.com/juseg>_`.
 - :py:meth:`Dataset.swap_dims` and :py:meth:`DataArray.swap_dims`
   now allow swapping to dimension names that don't exist yet. (:pull:`3636`)
   By `Justus Magin <https://github.com/keewis>`_.
diff --git a/xarray/backends/api.py b/xarray/backends/api.py
index 23d09ba5..eea1fb9d 100644
--- a/xarray/backends/api.py
+++ b/xarray/backends/api.py
@@ -718,6 +718,7 @@ def open_mfdataset(
     autoclose=None,
     parallel=False,
     join="outer",
+    attrs_file=None,
     **kwargs,
 ):
     """Open multiple files as a single dataset.
@@ -729,8 +730,8 @@ def open_mfdataset(
     ``combine_by_coords`` and ``combine_nested``. By default the old (now deprecated)
     ``auto_combine`` will be used, please specify either ``combine='by_coords'`` or
     ``combine='nested'`` in future. Requires dask to be installed. See documentation for
-    details on dask [1]_. Attributes from the first dataset file are used for the
-    combined dataset.
+    details on dask [1]_. Global attributes from the ``attrs_file`` are used
+    for the combined dataset.
 
     Parameters
     ----------
@@ -827,6 +828,10 @@ def open_mfdataset(
         - 'override': if indexes are of same size, rewrite indexes to be
           those of the first object with that dimension. Indexes for the same
           dimension must have the same size in all objects.
+    attrs_file : str or pathlib.Path, optional
+        Path of the file used to read global attributes from.
+        By default global attributes are read from the first file provided,
+        with wildcard matches sorted by filename.
     **kwargs : optional
         Additional arguments passed on to :py:func:`xarray.open_dataset`.
 
@@ -961,7 +966,15 @@ def open_mfdataset(
         raise
 
     combined._file_obj = _MultiFileCloser(file_objs)
-    combined.attrs = datasets[0].attrs
+
+    # read global attributes from the attrs_file or from the first dataset
+    if attrs_file is not None:
+        if isinstance(attrs_file, Path):
+            attrs_file = str(attrs_file)
+        combined.attrs = datasets[paths.index(attrs_file)].attrs
+    else:
+        combined.attrs = datasets[0].attrs
+
     return combined
 
 


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index a23527bd..4fccdf2d 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -2832,6 +2832,42 @@ class TestDask(DatasetIOBase):
                     with raises_regex(AttributeError, "no attribute"):
                         actual.test2
 
+    def test_open_mfdataset_attrs_file(self):
+        original = Dataset({"foo": ("x", np.random.randn(10))})
+        with create_tmp_files(2) as (tmp1, tmp2):
+            ds1 = original.isel(x=slice(5))
+            ds2 = original.isel(x=slice(5, 10))
+            ds1.attrs["test1"] = "foo"
+            ds2.attrs["test2"] = "bar"
+            ds1.to_netcdf(tmp1)
+            ds2.to_netcdf(tmp2)
+            with open_mfdataset(
+                [tmp1, tmp2], concat_dim="x", combine="nested", attrs_file=tmp2
+            ) as actual:
+                # attributes are inherited from the master file
+                assert actual.attrs["test2"] == ds2.attrs["test2"]
+                # attributes from ds1 are not retained, e.g.,
+                assert "test1" not in actual.attrs
+
+    def test_open_mfdataset_attrs_file_path(self):
+        original = Dataset({"foo": ("x", np.random.randn(10))})
+        with create_tmp_files(2) as (tmp1, tmp2):
+            tmp1 = Path(tmp1)
+            tmp2 = Path(tmp2)
+            ds1 = original.isel(x=slice(5))
+            ds2 = original.isel(x=slice(5, 10))
+            ds1.attrs["test1"] = "foo"
+            ds2.attrs["test2"] = "bar"
+            ds1.to_netcdf(tmp1)
+            ds2.to_netcdf(tmp2)
+            with open_mfdataset(
+                [tmp1, tmp2], concat_dim="x", combine="nested", attrs_file=tmp2
+            ) as actual:
+                # attributes are inherited from the master file
+                assert actual.attrs["test2"] == ds2.attrs["test2"]
+                # attributes from ds1 are not retained, e.g.,
+                assert "test1" not in actual.attrs
+
     def test_open_mfdataset_auto_combine(self):
         original = Dataset({"foo": ("x", np.random.randn(10)), "x": np.arange(10)})
         with create_tmp_file() as tmp1: