def test_merge_datasets(self):
        data = create_test_data()

        actual = xr.merge([data[['var1']], data[['var2']]])
        expected = data[['var1', 'var2']]
        assert actual.identical(expected)

        actual = xr.merge([data, data])
        assert actual.identical(data)

    def test_merge_no_conflicts_broadcast(self):
        datasets = [xr.Dataset({'x': ('y', [0])}), xr.Dataset({'x': np.nan})]
        actual = xr.merge(datasets)
        expected = xr.Dataset({'x': ('y', [0])})
        assert expected.identical(actual)

        datasets = [xr.Dataset({'x': ('y', [np.nan])}), xr.Dataset({'x': 0})]
        actual = xr.merge(datasets)
        assert expected.identical(actual)

    def test_merge_no_conflicts_single_var(self):
        ds1 = xr.Dataset({'a': ('x', [1, 2]), 'x': [0, 1]})
        ds2 = xr.Dataset({'a': ('x', [2, 3]), 'x': [1, 2]})
        expected = xr.Dataset({'a': ('x', [1, 2, 3]), 'x': [0, 1, 2]})
        assert expected.identical(xr.merge([ds1, ds2],
                                           compat='no_conflicts'))
        assert expected.identical(xr.merge([ds2, ds1],
                                           compat='no_conflicts'))
        assert ds1.identical(xr.merge([ds1, ds2],
                                      compat='no_conflicts',
                                      join='left'))
        assert ds2.identical(xr.merge([ds1, ds2],
                                      compat='no_conflicts',
                                      join='right'))
        expected = xr.Dataset({'a': ('x', [2]), 'x': [1]})
        assert expected.identical(xr.merge([ds1, ds2],
                                           compat='no_conflicts',
                                           join='inner'))

        with pytest.raises(xr.MergeError):
            ds3 = xr.Dataset({'a': ('x', [99, 3]), 'x': [1, 2]})
            xr.merge([ds1, ds3], compat='no_conflicts')

        with pytest.raises(xr.MergeError):
            ds3 = xr.Dataset({'a': ('y', [2, 3]), 'y': [1, 2]})
            xr.merge([ds1, ds3], compat='no_conflicts')

def rinexobs2(fn: Path,
              use: Sequence[str] = None,
              tlim: Tuple[datetime, datetime] = None,
              useindicators: bool = False,
              meas: Sequence[str] = None,
              verbose: bool = False,
              *,
              fast: bool = True,
              interval: Union[float, int, timedelta] = None) -> xarray.Dataset:

    if isinstance(use, str):
        use = [use]

    if use is None or not use[0].strip():
        use = ('C', 'E', 'G', 'J', 'R', 'S')

    obs = xarray.Dataset({}, coords={'time': [], 'sv': []})
    attrs: Dict[str, Any] = {}
    for u in use:
        o = rinexsystem2(fn, system=u, tlim=tlim,
                         useindicators=useindicators, meas=meas,
                         verbose=verbose,
                         fast=fast, interval=interval)
        if len(o) > 0:
            attrs = o.attrs
            obs = xarray.merge((obs, o))

    obs.attrs = attrs

    return obs

    def __init__(self, **kwds):
        node_at_cell = kwds.pop("node_at_cell", None)
        nodes_at_face = kwds.pop("nodes_at_face", None)

        update_node_at_cell(self.ds, node_at_cell)
        update_nodes_at_face(self.ds, nodes_at_face)

        rename = {
            "mesh": "dual",
            "node": "corner",
            "link": "face",
            "patch": "cell",
            "x_of_node": "x_of_corner",
            "y_of_node": "y_of_corner",
            "nodes_at_link": "corners_at_face",
            "links_at_patch": "faces_at_cell",
            "max_patch_links": "max_cell_faces",
        }
        self._ds = xr.merge([self._ds, self._dual.ds.rename(rename)])

        self._origin = (0.0, 0.0)

        self._frozen = False
        self.freeze()

        if kwds.get("sort", True):
            self.sort()

    def __init__(self, **kwds):
        node_at_cell = kwds.pop('node_at_cell', None)
        nodes_at_face = kwds.pop('nodes_at_face', None)

        update_node_at_cell(self.ds, node_at_cell)
        update_nodes_at_face(self.ds, nodes_at_face)

        rename = {
            'mesh': 'dual',
            'node': 'corner',
            'link': 'face',
            'patch': 'cell',
            'x_of_node': 'x_of_corner',
            'y_of_node': 'y_of_corner',
            'nodes_at_link': 'corners_at_face',
            'links_at_patch': 'faces_at_cell',
            'max_patch_links': 'max_cell_faces',
        }
        self._ds = xr.merge([self._ds, self._dual.ds.rename(rename)])

        self._origin = (0., 0.)

        self._frozen = False
        self.freeze()

        if kwds.get('sort', True):
            self.sort()

 def test_merge_returns_merged_array(self):
     test_array = self.array.copy()
     test_array.name = 'test'
     test_dataset = xr.merge([test_array, self.array])
     merged_array = test_dataset.to_array(name='merged_array')
     self.array.pp.grid = self.grid
     returned_array = self.array.pp.merge(test_array)
     np.testing.assert_equal(merged_array.values, returned_array.values)

    def test_merge_no_conflicts_multi_var(self):
        data = create_test_data()
        data1 = data.copy(deep=True)
        data2 = data.copy(deep=True)

        expected = data[['var1', 'var2']]
        actual = xr.merge([data1.var1, data2.var2], compat='no_conflicts')
        assert expected.identical(actual)

        data1['var1'][:, :5] = np.nan
        data2['var1'][:, 5:] = np.nan
        data1['var2'][:4, :] = np.nan
        data2['var2'][4:, :] = np.nan
        del data2['var3']

        actual = xr.merge([data1, data2], compat='no_conflicts')
        assert data.equals(actual)

    def test_merge_attr_retention(self):
        da1 = create_test_dataarray_attrs(var='var1')
        da2 = create_test_dataarray_attrs(var='var2')
        da2.attrs = {'wrong': 'attributes'}
        original_attrs = da1.attrs

        # merge currently discards attrs, and the global keep_attrs
        # option doesn't affect this
        result = merge([da1, da2])
        assert result.attrs == original_attrs

    def convert_species_bit_flags(data: Dict) -> xr.Dataset:
        """
        Convert the int32 species flags to a dataset of distinct flags

        Parameters
        ----------
        data
            Dictionary of input data as returned by `load_data`

        Returns
        -------
            Dataset of the index bit flags
        """
        flags = dict()
        flags['separation_method'] = [0, 1, 2]
        flags['one_chan_aerosol_corr'] = 3
        flags['no_935_aerosol_corr'] = 4
        flags['Large_1020_OD'] = 5
        flags['NO2_Extrap'] = 6
        flags['Water_vapor_ratio'] = [7, 8, 9, 10]
        flags['Cloud_Bit_1'] = 11
        flags['Cloud_Bit_2'] = 12
        flags['No_H2O_Corr'] = 13
        flags['In_Troposphere'] = 14

        separation_method = dict()
        separation_method['no_aerosol_method'] = 0
        separation_method['trans_no_aero_to_five_chan'] = 1
        separation_method['standard_method'] = 2
        separation_method['trans_five_chan_to_low'] = 3
        separation_method['four_chan_method'] = 4
        separation_method['trans_four_chan_to_three_chan'] = 5
        separation_method['three_chan_method'] = 6
        separation_method['extension_method'] = 7

        f = dict()
        for key in flags.keys():
            if hasattr(flags[key], '__len__'):
                if key == 'separation_method':
                    for k in separation_method.keys():
                        temp = data['ProfileInfVec'] & np.sum([2 ** k for k in flags[key]])
                        f[k] = temp == separation_method[k]
                else:
                    temp = data['ProfileInfVec'] & np.sum([2 ** k for k in flags[key]])
                    f[key] = temp >> flags[key][0]  # shift flag to save only significant bits
            else:
                f[key] = (data['ProfileInfVec'] & 2 ** flags[key]) > 0

        xr_data = []
        time = pd.to_timedelta(data['mjd'], 'D') + pd.Timestamp('1858-11-17')
        for key in f.keys():
            xr_data.append(xr.DataArray(f[key], coords=[time, data['Alt_Grid'][0:140]], dims=['time', 'Alt_Grid'],
                                        name=key))

        return xr.merge(xr_data)

def nan_ds():
    def make_data(a, b):
        return np.cos(a) + b, a + np.sin(b)
    r = Runner(make_data, ('a10sum', 'b10sum'))
    ds1 = r.run_combos((('a', np.linspace(1, 3, 10)),
                        ('b', np.linspace(1, 3, 10))))
    ds2 = r.run_combos((('a', np.linspace(1.5, 3.5, 10)),
                        ('b', np.linspace(1, 3, 10))))
    ds3 = r.run_combos((('a', np.linspace(4, 6, 10)),
                        ('b', np.linspace(4, 6, 10))))
    return xr.merge([ds1, ds2, ds3])

def merge(ds_1: DatasetLike.TYPE,
          ds_2: DatasetLike.TYPE,
          ds_3: DatasetLike.TYPE = None,
          ds_4: DatasetLike.TYPE = None,
          join: str = 'outer',
          compat: str = 'no_conflicts') -> xr.Dataset:
    """
    Merge up to four datasets to produce a new dataset with combined variables from each input dataset.

    This is a wrapper for the ``xarray.merge()`` function.

    For documentation refer to xarray documentation at
    http://xarray.pydata.org/en/stable/generated/xarray.Dataset.merge.html#xarray.Dataset.merge

    The *compat* argument indicates how to compare variables of the same name for potential conflicts:

    * "broadcast_equals": all values must be equal when variables are broadcast
      against each other to ensure common dimensions.
    * "equals": all values and dimensions must be the same.
    * "identical": all values, dimensions and attributes must be the same.
    * "no_conflicts": only values which are not null in both datasets must be equal.
      The returned dataset then contains the combination of all non-null values.

    :param ds_1: The first input dataset.
    :param ds_2: The second input dataset.
    :param ds_3: An optional 3rd input dataset.
    :param ds_4: An optional 4th input dataset.
    :param join: How to combine objects with different indexes.
    :param compat: How to compare variables of the same name for potential conflicts.
    :return: A new dataset with combined variables from each input dataset.
    """

    ds_1 = DatasetLike.convert(ds_1)
    ds_2 = DatasetLike.convert(ds_2)
    ds_3 = DatasetLike.convert(ds_3)
    ds_4 = DatasetLike.convert(ds_4)

    datasets = []
    for ds in (ds_1, ds_2, ds_3, ds_4):
        if ds is not None:
            included = False
            for ds2 in datasets:
                if ds is ds2:
                    included = True
            if not included:
                datasets.append(ds)

    if len(datasets) == 0:
        raise ValidationError('At least two different datasets must be given')
    elif len(datasets) == 1:
        return datasets[0]
    else:
        return xr.merge(datasets, compat=compat, join=join)

 def test_merge_fill_value(self, fill_value):
     ds1 = xr.Dataset({'a': ('x', [1, 2]), 'x': [0, 1]})
     ds2 = xr.Dataset({'b': ('x', [3, 4]), 'x': [1, 2]})
     if fill_value == dtypes.NA:
         # if we supply the default, we expect the missing value for a
         # float array
         fill_value = np.nan
     expected = xr.Dataset({'a': ('x', [1, 2, fill_value]),
                            'b': ('x', [fill_value, 3, 4])},
                           {'x': [0, 1, 2]})
     assert expected.identical(ds1.merge(ds2, fill_value=fill_value))
     assert expected.identical(ds2.merge(ds1, fill_value=fill_value))
     assert expected.identical(xr.merge([ds1, ds2], fill_value=fill_value))

def __get_maximum_storage_and_corresponding_dates(start_year:int, end_year:int, data_manager:DataManager, storage_varname=""):
    cache_file_current = "cache_{}-{}_calculate_flood_storage_{}.nc".format(start_year, end_year, storage_varname)
    cache_file_current = Path(cache_file_current)

    # if the variables were calculated already
    if cache_file_current.exists():
        ds = xarray.open_dataset(str(cache_file_current))
    else:
        data_current = data_manager.get_min_max_avg_for_period(
            start_year=start_year, end_year=end_year, varname_internal=storage_varname
        )

        ds = xarray.merge([da for da in data_current.values()])
        ds.to_netcdf(str(cache_file_current))

    return ds

def merge_sync_conflict_datasets(base_name, engine='h5netcdf',
                                 combine_first=False):
    """Glob files based on `base_name`, merge them, save this new dataset if
    it contains new info, then clean up the conflicts.

    Parameters
    ----------
        base_name : str
            Base file name to glob on - should include '*'.
        engine : str , optional
            Load and save engine used by xarray.
    """
    fnames = glob(base_name)
    if len(fnames) < 2:
        print('Nothing to do - need multiple files to merge.')
        return

    # make sure first filename is the shortest -> assumed original
    fnames.sort(key=len)

    print("Merging:\n{}\ninto ->\n{}\n".format(fnames, fnames[0]))

    def load_dataset(fname):
        return load_ds(fname, engine=engine)

    datasets = list(map(load_dataset, fnames))

    # combine all the conflicts
    if combine_first:
        full_dataset = datasets[0]
        for ds in datasets[1:]:
            full_dataset = full_dataset.combine_first(ds)
    else:
        full_dataset = xr.merge(datasets)

    # save new dataset?
    if full_dataset.identical(datasets[0]):
        # nothing to do
        pass
    else:
        save_ds(full_dataset, fnames[0], engine=engine)

    # clean up conflicts
    for fname in fnames[1:]:
        os.remove(fname)

def _epoch(data: xarray.Dataset, raw: str,
           hdr: Dict[str, Any],
           time: datetime,
           sv: List[str],
           useindicators: bool,
           verbose: bool) -> xarray.Dataset:
    """
    block processing of each epoch (time step)
    """
    darr = np.atleast_2d(np.genfromtxt(io.BytesIO(raw.encode('ascii')),
                                       delimiter=(14, 1, 1) * hdr['Fmax']))
# %% assign data for each time step
    for sk in hdr['fields']:  # for each satellite system type (G,R,S, etc.)
        # satellite indices "si" to extract from this time's measurements
        si = [i for i, s in enumerate(sv) if s[0] in sk]
        if len(si) == 0:  # no SV of this system "sk" at this time
            continue

        # measurement indices "di" to extract at this time step
        di = hdr['fields_ind'][sk]
        garr = darr[si, :]
        garr = garr[:, di]

        gsv = np.array(sv)[si]

        dsf: Dict[str, tuple] = {}
        for i, k in enumerate(hdr['fields'][sk]):
            dsf[k] = (('time', 'sv'), np.atleast_2d(garr[:, i*3]))

            if useindicators:
                dsf = _indicators(dsf, k, garr[:, i*3+1:i*3+3])

        if verbose:
            print(time, '\r', end='')

        epoch_data = xarray.Dataset(dsf, coords={'time': [time], 'sv': gsv})
        if len(data) == 0:
            data = epoch_data
        elif len(hdr['fields']) == 1:  # one satellite system selected, faster to process
            data = xarray.concat((data, epoch_data), dim='time')
        else:  # general case, slower for different satellite systems all together
            data = xarray.merge((data, epoch_data))

    return data

    def convert_index_bit_flags(data: Dict) -> xr.Dataset:
        """
        Convert the int32 index flags to a dataset of distinct flags

        Parameters
        ----------
        data
            Dictionary of input data as returned by ``load_data``

        Returns
        -------
            Dataset of the index bit flags
        """
        flags = dict()
        flags['pmc_present'] = 0
        flags['h2o_zero_found'] = 1
        flags['h2o_slow_convergence'] = 2
        flags['h2o_ega_failure'] = 3
        flags['default_nmc_temp_errors'] = 4
        flags['ch2_aero_model_A'] = 5
        flags['ch2_aero_model_B'] = 6
        flags['ch2_new_wavelength'] = 7
        flags['incomplete_nmc_data'] = 8
        flags['mirror_model'] = 15
        flags['twomey_non_conv_rayleigh'] = 19
        flags['twomey_non_conv_386_Aero'] = 20
        flags['twomey_non_conv_452_Aero'] = 21
        flags['twomey_non_conv_525_Aero'] = 22
        flags['twomey_non_conv_1020_Aero'] = 23
        flags['twomey_non_conv_NO2'] = 24
        flags['twomey_non_conv_ozone'] = 25
        flags['no_shock_correction'] = 30

        f = dict()
        for key in flags.keys():
            f[key] = (data['InfVec'] & 2 ** flags[key]) > 0

        xr_data = []
        time = pd.to_timedelta(data['mjd'], 'D') + pd.Timestamp('1858-11-17')
        for key in f.keys():
            xr_data.append(xr.DataArray(f[key], coords=[time], dims=['time'], name=key))

        return xr.merge(xr_data)

def merge_data(data, merge_attrs=True):
    """ Combine multiple Datasets or DataArrays into one Dataset

    Args:
        data (dict[name, xr.DataArray or xr.Dataset]): xr.DataArray or
            xr.Dataset objects to merge
        merge_attrs (bool): Attempt to merge DataArray attributes. In order for
            these attributes to be able to merge, they must be pd.Series
            and have compatible indexes.

    Returns:
        xr.Dataset: Merged xr.DataArray objects in one xr.Dataset
    """
    datasets = [dat.to_dataset(dim='band') if isinstance(dat, xr.DataArray)
                else dat for dat in data.values()]

    ds = xr.merge(datasets, compat='minimal')

    # Overlapping but not conflicting variables can't be merged for now
    # https://github.com/pydata/xarray/issues/835
    # In meantime, check for dropped variables
    dropped = set()
    for _ds in datasets:
        dropped.update(set(_ds.data_vars).difference(set(ds.data_vars)))

    # TODO: refactor this and repeat for coords and vars
    if dropped:
        # dropped_vars = {}
        for var in dropped:
            dims = [_ds[var].dims for _ds in datasets]
            if all([dims[0] == d for d in dims[1:]]):
                dfs = [_ds.reset_coords()[var].to_series()
                       for _ds in datasets]
                # dropped_vars[var] = (dims[0], pd.concat(dfs).sort_index())
                da = xr.DataArray(pd.concat(dfs).sort_index())
                ds[var] = da
            else:
                logger.debug("Cannot restore dropped coord {} because "
                             "dimensions are inconsistent across datasets")
        # ds = ds.assign_coords(**ds.coords.merge(dropped_vars))

    return ds

def xr2mat(fields, sample_dims, feature_dims,
           scale=True):
    """Prepare list of data arrays for input to Machine Learning

    Parameters
    ----------
    fields: Dataset object
        input dataset
    sample_dims: tuple of string
        Dimensions which will be considered samples
    feature_dims: tuple of strings
        dimensions which be considered features
    scale: Bool
        center and scale the output. [default: True]

    Returns
    -------
    data: DataArray
    scaling: DataArray or None

    See Also
    --------
    gnl.data_matrix.DataMatrix : better version of this function
    """
    raise DeprecationWarning("Please switch to gnl.data_matrix.DataMatrix")
    normalize_dim = 'z'   # this could be an argument

    if not isinstance(fields, xr.Dataset):
        fields = xr.merge(fields)

    dat = fields.to_array()

    if scale:
        mu = dat.mean(sample_dims)
        V = np.sqrt(((dat-mu)**2).sum(normalize_dim)).mean(sample_dims)
        dat = (dat-mu)/V
    else:
        V = None

    stacked = dat.stack(features=('variable',)+tuple(feature_dims),
                        samples=sample_dims)
    return stacked.transpose('samples', 'features'), V

 def test_merge_alignment_error(self):
     ds = xr.Dataset(coords={'x': [1, 2]})
     other = xr.Dataset(coords={'x': [2, 3]})
     with raises_regex(ValueError, 'indexes .* not equal'):
         xr.merge([ds, other], join='exact')