Scatter vs. Lines:

Just wanted to link to a somewhat related discussion happening in brian-rose/climlab#50.

Hi all, I just created a simple workaround, which might be useful for others:
https://gist.github.com/maawoo/0b34d371c3cc1960a1589ccaded868c2

It uses the _nan_quantile method of xclim and works fine for my applications. Here is a quick comparison using the same example data as in my initial post:

EDIT: I've updated the code to use numbagg instead of xclim.

Great observation @mathause . I think there are two parts of this: - Do we want other libraries which do da.longitude to raise a mypy error? That may be a tradeoff with raising the true error on da.isel - How do we design the type hierarchy? We could add methods to DataWithCoords or add some Dataset_Or_DataArray-like type

Having methods like isel typed would be a win I think

Good point re accessors, I hadn't considered those. So sounds like raising an error on da.isel isn't possible regardless...

Following thoughts and discussions in various issues (e.g., #6836), I'd like to suggest another section to the ones in the top comment:

Deprecate pandas.MultiIndex special cases in Xarray

• remove the multi-index “dimension” coordinate (tuple elements)
• do not automatically promote pandas.MultiIndex objects as dimension + level coordinates, e.g., like in xr.Dataset(coords={“x”: pd_midx}) but instead treat it as a single duck-array.
• do not accept pandas.MultiIndex as dim argument in xarray.concat() (#7148)
• remove obj.to_index() for all xarray objects?
• (EDIT) remove Dataset.reset_index() and DataArray.reset_index()

They are source of many problems and complexities in Xarray internals (many regressions reported since the index refactor were related to those special cases) and I'm not sure that the value they add is really worth the trouble. Also, in the long term the special treatment of PandasMultiIndex vs. other Xarray multi-indexes may add some confusion.

Some of those features are widely used (e.g., the creation of Dataset / DataArray from pandas multi-indexes is used in many places in unit tests), so we would need convenient alternatives and a smooth transition.

Actually it looks like this example is relevant: https://xarray.pydata.org/en/stable/examples/apply_ufunc_vectorize_1d.html

Hi dcherian,

I had a look at the apply_ufunc() example you linked and have re-implemented my code. The example helped me understand apply_ufunc() usage better but is very different from my use case and I still am unable to parallelize using dask.

The key difference is apply_ufunc() as described in the docs and the example, applys a function to a vector of data of a single type (in the example case it is air temperature across the 3 dimensions lat,long,time).

Where as I need to apply an operation using heterogeneous data (depth_bins, lower_limit, upper_limit) over a single dimension (time) to produce a new array of depths over time (which is why I tried groupby/map initially).

Anyhow, I have an implementation using apply_ufunc() that works using xarray and numpy arrays with apply_ufunc(), but when I try to parallelize it using dask my ufunc is called with empty arrays by xarray and it fails.

I.e. You can see when running the code below it logs the following when entering the ufunc: args: (array([], shape=(0, 0), dtype=int32), array([], dtype=int32), array([], dtype=int32), array([], dtype=int32)), kwargs: {}

I was expecting this to be called once for each chunk with 1000 items for each array.

Have I done something wrong in this work-around for the groupby/map code?

Thanks, Brendon

```python import sys import math import logging import dask import xarray import numpy

logger = logging.getLogger('main')

if name == 'main': logging.basicConfig( stream=sys.stdout, format='%(asctime)s %(levelname)-8s %(message)s', level=logging.INFO, datefmt='%Y-%m-%d %H:%M:%S')

logger.info('Starting dask client')
client = dask.distributed.Client()

SIZE = 3000
SONAR_BINS = 2000
upper_limit = numpy.random.randint(0, 10, (SIZE))
lower_limit = numpy.random.randint(20, 30, (SIZE))
sonar_data = numpy.random.randint(0, 255, (SIZE, SONAR_BINS))

time = range(0, SIZE)

channel = xarray.Dataset({
        'upper_limit': (['time'], upper_limit, {'units': 'depth meters'}),
        'lower_limit': (['time'],  lower_limit, {'units': 'depth meters'}),
        'data': (['time', 'depth_bin'], sonar_data, {'units': 'amplitude'}),
    },
    coords={
        'time': (['time'], time),
        'depth_bin': (['depth_bin'], range(0,SONAR_BINS)),
    })

logger.info('get overall min/max radar range we want to normalize to called the adjusted range')
adjusted_min, adjusted_max = channel.upper_limit.min().values.item(), channel.lower_limit.max().values.item()
adjusted_min = math.floor(adjusted_min)
adjusted_max = math.ceil(adjusted_max)
logger.info('adjusted_min: %s, adjusted_max: %s', adjusted_min, adjusted_max)

bin_count = len(channel.depth_bin)
logger.info('bin_count: %s', bin_count)

adjusted_depth_per_bin = (adjusted_max - adjusted_min) / bin_count
logger.info('adjusted_depth_per_bin: %s', adjusted_depth_per_bin)

adjusted_bin_depths = [adjusted_min + (j * adjusted_depth_per_bin) for j in range(0, bin_count)]
logger.info('adjusted_bin_depths[0]: %s ... [-1]: %s', adjusted_bin_depths[0], adjusted_bin_depths[-1])

def InterpSingle(unadjusted_depth_amplitudes, unadjusted_min, unadjusted_max, time):

    if (time % 1000) == 0:
        total = len(channel.time)
        perc = 100.0 * time / total
        logger.info('%s : %s of %s', perc, time, total)

    unadjusted_depth_per_bin = (unadjusted_max - unadjusted_min) / bin_count

    min_index = (adjusted_min - unadjusted_min) / unadjusted_depth_per_bin
    max_index = ((adjusted_min + ((bin_count - 1) * adjusted_depth_per_bin)) - unadjusted_min) / unadjusted_depth_per_bin
    index_mapping = numpy.linspace(min_index, max_index, bin_count)
    adjusted_depth_amplitudes = numpy.interp(index_mapping, range(0, len(unadjusted_depth_amplitudes)), unadjusted_depth_amplitudes, left=0, right=0)
    return adjusted_depth_amplitudes

def Interp(*args, **kwargs):
    logger.info('args: %s, kwargs: %s', args, kwargs)

    data = args[0]
    upper_limit = args[1]
    lower_limit = args[2]
    time = args[3]
    #logger.info('data: %s len(data[0]): %s data.shape: %s', data, len(data[0]), data.shape)

    adjusted = []
    for i in range(0, len(upper_limit)):
        d = data[i]
        u = upper_limit[i]
        l = lower_limit[i]
        t = time[i]
        result = InterpSingle(d, u, l, t)
        adjusted.append(result)

    #logger.info('adjusted: %s', adjusted)
    return adjusted



channel = channel.chunk({'time':1000}) # Comment this line out to disable dask
#logger.info('Channel: %s', channel)
#logger.info('shape of data: %s len(data[0]): %s', channel.data.shape, len(channel.data[0]))
m2_lazy = xarray.apply_ufunc(
    Interp, 
    channel.data, 
    channel.upper_limit, 
    channel.lower_limit, 
    channel.time, 
    input_core_dims=[['depth_bin'], [], [], []], 
    output_core_dims=[['depth']],
    dask='parallelized',   # Comment this line out to disable dask
    output_dtypes=[numpy.dtype(numpy.int32)],   # Comment this line out to disable dask
    output_sizes={'depth':len(adjusted_bin_depths)},   # Comment this line out to disable dask
    )
m2 = m2_lazy.compute() # Comment this line out to disable dask
m2 = m2.assign_coords({'depth':adjusted_bin_depths})

logger.info('Closing dask client')
client.close()

logger.info('Exit process')

```

https://github.com/pydata/xarray/issues/2233#issuecomment-397602084 Would the new xarray index/coordinate internal refactoring now allow us to address this issue?

cc @kthyng

If you want you can leave a comment But the mypy CI should fail on unused ignores, so we will notice :)

Might be worth an issue over at numpy with the example from the test.

numpy/numpy#23886

The issue is already resolved over at numpy which is really great! It was also marked as backport. @headtr1ck How are these issues resolved currently or how do we track removing the ignore?

Might be worth an issue over at numpy with the example from the test.

https://github.com/numpy/numpy/issues/23886

Should we delete the cfgrib example instead?

https://github.com/corteva/rioxarray/issues/676

```

input array

array = xr.DataArray([1,3,6,np.nan,19,20,13], dims=['time'], coords=[pd.date_range('2023-06-05 00:00','2023-06-05 06:00',freq='H')])

array xarray.DataArray(time: 7 array([ 1., 3., 6., nan, 19., 20., 13.]) Coordinates: time (time) datetime64[ns] 2023-06-05 ... 2023-06-05T06: Indexes: (1) Attributes: (0)

however the integrated value ends up as a NaN

array.integrate('time') xarray.DataArray array(nan) Coordinates: (0) Indexes: (0) Attributes: (0)

if one still wanted to know the integrated values for where there were values it would essentially by like integrating the separate chunks for where the valid values existed

first chunk

array.isel(time=slice(0,3)).integrate('time') xarray.DataArray array(2.34e+13) Coordinates: (0) Indexes: (0) Attributes: (0)

second chunk

array.isel(time=slice(4,7)).integrate('time') xarray.DataArray array(1.296e+14) Coordinates: (0) Indexes: (0) Attributes: (0)

and then the sum would be the fully integrated area

``` @dcherian I essentially was wondering whether it was possible for a skipna argument or some kind of NaN handling to be implemented that would allow users to avoid integrating in chunks due to the presence of NaNs. I do not work in dev so I would not know how to implement this, but I thought I'd see if others had thoughts.

but is it not possible for it to calculate the integrated values where there were regular values?

@chfite Can you provide an example of what you would want it to do please

You could use DataArray.round to round to significant decimals.

Better to use a tolerance in the assterion testing.

@vrishk Sorry for the delay here and thanks for bringing this to attention. We now have at least two requests which might move this forward (moving ensure_dtype_not_object into the backends). But this would need some discussion first, how to do this.

@mktippett Thanks for raising this. The issue should be cleared after #7888 is merged.

You could use DataArray.round to round to significant decimals.

Oh no, the doctest failure is because the test is flaky, this was introduced by me in #7821, see here: https://github.com/pydata/xarray/pull/7821#issuecomment-1537142237 and here: https://github.com/pydata/xarray/pull/7821/commits/a0e6659ca01188378f29a35b418d6f9e2b889d2e. I'll submit another patch to fix it soon, although I'm not sure how. If you have any tips to avoid this problem, let me know.

Thank you @keewis

Hi! I would also like to see this implemented, so I rebased this branch and added a test in #7891.

I think the only other one is dask, which should also work.

@dcherian : agreed! But I am afraid it might break other components. Although numpy seems to be able to handle both tuple and list in normalize_axis_tuple and I cannot see any other issues rising from this: https://github.com/numpy/numpy/blob/f67467c21a1797becde3097661996f60df4080ff/numpy/core/numeric.py#L1328

This seems like a real easy fix? axis = tuple(self.get_axis_num(d) for d in dim)

EDIT: the Array API seems to type axis as Optional[Union[int, Tuple[int, ...]]] pretty consistently, so it seems like we should always pass tuples down to the array library

Thanks for opening your first issue here at xarray! Be sure to follow the issue template! If you have an idea for a solution, we would really welcome a Pull Request with proposed changes. See the Contributing Guide for more. It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better. Thank you!

This seems to be a numpy issue, mypy thinks that you cannot call np.dtype like you do.

Might be worth an issue over at numpy with the example from the test.

For now we can simply ignore this error.

xarray/tests/test_coding_strings.py:36: error: No overload variant of "dtype" matches argument types "str", "Dict[str, Type[str]]" [call-overload]

cc @Illviljan @headtr1ck

@kmuehlbauer thanks for adding tests! I'm not sure what the mypy error is either, I'm afraid...

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!)

that appears to work on both my laptop and my local HPC, and is arguably a lot easier to implement / understand as we don't need to make sure all the globals we use are still available (which in this case would be acquire, OPTIONS, warnings, and RuntimeWarning).

Edit: let me change the PR to do that instead

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

I thought I did, but apparently something changed: right now it fails because OPTIONS is not available anymore, so we might have to add a reference to that, as well. Additionally, for the warning we use warnings.warn and RuntimeWarning that might also have disappeared already (but those are standard library / builtins, so hopefully not?)

In any case, you can verify this, too: - create a new environment using mamba create -n test python=3.11 numpy pandas packaging pooch netcdf4 and activate it - run pip install 'dask[array]' to install dask without distributed (this appears to make a difference for me, not sure if that's the same elsewhere) - editable-install xarray so we can easily switch between branches - run python -c 'import xarray as xr; ds = xr.tutorial.open_dataset("air_temperature", chunks={})'

This should print an error for main and shouldn't with this branch (confirmed on my local HPC).

the issue is that this doesn't occur on normal garbage collection but only on interpreter shutdown. So really, I don't think we have any way to test this using pytest as that itself is written in python (unless of course we can make use of sub-interpreters, but that might be more trouble than it's worth).

Maybe you can add a test that creates a cachingFilemanager object, then deletes it, then run gc.collect() and check somehow if it works? But no idea how pytest interferes with this or how you can ensure that there are no references to the module anymore?

OK, I won't notify you again about this release, but will get in touch when a new version is available. If you'd rather skip all updates until the next major or minor version, let me know by commenting @dependabot ignore this major version or @dependabot ignore this minor version. You can also ignore all major, minor, or patch releases for a dependency by adding an ignore condition with the desired update_types to your config file.

If you change your mind, just re-open this PR and I'll resolve any conflicts on it.

I agree that this seems very hard to test!

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

Hmmm ndim is in the array api so potentially we could just update the test.

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.

Please ask over at the cfgrib repo. But it does look like a bad environment / bad install.

The cfgrib notebook in the documentation is broken. I guess it's related to this PR. See: https://docs.xarray.dev/en/stable/examples/ERA5-GRIB-example.html

Same problem with https://github.com/pydata/xarray/issues/7841

The cfgrib notebook in the documentation is broken. I guess it's related to this PR. See: https://docs.xarray.dev/en/stable/examples/ERA5-GRIB-example.html

@tomwhite I've added tests to check the backend code for vlen string dtype metadadata. Also had to add specific check for the h5py vlen string metadata. I think we've covered everything for the proposed change to allow empty vlen strings dtype metadata.

I'm looking at the mypy error and do not have the slightest clue what and where to change. Any help appreciated.

this is #7879 (and thus probably #7079). I suspect our locks are not working properly, but in any case we really should try to fix this.

Just running import cfgrib returns the error below. But this works for you?

Is this an issue that cannot be solved?

Thanks again for all your help!

```

ModuleNotFoundError Traceback (most recent call last) Cell In[32], line 1 ----> 1 import cfgrib

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib__init__.py:20 18 # cfgrib core API depends on the ECMWF ecCodes C-library only 19 from .abc import Field, Fieldset, Index, MappingFieldset ---> 20 from .cfmessage import COMPUTED_KEYS 21 from .dataset import ( 22 Dataset, 23 DatasetBuildError, (...) 27 open_from_index, 28 ) 29 from .messages import FieldsetIndex, FileStream, Message

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib\cfmessage.py:29 26 import attr 27 import numpy as np ---> 29 from . import abc, messages 31 LOG = logging.getLogger(name) 33 # taken from eccodes stepUnits.table

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib\messages.py:28 25 import typing as T 27 import attr ---> 28 import eccodes # type: ignore 29 import numpy as np 31 from . import abc

File ~\Anaconda3\envs\doom_test\lib\site-packages\eccodes__init__.py:13 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 13 from .eccodes import * # noqa 14 from .highlevel import *

File ~\Anaconda3\envs\doom_test\lib\site-packages\eccodes\eccodes.py:12 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 12 from gribapi import ( 13 CODES_PRODUCT_ANY, 14 CODES_PRODUCT_BUFR, 15 CODES_PRODUCT_GRIB, 16 CODES_PRODUCT_GTS, 17 CODES_PRODUCT_METAR, 18 ) 19 from gribapi import GRIB_CHECK as CODES_CHECK 20 from gribapi import GRIB_MISSING_DOUBLE as CODES_MISSING_DOUBLE

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi__init__.py:13 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 13 from .gribapi import * # noqa 14 from .gribapi import version, lib 16 # The minimum recommended version for the ecCodes package

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\gribapi.py:34 30 from functools import wraps 32 import numpy as np ---> 34 from gribapi.errors import GribInternalError 36 from . import errors 37 from .bindings import ENC

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\errors.py:16 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 9 # does it submit to any jurisdiction. 10 # 12 """ 13 Exception class hierarchy 14 """ ---> 16 from .bindings import ENC, ffi, lib 19 class GribInternalError(Exception): 20 """ 21 @brief Wrap errors coming from the C API in a Python exception object. 22 23 Base class for all exceptions 24 """

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\bindings.py:40 37 # default encoding for ecCodes strings 38 ENC = "ascii" ---> 40 ffi = cffi.FFI() 41 CDEF = pkgutil.get_data(name, "grib_api.h") 42 CDEF += pkgutil.get_data(name, "eccodes.h")

File ~\Anaconda3\envs\doom_test\lib\site-packages\cffi\api.py:48, in FFI.init(self, backend) 42 """Create an FFI instance. The 'backend' argument is used to 43 select a non-default backend, mostly for tests. 44 """ 45 if backend is None: 46 # You need PyPy (>= 2.0 beta), or a CPython (>= 2.6) with 47 # _cffi_backend.so compiled. ---> 48 import _cffi_backend as backend 49 from . import version 50 if backend.version != version: 51 # bad version! Try to be as explicit as possible.

ModuleNotFoundError: No module named '_cffi_backend' ```

No module named '_cffi_backend'

Does simply import cfgrib work for you? I suspect it doesn't, which would explain the issue. It's unfortunate that the error is rewritten to "unknown engine", but I'm not sure how we would detect that it's a dependency that fails.

The version of cfgrib i have is: cfgrib=0.9.10.4

I tried to update my own environment as you advised and have this error:

```

ModuleNotFoundError Traceback (most recent call last) File ~\Anaconda3\envs\doom_test\lib\site-packages\xarray\tutorial.py:151, in open_dataset(name, cache, cache_dir, engine, **kws) 150 try: --> 151 import cfgrib # noqa 152 except ImportError as e:

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib__init__.py:20 19 from .abc import Field, Fieldset, Index, MappingFieldset ---> 20 from .cfmessage import COMPUTED_KEYS 21 from .dataset import ( 22 Dataset, 23 DatasetBuildError, (...) 27 open_from_index, 28 )

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib\cfmessage.py:29 27 import numpy as np ---> 29 from . import abc, messages 31 LOG = logging.getLogger(name)

File ~\Anaconda3\envs\doom_test\lib\site-packages\cfgrib\messages.py:28 27 import attr ---> 28 import eccodes # type: ignore 29 import numpy as np

File ~\Anaconda3\envs\doom_test\lib\site-packages\eccodes__init__.py:13 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 13 from .eccodes import * # noqa 14 from .highlevel import *

File ~\Anaconda3\envs\doom_test\lib\site-packages\eccodes\eccodes.py:12 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 12 from gribapi import ( 13 CODES_PRODUCT_ANY, 14 CODES_PRODUCT_BUFR, 15 CODES_PRODUCT_GRIB, 16 CODES_PRODUCT_GTS, 17 CODES_PRODUCT_METAR, 18 ) 19 from gribapi import GRIB_CHECK as CODES_CHECK

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi__init__.py:13 1 # 2 # (C) Copyright 2017- ECMWF. 3 # (...) 10 # 11 # ---> 13 from .gribapi import * # noqa 14 from .gribapi import version, lib

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\gribapi.py:34 32 import numpy as np ---> 34 from gribapi.errors import GribInternalError 36 from . import errors

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\errors.py:16 12 """ 13 Exception class hierarchy 14 """ ---> 16 from .bindings import ENC, ffi, lib 19 class GribInternalError(Exception):

File ~\Anaconda3\envs\doom_test\lib\site-packages\gribapi\bindings.py:40 38 ENC = "ascii" ---> 40 ffi = cffi.FFI() 41 CDEF = pkgutil.get_data(name, "grib_api.h")

File ~\Anaconda3\envs\doom_test\lib\site-packages\cffi\api.py:48, in FFI.init(self, backend) 45 if backend is None: 46 # You need PyPy (>= 2.0 beta), or a CPython (>= 2.6) with 47 # _cffi_backend.so compiled. ---> 48 import _cffi_backend as backend 49 from . import version

ModuleNotFoundError: No module named '_cffi_backend'

The above exception was the direct cause of the following exception:

ImportError Traceback (most recent call last) Cell In[30], line 2 1 import xarray as xr ----> 2 xr.tutorial.open_dataset("era5-2mt-2019-03-uk.grib")

File ~\Anaconda3\envs\doom_test\lib\site-packages\xarray\tutorial.py:153, in open_dataset(name, cache, cache_dir, engine, **kws) 151 import cfgrib # noqa 152 except ImportError as e: --> 153 raise ImportError( 154 "Reading this tutorial dataset requires the cfgrib package." 155 ) from e 157 url = f"{base_url}/raw/{version}/{path.name}" 159 # retrieve the file

ImportError: Reading this tutorial dataset requires the cfgrib package. ```

Thanks @mgunyho !

Thanks you @keewis and @Illviljan! I made comment to deprecate cdms2 in xarray in another issue/PR last year and didn't get around to it. I linked this PR in an xCDAT discussion post here.

No, this should still work: sh conda create -n test -c conda-forge xarray ipython python=3.11 cfgrib pooch conda activate test ipython ```python import xarray as xr

xr.tutorial.open_dataset("era5-2mt-2019-03-uk.grib") `` We somewhat recently dropped the builtincfgribengine in favor of the one provided by thecfgribpackage (and that is also the reason why the example in the docs fails:cfgrib` is not installed into the docs environment anymore, which is definitely an oversight).

Which version of cfgrib do you have? For reference, in the environment built with the command above I have cfgrib=0.9.10.3

Thanks for getting back to me!

cfgrib is installed. I believe it might be a case that grib files are no longer readable in this way that I am following from: https://docs.xarray.dev/en/stable/examples/ERA5-GRIB-example.html

As there are error messages on this example page.

great, thanks for the confirmation!

as stated by the exception, the cfgrib engine is unknown, which usually means you're missing the cfgrib package (or this is a environment issue). If you did indeed install it, can you post the output of either conda list or pip list (in case you're not using conda)?

Thanks for opening your first issue here at xarray! Be sure to follow the issue template! If you have an idea for a solution, we would really welcome a Pull Request with proposed changes. See the Contributing Guide for more. It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better. Thank you!

I am closing this. Feel free to re-open/ or open a new issue.

I think this is intended (though certainly not very easy to get right): see the second part of the warning in the __del__ documentation.

You are right, that warning is exactly what is causing the issues.

I would have thought that the global (or module level here) variable/function aquire should have at least one reference until after the deletion of the object.

I think this is intended (though certainly not very easy to get right): see the second part of the warning in the __del__ documentation.

That's quite a weird bug. I would have thought that the global (or module level here) variable/function aquire should have at least one reference until after the deletion of the object. Is that a bug in pythons garbage collection?

Or does the garbage collection already start when calling del and does not wait for the completion of the __del__ method?

If I'm reading the different issues correctly, that means this is a duplicate of #7079

There are similar segfaults in an xncml PR: https://github.com/xarray-contrib/xncml/pull/48

Googling around suggest it is related to netCDF not being thread-safe and recent python-netcdf4 releasing the GIL.

closing, since anything still missing should be feature requests for xncml

I think this has been fixed by xncml and/or kerchunk.

One other tricky edge case that occurs to me is one where an extreme fill value (e.g. 1e30) is used for floating point fields. If we decode the times first, it might appear that the dates cannot be represented as nanosecond-precision values, but in reality they would be. We may need to think more about how to handle this edge case in addition to #7817.

Actually, it goes away with pip install jinja2. We don't use jinja2 at all, so either this is some kind of weird effect on garbage collection (a timing issue?), or dask is doing something differently as soon as jinja2 is available.

Edit: most likely this is a timing issue... the offending line tries to make use of the internal acquire function, which I think at that point has already been destroyed. To fix that, I think we need to somehow store a reference on the file manager?

I can reproduce this locally: - download and unpack the files from https://github.com/pydata/xarray/issues/7814#issuecomment-1535168128 - use mamba create -n test python=3.11 xarray netcdf4 to create the environment (note: no dask) - use pip install "dask[array]" to install dask (does not pull distributed like the package from conda-forge) - put the code into a script and execute it

For reference, the full traceback is: pytb Exception ignored in: <function CachingFileManager.__del__ at 0x7fbdb237b430> Traceback (most recent call last): File "/home/jmagin/.local/opt/mambaforge/envs/test/lib/python3.9/site-packages/xarray/backends/file_manager.py", line 246, in __del__ TypeError: 'NoneType' object is not callable

As far as I can tell, this means we're using something from distributed with a broken fallback, since the error goes away as soon as I install distributed.

Congratulations on completing your first pull request! Welcome to Xarray! We are proud of you, and hope to see you again!

cos is a float operation so I would lean towards using a isclose-check: xr.testing.assert_allclose(a + 1, np.cos(a)).

Same issue here. I installed xarray with conda/mamba (not a dev install).

``` INSTALLED VERSIONS

commit: None python: 3.11.3 | packaged by conda-forge | (main, Apr 6 2023, 08:57:19) [GCC 11.3.0] python-bits: 64 OS: Linux OS-release: 3.10.0-1160.42.2.el7.x86_64 machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: ('en_US', 'UTF-8') libhdf5: 1.14.0 libnetcdf: 4.9.2

xarray: 2023.4.2 pandas: 2.0.1 numpy: 1.24.3 scipy: 1.10.1 netCDF4: 1.6.3 h5netcdf: None h5py: None zarr: 2.14.2 dask: 2023.4.1 distributed: None pip: 23.1.2 IPython: 8.13.1 ```

Edit: downgrading to 2023.4.0 solved the issue.

how did you set up your environment? This works for me: sh mamba create -n test python=3.11 xarray dask netcdf4 pooch ipython mamba activate test ipython python xr.tutorial.open_dataset("rasm", chunks={})

Interestingly enough, though, is that you should only see this with xarray=2023.5.0, while your environment claims to have xarray=2023.4.2. It seems there is something wrong with your environment?

Yes float64 should cause less imprecision. You can convert using astype:

```python import numpy as np import xarray as xr

da = xr.DataArray(np.array([1, 2], dtype=np.float32))

da = da.astype(float) ```

As for the other problems I think you are better of asking the people over at rioxarray. However, you should first gather all the steps you did to convert the data as code. This way it is easier to see what you are actually doing.

Thanks for getting back. I did post in rioxarray and yet, the last step I mentioned isn't successful there too. I'll post the code maybe 8hrs from here(can reach out to my sys then). Thanks for all the helpful suggestions so far. Really helpful.

Yes float64 should cause less imprecision. You can convert using astype:

```python import numpy as np import xarray as xr

da = xr.DataArray(np.array([1, 2], dtype=np.float32))

da = da.astype(float) ```

As for the other problems I think you are better of asking the people over at rioxarray. However, you should first gather all the steps you did to convert the data as code. This way it is easier to see what you are actually doing.

xarray handles nan values and ignores them per default - so you don't need to remove them. For example:

```python import numpy as np import xarray as xr

da = xr.DataArray([1, 2, 3, np.nan]) da.mean() ```

This is really helpful as I didn't know this before.

If you have precision problems - that might be because you have float32 values.

Which format would not cause the issue in that case float 64? If yes, can we manually convert?

I don't know what goes wrong with your lon values - that is an issue in the reprojection. You could convert them to 0...360 by using

python lon_dim = "x" new_lon = np.mod(da[lon_dim], 360) da = da.assign_coords(**{lon_dim: new_lon}) da.reindex(**{lon_dim : np.sort(da[lon_dim])})

Yeah. I have done the 180 to 360 deg conversions before. But the issue is more of with rioxarray reprojection I feel The internet data is in meters, as I wanted in degrees/lat-lon format, I converted the data from polar stereographic to wgs84. This converted the datas coordinates to degrees, latitudes are perfect. But longitude are arranged to -180 to +180 instead of 160E to 199W. I as well tried wrapping longitude to 0-360, but it should technically fall in 160-200 range while the long show all 0-360 and stretch throughout, which isn't right.

So, converting the existing gridded data (in meters) to lat-lon projection without affecting the resolution and without nan is my ultimate aim/objective. I successfully converted data to lat-lon and clipped to region but, it drastically changed the resolution like around 20 times maybe. Preserving the resolution is very imp for my work. So, that's the issue with longitudes

Thanks for your time if you went through this.

Thank you for opening this pull request! It may take us a few days to respond here, so thank you for being patient. If you have questions, some answers may be found in our contributing guidelines.

I agree, this change should be fine.

Thank you @TomNicholas, that is encouraging to hear. I will wait for @keewis to respond before filing a PR.

FWIW, I have tested the modification I suggest in my fork of xarray, and it works well for our purposes. It just generalizes the exception catch.

xarray handles nan values and ignores them per default - so you don't need to remove them. For example: ```python import numpy as np import xarray as xr

da = xr.DataArray([1, 2, 3, np.nan]) da.mean() `` If you have precision problems - that might be because you havefloat32` values.

I don't know what goes wrong with your lon values - that is an issue in the reprojection. You could convert them to 0...360 by using

```python

lon_dim = "x" new_lon = np.mod(da[lon_dim], 360) da = da.assign_coords({lon_dim: new_lon}) da.reindex({lon_dim : np.sort(da[lon_dim])})

```

I want to add a +1 to disable it by default. It's pretty common to be using float32 precision arrays. I have a rolling mean operation early on in some code and the errors balloon over time in subsequent processes. This was a super obscure bug to track down as well.

Well, that does makes sense. I want to calculate anomalies along x-y grids and I'm guessing the nan values are interfering with the results. Also, I have another question which isn't regarding Nan's. if it is right here, I may proceed. (else tag/link to other places/forums relevant). Assuming you must be knowing: I reprojected my nc file from meters to degrees Now, although the projection is right, the values of longitude aren't. python x (x) float64 -179.2 -177.7 ... 177.7 179.2 array([-179.217367, -177.65215 , -176.086933, -174.521715, -172.956498, -171.391281, -169.826063, -168.260846, -166.695629, -165.130412, -163.565194, -161.999977, -160.43476 , -158.869542, 163.565218, 165.130436, 166.695653, 168.26087 , 169.826088, 171.391305, 172.956522, 174.521739, 176.086957, 177.652174, 179.217391]) This is not how it is supposed to be: It should fall with 160-200 longitudes (post wrapping 360)

Is there a way xarray can sort this automatically or do I need to manually reset the cordinates?

Thanks @ghiggi for your comment.

The problem is we have at least two contradicting user requests here, see #7328 and #7862.

I'm sure there is a solution to accommodate both sides.

Thanks for raising this @vhaasteren ! We want to do what we can to support users from all fields of science :)

I would be okay with that change (especially as it's not really special-casing pint-pulsar, so much as generalizing an existing error-catching mechanism), but would defer to the opinion of @keewis on this.

Solution for those who just found this issue:

Just re-install xarray. pip install -e . is sufficient. Re-installing any way through pip/conda should register the dask chunkmanager entrypoint.

@Illviljan I brought this up in the xarray team call today and we decided that since this only affects people who have previously cloned the xarray repository, are using a development install, and then updated by pulling changes from main; this problem only affects maybe ~10-20 people worldwide, all of whom are developers who are equipped to quickly solve it.

I'm going to add a note into the what's new entry for this version now - if you think we need to do more then let me know.

EDIT: I added a note to whatsnew in https://github.com/pydata/xarray/commit/69445c62953958488a6b35fafd8b9cfd6c0374a5, and updated the release notes.

If you're curious what happened, we had the same problem: https://github.com/SciTools/iris/issues/5280#issuecomment-1525802077

Just wish I'd spotted this sooner but it's quite hard to follow two organisations' repos 😆

Dask array with dtype object can contain whatever python object (i.e. I saw examples of geometry and matplotlib collections within dask arrays with object dtype). As a consequence, dask do not try the conversion to i.e. str to estimate the array size, since there is no clean way AFAIK to attach an attribute to dtype suggesting that the object is actually a string.

With your PR, the dtype is not anymore object when creating the dask.array and this solves the issue I guess. Did I overlooked something?

Do you know where the in-place modification is happening? We could just copy there and fix this particular issue.

Not sure, but I'll take a look!

We dropped Python 3.8 support prior to the Pandas 2 release and have no plans to backport support at this time.

xref: #7765

Thanks for the answer!

So it looks like the changes here with the fix in my branch will get your issue resolved @tomwhite, right?

Yes - thanks!

I'm a bit worried, that this might break other users workflows, if they depend on the current conversion to floating point for some reason.

The floating point default is preserved if you do e.g. xr.Dataset({"a": np.array([], dtype=object)}). The change here will only convert to string if there is extra metadata present that says it is a string.

We dropped Python 3.8 support prior to the Pandas 2 release and have no plans to backport support at this time.

xref: #7765

Thanks for trying. I can't think of any downsides for the netcdf4-fix, as it just adds the needed metadata to the object-dtype. But you never know, so it would be good to get another set of eyes on it.

So it looks like the changes here with the fix in my branch will get your issue resolved @tomwhite, right?

I'm a bit worried, that this might break other users workflows, if they depend on the current conversion to floating point for some reason. Also other backends might rely on this feature. Especially because this has been there since the early days when xarray was known as xray.

@dcherian What would be the way to go here?

There is also a somehow contradicting issue in #7868.

Thanks for opening your first issue here at xarray! Be sure to follow the issue template! If you have an idea for a solution, we would really welcome a Pull Request with proposed changes. See the Contributing Guide for more. It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better. Thank you!

Could you verify the above example, please?

The code looks fine, and I get the same result when I run it with this PR.

Your fix in https://github.com/kmuehlbauer/xarray/tree/preserve-vlen-string-dtype changes the metadata so it is correctly preserved as metadata: {'element_type': <class 'str'>}.

I feel less qualified to evaluate the impact of the netcdf4 fix.

My main question here is, why is dask not trying to retrieve the object types from dtype.metadata? Or does it and fail for some reason?.

@tomwhite I've put a commit with changes to zarr/netcdf4-backends which should preserve the dtype metadata here: https://github.com/kmuehlbauer/xarray/tree/preserve-vlen-string-dtype.

I'm not really sure if that is the right location, but as it was already present that location at netcdf4-backend I think it will do.

Do you know where the in-place modification is happening? We could just copy there and fix this particular issue.

@tomwhite Special casing on netcdf4 backend should be possible, too.

But it might need fixing at zarr backend, too:

python ds = xr.Dataset({"a": np.array([], dtype=xr.coding.strings.create_vlen_dtype(str))}) print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ds.to_zarr("a.zarr") print("\n### Loading ###") with xr.open_dataset("a.zarr", engine="zarr") as ds: print(f"dtype: {ds['a'].dtype}") print(f"metadata: {ds['a'].dtype.metadata}") ```python dtype: object metadata: {'element_type': <class 'str'>}

Loading

dtype: object metadata: None ```

Could you verify the above example, please? I'm relatively new to zarr :grimacing:

Thanks for taking a look @kmuehlbauer and for the useful example code. I hadn't considered the netcdf cases, so thanks for pointing those out.

Engine netcdf4 does not roundtrip here, losing the dtype metadata information. There is special casing for h5netcdf backend, though.

Could netcdf4 do the same special-casing as h5netcdf?

You can always reach out to the creator of the original PR and comment in the PR.

But it looks like this particular PR was reaching a dead end and should be completely rewritten. But anyway the reviewers left helpful remarks on how to proceed.

This is currently stuck waiting until the problems with the last netcdf-c versions are fixed in a new release. See the issues (https://github.com/pydata/xarray/issues/7388).

When they are fixed I will write the tests If I have time. But of course any help and suggestions are welcomed.

Thanks for opening your first issue here at xarray! Be sure to follow the issue template! If you have an idea for a solution, we would really welcome a Pull Request with proposed changes. See the Contributing Guide for more. It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better. Thank you!

@markelg Thanks a lot for adding this! Do you have time to finalize it in the near future? If not, I could also take a look at the tests if you like.

Yes but there are less - so as mentioned it removes all columns/ rows with only nans, if there is at least one non-nan value the row is kept.

What is the reason that you want to get rid of the nan values?

@ghiggi Glad it works, but we still have to check if that is the correct location for the fix, as it's not CF specific.