The xr.Dataset.expand_dims() method can be used to add new dimensions to a dataset. The axis parameter is used to specify where to insert the new dimension in the dataset. However, it's worth noting that the axis parameter only works when expanding along a 1D coordinate, not when expanding along a multi-dimensional array.

Here's an example to illustrate how to use the axis parameter to expand a dataset along a 1D coordinate:

import xarray as xr

create a sample dataset

data = xr.DataArray([[1, 2], [3, 4]], dims=('x', 'y')) ds = xr.Dataset({'foo': data})

add a new dimension along the 'x' coordinate using the 'axis' parameter

ds_expanded = ds.expand_dims({'z': [1]}, axis='x')

In this example, we create a 2D array with dimensions x and y, and then add a new dimension along the x coordinate using the axis='x' parameter.

However, if you try to use the axis parameter to expand a dataset along a multi-dimensional array, you may encounter an error. This is because expanding along a multi-dimensional array would result in a dataset with non-unique dimension names, which is not allowed in xarray.

Here's an example to illustrate this issue:

import xarray as xr

create a sample dataset with a 2D array

data = xr.DataArray([[1, 2], [3, 4]], dims=('x', 'y')) ds = xr.Dataset({'foo': data})

add a new dimension along the 'x' and 'y' coordinates using the 'axis' parameter

ds_expanded = ds.expand_dims({'z': [1]}, axis=('x', 'y'))

In this example, we try to use the axis=('x', 'y') parameter to add a new dimension along both the x and y coordinates. However, this results in a ValueError because the resulting dataset would have non-unique dimension names.

To add a new dimension along a multi-dimensional array, you can instead use the xr.concat() function to concatenate the dataset with a new data array along the desired dimension:

import xarray as xr

create a sample dataset with a 2D array

data = xr.DataArray([[1, 2], [3, 4]], dims=('x', 'y')) ds = xr.Dataset({'foo': data})

add a new dimension along the 'x' and 'y' coordinates using xr.concat

ds_expanded = xr.concat([ds, xr.DataArray([1], dims=('z'))], dim='z')

In this example, we use the xr.concat() function to concatenate the original dataset with a new data array that has a single value along the new dimension z. The dim='z' parameter is used to specify that the new dimension should be named z.

If you are encountering an error message that says "Plotting with time-zone-aware pd.Timestamp axis not possible", it means that you are trying to plot a Pandas DataFrame or Series that has a time-zone-aware pd.Timestamp axis using a plotting library that does not support time zones.

To fix this error, you can convert the time-zone-aware pd.Timestamp axis to a time-zone-naive datetime object. This can be done using the tz_localize() method to set the time zone, followed by the tz_convert() method to convert to a new time zone or remove the time zone information altogether.

Here is an example:

import pandas as pd import matplotlib.pyplot as plt

Create a time-series DataFrame with a time-zone-aware pd.Timestamp axis

data = pd.DataFrame({'value': [1, 2, 3, 4]}, index=pd.date_range('2022-03-01 00:00:00', periods=4, freq='H', tz='US/Eastern'))

Convert the time-zone-aware pd.Timestamp axis to a time-zone-naive datetime object

data.index = data.index.tz_localize(None)

Plot the DataFrame using Matplotlib

data.plot() plt.show()

In this example, we create a time-series DataFrame with a time-zone-aware pd.Timestamp axis using the pd.date_range() function with the tz parameter set to 'US/Eastern'. We then use the tz_localize() method to set the time zone to None to convert the axis to a time-zone-naive datetime object. Finally, we plot the DataFrame using Matplotlib and the plot() method.

Note that converting the time-zone-aware pd.Timestamp axis to a time-zone-naive datetime object means that the time zone information is lost, so make sure that this is acceptable for your use case before making this conversion.

When using NumPy arrays, the np.multiply() function and the * operator behave the same way and perform element-wise multiplication on the arrays. Similarly, the np.add() function and the + operator perform element-wise addition.

However, when using Dask arrays, there is a difference between using the * and + operators and using the dask.array.multiply() and dask.array.add() functions. This is because Dask arrays are lazy and do not compute the result of an operation until it is explicitly requested. When you use the * or + operators, Dask constructs a task graph that describes the computation, but does not actually execute it until you explicitly call a computation method like dask.compute() or dask.persist().

On the other hand, when you use the dask.array.multiply() or dask.array.add() functions, Dask immediately constructs a task graph and adds it to the computation graph, triggering the computation to begin.

Here's an example to illustrate the difference:

import dask.array as da

x = da.ones((1000, 1000), chunks=(100, 100)) y = da.ones((1000, 1000), chunks=(100, 100))

using the * operator

z = x * y

no computation is triggered yet

using dask.array.multiply()

z = da.multiply(x, y)

computation is immediately triggered

In this example, the * operator creates a task graph for the multiplication but does not execute it, whereas the dask.array.multiply() function immediately adds the task graph to the computation graph and triggers the computation to begin.

It's worth noting that using the * and + operators can be more convenient and can lead to cleaner code, especially for simple operations. However, if you need more control over when computations are executed or want to avoid unnecessary computations, you should use the dask.array.multiply() and dask.array.add() functions.

@jonas-constellr

It's possible that the failure you're experiencing is due to an issue with how the h5netcdf library is interacting with Dask.

One potential solution to this issue is to try using the netCDF4 library instead of h5netcdf. netCDF4 is another popular library for reading and writing netCDF files, and it has built-in support for parallel I/O through Dask.

To use netCDF4 with xarray, you can simply pass the 'netcdf4' engine to the xr.open_mfdataset function:

python import xarray as xr

Open multiple netCDF files with netCDF4 engine and parallel I/O

ds = xr.open_mfdataset('path/to/files/*.nc', engine='netcdf4', parallel=True) If you need to use h5netcdf for some reason, another potential solution is to use the dask.array.from_delayed function to manually create a Dask array from the h5netcdf data. This can be done by first reading in the data using h5netcdf, and then using dask.delayed to parallelize the data loading across multiple chunks. Here's an example:

python

import h5netcdf import dask.array as da from dask import delayed

Define function to read in a single chunk of data from the netCDF file

@delayed def read_chunk(filename, varname, start, count): with h5netcdf.File(filename, 'r') as f: var = f[varname][start[0]:start[0]+count[0], start[1]:start[1]+count[1]] return var

Define function to read in the entire dataset using dask.array.from_delayed

def read_data(files, varname): chunks = (1000, 1000) # Define chunk size data = [read_chunk(f, varname, start, chunks) for f in files] data = [da.from_delayed(d, shape=chunks, dtype='float64') for d in data] data = da.concatenate(data, axis=0) return data

Open multiple netCDF files with h5netcdf engine and parallel I/O

files = ['path/to/files/file1.nc', 'path/to/files/file2.nc', ...] varname = 'my_variable' data = read_data(files, varname) This code reads in the data from each file in chunks, and returns a Dask array that is a concatenation of all the chunks. The read_chunk function uses h5netcdf.File to read in a single chunk of data from a file, and returns a delayed object that represents the loading of that chunk. The read_data function uses dask.delayed to parallelize the loading of the chunks across all the files, and then uses dask.array.from_delayed to create a Dask array from the delayed objects. Finally, the function returns the concatenated Dask array.

What happened?

I was trying to read multiple byte netcdf (requires h5netcdf engine) file with xr.open_mfdataset with parallel=True to leverage dask.delayed capabilities (parallel=False works though) but it failed.

The netcdf files were noaa-goes16 satellite images, but I can't tell if it matters.

What did you expect to happen?

It should have loaded all the netcdf files into a xarray.DataSet object

Minimal Complete Verifiable Example

```

Python

import fsspec

import xarray as xr

paths = [

's3://noaa-goes16/ABI-L2-LSTC/2022/185/03/OR_ABI-L2-LSTC-M6_G16_s20221850301180_e20221850303553_c20221850305091.nc',

's3://noaa-goes16/ABI-L2-LSTC/2022/185/02/OR_ABI-L2-LSTC-M6_G16_s20221850201180_e20221850203553_c20221850205142.nc'

]

fs = fsspec.filesystem('s3')

xr.open_mfdataset(

[fs.open(path, mode="rb") for path in paths],

engine="h5netcdf",

combine="nested",

concat_dim="t",

parallel=True

)

```

MVCE confirmation

• [ ] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.

• [X] Complete example — the example is self-contained, including all data and the text of any traceback.

• [ ] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.

• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

```Python

---

KeyError Traceback (most recent call last)

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/xarray/backends/file_manager.py:210, in CachingFileManager._acquire_with_cache_info(self, needs_lock)

209 try:

--> 210 file = self._cache[self._key]

211 except KeyError:

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.getitem(self, key)

 55 with self._lock:

---> 56 value = self._cache[key]

 57     self._cache.move_to_end(key)

?[0;31mKeyError?[0m: [<class 'h5netcdf.core.File'>, ((b'\x89HDF\r\n', b'\x1a\n', b'\x02\x08\x08\x00\x00\x00 ... EXTREMELY STRING ... 00\x00\x00\x00\x00\x00\x0ef']

During handling of the above exception, another exception occurred:

TypeError Traceback (most recent call last)

Cell In[9], line 11

  4 paths = [

  5     's3://noaa-goes16/ABI-L2-LSTC/2022/185/03/OR_ABI-L2-LSTC-M6_G16_s20221850301180_e20221850303553_c20221850305091.nc',

  6     's3://noaa-goes16/ABI-L2-LSTC/2022/185/02/OR_ABI-L2-LSTC-M6_G16_s20221850201180_e20221850203553_c20221850205142.nc'

  7     ]

  9 fs = fsspec.filesystem('s3')

---> 11 xr.open_mfdataset(

 12     [fs.open(path, mode="rb") for path in paths],

 13     engine="h5netcdf",

 14     combine="nested",

 15     concat_dim="t",

 16     parallel=True

 17 ).LST

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/xarray/backends/api.py:991, in open_mfdataset(paths, chunks, concat_dim, compat, preprocess, engine, data_vars, coords, combine, parallel, join, attrs_file, combine_attrs, **kwargs)

986     datasets = [preprocess(ds) for ds in datasets]

988 if parallel:

989     # calling compute here will return the datasets/file_objs lists,

990     # the underlying datasets will still be stored as dask arrays

--> 991 datasets, closers = dask.compute(datasets, closers)

993 # Combine all datasets, closing them in case of a ValueError

994 try:

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/base.py:599, in compute(traverse, optimize_graph, scheduler, get, args, *kwargs)

596     keys.append(x.__dask_keys__())

597     postcomputes.append(x.__dask_postcompute__())

--> 599 results = schedule(dsk, keys, **kwargs)

600 return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/threaded.py:89, in get(dsk, keys, cache, num_workers, pool, **kwargs)

 86     elif isinstance(pool, multiprocessing.pool.Pool):

 87         pool = MultiprocessingPoolExecutor(pool)

---> 89 results = get_async(

 90     pool.submit,

 91     pool._max_workers,

 92     dsk,

 93     keys,

 94     cache=cache,

 95     get_id=_thread_get_id,

 96     pack_exception=pack_exception,

 97     **kwargs,

 98 )

100 # Cleanup pools associated to dead threads

101 with pools_lock:

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/local.py:511, in get_async(submit, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, chunksize, **kwargs)

509         _execute_task(task, data)  # Re-execute locally

510     else:

--> 511 raise_exception(exc, tb)

512 res, worker_id = loads(res_info)

513 state["cache"][key] = res

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/local.py:319, in reraise(exc, tb)

317 if exc.__traceback__ is not tb:

318     raise exc.with_traceback(tb)

--> 319 raise exc

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/local.py:224, in execute_task(key, task_info, dumps, loads, get_id, pack_exception)

222 try:

223     task, data = loads(task_info)

--> 224 result = _execute_task(task, data)

225     id = get_id()

226     result = dumps((result, id))

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/core.py:119, in _execute_task(arg, cache, dsk)

115     func, args = arg[0], arg[1:]

116     # Note: Don't assign the subtask results to a variable. numpy detects

117     # temporaries by their reference count and can execute certain

118     # operations in-place.

--> 119 return func(*(_execute_task(a, cache) for a in args))

120 elif not ishashable(arg):

121     return arg

File ~/miniconda3/envs/rxr/lib/python3.11/site-packages/dask/utils.py:73, in apply(func, args, kwargs)

 42 """Apply a function given its positional and keyword arguments.

 43

 44 Equivalent to ``func(*args, **kwargs)``

(...)

...

---> 19 filename = fspath(filename)

 20     if sys.platform == "win32":

 21         if isinstance(filename, str):

TypeError: expected str, bytes or os.PathLike object, not tuple

```

Anything else we need to know?

No response

Environment

INSTALLED VERSIONS

---

commit: None

python: 3.11.0 | packaged by conda-forge | (main, Jan 15 2023, 05:44:48) [Clang 14.0.6 ]

python-bits: 64

OS: Darwin

OS-release: 21.6.0

machine: x86_64

processor: i386

byteorder: little

LC_ALL: None

LANG: None

LOCALE: (None, 'UTF-8')

libhdf5: 1.12.2

libnetcdf: None

xarray: 2023.2.0

pandas: 1.5.3

numpy: 1.24.2

scipy: 1.10.1

netCDF4: None

pydap: None

h5netcdf: 1.1.0

h5py: 3.8.0

Nio: None

zarr: None

cftime: None

nc_time_axis: None

PseudoNetCDF: None

rasterio: 1.3.6

cfgrib: None

iris: None

bottleneck: None

dask: 2023.2.1

distributed: 2023.2.1

matplotlib: 3.7.0

cartopy: 0.21.1

seaborn: 0.12.2

numbagg: None

fsspec: 2023.1.0

cupy: None

pint: None

sparse: None

flox: None

numpy_groupies: None

setuptools: 67.4.0

pip: 23.0.1

conda: None

pytest: 7.2.1

mypy: None

IPython: 8.10.0

sphinx: None

/Users/jo/miniconda3/envs/rxr/lib/python3.11/site-packages/_distutils_hack/init.py:33: UserWarning: Setuptools is replacing distutils.

warnings.warn("Setuptools is replacing distutils.")

Time offset arithmetic involves adding or subtracting a duration from a specific time to obtain a new time. This is commonly used when dealing with time zones or calculating time differences between two events.

For example, suppose the current time is 3:30 PM in New York City, which is in the Eastern Time Zone (ET). We want to calculate what time it is in Los Angeles, which is in the Pacific Time Zone (PT), considering the 3-hour time difference between the two zones.

To do this, we can use time offset arithmetic by subtracting 3 hours from the current time in ET:

3:30 PM ET - 3 hours = 12:30 PM PT

Therefore, the current time in Los Angeles is 12:30 PM.

Another example of time offset arithmetic is when calculating the duration between two events. Suppose an event starts at 9:00 AM and ends at 10:30 AM. We can calculate the duration of the event by subtracting the start time from the end time:

10:30 AM - 9:00 AM = 1 hour 30 minutes

Therefore, the event lasted for 1 hour and 30 minutes.

When using xr.merge with compat='minimal', the resulting Dataset may have unexpected behavior, including causing len to return wrong and possibly negative values.

This is because compat='minimal' mode allows for the merging of datasets with non-matching dimensions and coordinates. When this happens, the merged dataset may contain "empty" dimensions or coordinates that have lost their original values and attributes.

To avoid this issue, you can either use compat='override' mode, which will overwrite conflicting dimensions and coordinates, or manually align and concatenate the datasets before merging them.

Here's an example of how to manually align and concatenate two datasets before merging them:

import xarray as xr

create two sample datasets with different dimensions

ds1 = xr.Dataset({'foo': (['x', 'y'], [[1, 2], [3, 4]])}, coords={'x': [0, 1], 'y': [0, 1]}) ds2 = xr.Dataset({'bar': (['x', 'z'], [[5, 6], [7, 8]])}, coords={'x': [0, 1], 'z': [0, 1]})

align and concatenate the datasets along the 'x' dimension

ds1_aligned = ds1.sel(x=slice(None, 1)) ds2_aligned = ds2.sel(x=slice(1, None)) merged = xr.merge([ds1_aligned, ds2_aligned])

In this example, we first slice each dataset to remove any non-matching dimensions and coordinates, and then concatenate them along the 'x' dimension. Finally, we merge the aligned datasets to create a single merged dataset without any corrupted dimensions or coordinates.

The interpolate_na function is typically used to fill missing values (NAs) in a data frame or array by interpolating between existing values. It has an optional argument called max_gap which specifies the maximum number of consecutive NAs that can be filled in a single interpolation step.

However, the max_gap argument may not work as expected at the boundaries of an array, as there may not be enough data points available to fill the maximum gap. For example, if the max_gap is set to 3 and there are only two consecutive NAs at the boundary of an array, the function will not be able to fill those NAs.

One way to handle this issue is to reduce the max_gap value near the boundaries of the array. For example, you could set the max_gap to 1 for the first and last few rows or columns of the array, depending on the structure of your data. Alternatively, you could use a different interpolation method (e.g., linear interpolation) that does not require a fixed max_gap value.

It's also worth noting that the interpolate_na function may not always be the best approach for filling missing values, as it assumes that the data has a smooth, continuous structure. If your data has a more complex structure (e.g., sharp discontinuities), other methods such as regression or machine learning models may be more appropriate.

The interpolate_na method in xarray can be used to interpolate missing values in a dataset or data array. The max_gap argument is used to specify the maximum number of consecutive NaN values that can be interpolated. The max_map argument is used to specify the maximum number of interpolated values that can be used for each NaN value.

It's worth noting that the max_map argument only limits the number of interpolated values that can be used for each NaN value, but it does not limit the total number of interpolated values that can be used in the dataset. This means that if there are multiple consecutive NaN values, the max_map argument may not work as expected at the boundaries of the array.

Here's an example to illustrate this issue:

import xarray as xr import numpy as np

create a sample data array with a missing value at the beginning and end

data = np.array([np.nan, 1, 2, 3, 4, np.nan])

create a dataset with the sample data array

ds = xr.Dataset({'foo': (['x'], data)}, coords={'x': np.arange(6)})

interpolate missing values with a max_map of 2

ds_interp = ds.interpolate_na(max_gap=1, max_map=2)

In this example, we have a data array with missing values at the beginning and end, and we interpolate the missing values using a max_map of 2. However, the resulting dataset still has 4 interpolated values, which is more than the max_map of 2. This is because the max_map argument is only limiting the number of interpolated values that can be used for each NaN value, but it is not limiting the total number of interpolated values that can be used in the dataset.

To limit the total number of interpolated values in the dataset, you can use the limit argument, which specifies the maximum number of interpolated values that can be used in the entire dataset. Here's an example:

interpolate missing values with a max_map of 2 and a limit of 2

ds_interp = ds.interpolate_na(max_gap=1, max_map=2, limit=2)

In this example, we add a limit argument of 2, which limits the total number of interpolated values in the dataset to 2. This results in only 2 interpolated values in the resulting dataset, which is consistent with the limit argument.

@Ockenfuss i said you should try this three point I listed below and see if that could resolve the problem you raised.

1. Try adjusting the max_gap argument to a smaller value to see if that resolves the issue. For example, if max_gap is currently set to 10, try reducing it to 5 or even 1.

2. Consider using a different interpolation method that is better suited for the specific dataset and boundaries. For example, if linear interpolation is not working well at the array boundaries, try a cubic or spline interpolation method.

3. Check the data at the array boundaries to ensure that it is valid and not causing issues with the interpolation. For example, if there are NaN values or outliers at the boundaries, this could be affecting the interpolation.

https://github.com/pydata/xarray/issues/7593#issuecomment-1460894580