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

```

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.

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!

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.

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' ```

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 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.

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!

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

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.

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 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?

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?

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!

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!

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.

Thanks @kmuehlbauer ! https://github.com/pydata/xarray/pull/7869 solve the issues !

Summarizing: - With #7869, netCDF4 with NC_STRING variable arrays are now read into xarray as Unicode dtype (instead of object) - As a consequence dask can estimate the array's size and xr.open_dataset(fpath, chunks="auto") does not raise anymore the NotImplementedError. - NC_CHAR variable arrays continue to be read into xarray as fixed-length byte-string dtype. Maybe something more could be done to deserialize also NC_CHAR to Unicode. However, this might cause some backward incompatibilities and might be better to address this in a separate PR.

Thanks again @kmuehlbauer for having resolved the problem in less than 2 hours :1st_place_medal:

Can you try notnull instead of isnull - I often get the boolean array wrong in where:

python da = ds['z'] da = da.where(da.notnull(), drop=True)

Yes, I did.

As we can see, the nan values are not completely gone

Thanks alot for responding, but,

python da = ds['z'] da = da.where(da.isnull(), drop=True) is for pciking nan values? Because the data array has all 'nan' values

when I plot: I get,

I need to use data that has no empty cells for further analysis.

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!

Another cause for this issue could be an unclosed netcdf file. Specifically if you did something like:

my_file = Dataset('my_file.nc','w')

and did not close it with

my_file.close()

The ncfile is still open and the HDF library still has a lock on it.

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

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!

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!

cc @spencerkclark @znicholls

Legend thanks

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.

@martinfleis already volunteered to do it using the github action we're working on in the scientific-python space.

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!

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!

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

I'm sorry, but just to reiterate, "x" is referred to as a dimension several times in the document. It is then, one time, mistakenly referred to as a coordinate. This pull request, which should be accepted, is to make the document consistent in how it refers to "x".

+1 https://stackoverflow.com/questions/76239626/xarray-preprocess-for-nearest-lat-lon-non-nan-variable

Thanks, yes. The remote access and local access are different with xarray v0.20.2. remote:

local access:

Yeah, as a total new comer, this explanation seems not only unintuitive but also like we're fighting with common terminology used in math and everywhere else.

When you define the coordinates,

data = xr.DataArray(np.random.randn(2,3)*10, dims=("x", "y"), coords={"x": [10, 20]})

It looks like you are defining a specific coordinate system [10,20] on a specific dimension "x". The key of the object-value pair set in coords (ie. "x") MUST refer to a dimension, because "x" is a dimension that we are defining a coordinate system for. As we don't define any coordinate system on "y", asking xarray for data.coords["y"] gives us a array([0,1,2]) which we'd consider to be an intuitive default coordinate system.

Intuitively, we think: An xarray contains dimensions along which specific coordinates can be defined: xarray.dimension.coordinate

In the attribute definition, we say:

data.x.attrs["units"] = "x units" which feels like we are defining an attribute on the entire dimension. As opposed to data.x.10 (I realize the problem here) which would be definiting an attribute on a single coordinate within a dimension.

I'm done, thank you for your time. I will continue reading and learning about xarray.

Well, it's not consistent within this particular document, which is reasonably the first anyone reads. The dimension, as its refered to, in this example is called "x", we are assigning an attribute to "x", the coordinates of "x" are [10,20], we are not assigning an attribute to a coordinate within the dimension "x".

v 2022.11.0 v0.20.2

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.

Since folks have asked questions since this was closed (and I had a similar issue), this ValueError exception can be raised in a misleading way if the file being specified does not exist. The exception implies that it is a file format error, but it may just not exist.

An empty file will raise the same ValueError.

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!

Thank you all so much for the feedback and resources! I agree (1) testing the limits of xArray's API compatibility with sparse and (2) developing some documentation for what is/isn't supported are great places to start, so I'll get on that while I think about the other I/O issues (serialization, etc.)

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!

Thanks for the replies. Yes, that second suggestion sounds good @kmuehlbauer!

I realise it's not practical to add specific checks / messages for all engines, so something like this that links to a webpage that describes potential solutions seems like an excellent compromise. Your earlier solution (rephasing the error) I think would not help, however, as it still doesn't show users what the actual missing package is rioxarray vs rasterio.

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!

:heavy_plus_sign: 1

``` ❯ pip freeze |grep -Ei "xarray|numpy|netcdf" netCDF4==1.6.3 numpy==1.23.3 xarray==2023.4.2

❯ python -c "import xarray;import warnings;warnings.filterwarnings('error');import netCDF4"
Traceback (most recent call last): File "<string>", line 1, in <module> File "/home/z3538708/projects/gitlab.com/sqc-eng/se/core/.venv-core/lib/python3.10/site-packages/netCDF4/init.py", line 3, in <module> from ._netCDF4 import * File "src/netCDF4/_netCDF4.pyx", line 1, in init netCDF4._netCDF4 RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility. Expected 16 from C header, got 96 from PyObject

❯ python -c "import netCDF4;import xarray;import warnings;warnings.filterwarnings('error');import netCDF4" ❯ `` No errors when the order isnetCD4` import, before import xarray

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.

Thanks Kai, I have attached the output of xr.show_versions() and of pip list.

Does this happen only with these special files, or do you experience this every time?

I tried with two netcdf files from https://www.northwestknowledge.net/metdata/data/ and got the same TypeError: 'NoneType' object is not callable

Does the problem persists when specifying engine="netcdf4" or engine="h5netcdf" in the call to open_mfdataset?

Only engine="netcdf4" works as a parameter but I get the same error, h5netcdf is not recognized as an engine.

Does this also happen, if you open the files one-by-one (with xr.open_dataset) and combine the Datasets with xr.concat?

Yes, opening files with xr.open_dataset and combining them with xr.concat gives the same error.

Please note the following that I just realized, the resulting combined netCDF file is actually created even though the errors are displayed. Also, as I mentioned in my first message, if I run the commands line by line in the python console no error message is displayed, this happens only when the code is ran as a script.

Paul

pip_list.txt xr.show_versions.txt

Might be related to https://github.com/dask/distributed/issues/6402

Hi,

I forgot to rebuild the package after removing the BACKEND_... line. With only the line in pyproject.toml it works as it should! My mistake. Thanks for the patience.

Regards, Gaute

Hi,

Yes, I tried that, but I then got the same error as if I kept that line in the old format. I'll do a few tests and post the proper error here.

Gaute

If I do not manually add the backend to the array, but only have this line in https://github.com/gauteh/hidefix/blob/main/pyproject.toml#L29:

[project.entry-points."xarray.backends"] hidefix = "hidefix.xarray:HidefixBackendEntrypoint"

which is only what is supported by pyproject.toml/maturin I get an error where xarray expected a tuple and cannot parse the entrypoint, not just the adderss to the entrypoint - as it used to be (back in January at least).

As well as libgdal-dev I noticed that g++ was also required for the Fiona wheel build.

Deepak, thanks for replying, I have attached a couple of my netCDF files. Hopefully these will help for reproducing the issue.

Archive.zip

That's a totally valid scope limitation for the sparse package, and I understand the motivation.

I'm just saying that the principle of least astonishment is not being followed: the user cannot at the moment read either the xarray or sparse docs and know which portions of the xarray API will work when giving …, sparse=True, and which instead require a deliberate choice to densify, or see examples of how best to mix the two. It would be helpful to clarify—that's all.

@jbbutler please also see this comment et seq. https://github.com/pydata/sparse/issues/1#issuecomment-792342987 and related pydata/sparse#438.

To add to @rabernat's point about sparse support being "not well documented", I suspect (but don't know, as I'm just a user of xarray, not a developer) that it's also not thoroughly tested. I expected to be able to use e.g. DataArray.cumprod when the underlying data was sparse, but could not.

IMHO, I/O to/from sparse-backed objects is less valuable if only a small subset of xarray functionality is available on those objects. Perhaps explicitly testing/confirming which parts of the API do/do not currently work with sparse would support the improvements to the docs that Ryan mentioned, and reveal the work remaining to provide full(er) support.

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!

Hi all! As part of a research project, I'm looking to contribute to xArray's sparse capabilities, with an emphasis on sparse support for use-cases in the geosciences. I'm wondering if anyone in the geosciences (or adjacent disciplines!) has encountered problems with xArray's current level of sparse support, and what kinds of improvements they'd like to see to address those issues. From playing around, it seems the current strategy of backing DataArrays with COO sparse arrays takes care of a lot of use cases, but I have the following ideas that may (or may not) be useful to implement further:

• Functions/methods to convert from geopandas GeoDataframes of vector data to rasterized, potentially sparse ndarrays in an xArray Dataset/DataArray (reverse direction too); this is related to the issue of converting from sparse arrays back to multi-indexed pandas objects at the top of this issue (which I believe has yet to be solved)
• Loading sparse data from a netcdf() file directly into a Dataset/Array backed by sparse ndarray(s) (seems like the only way to get sparse backings is to either unstack or call '.from_dataframe()/series()' with the sparse flag set to True?)
• Support for other sparse array conventions (for ex, GCXS in the sparse package for better memory efficiency; I can't find any improvements to make on the current COO backing in terms of supported arithmetic operations, merges/joins, etc.)

I'd appreciate any feedback on these ideas, as well as any other things that would be nice to have implemented!

@kmuehlbauer - genius! Yes. That pull request should fix this issue exactly! And it explains why I see this issue and you don't - with undefined behavior anything can happen. Since we are on different OSes, our systems behave differently.

I just double checked with pandas and this fix will do the right thing: python import pandas as pd print(pd.to_timedelta([np.nan, 0],"ns") + np.datetime64('1970-01-01')) DatetimeIndex(['NaT', '1970-01-01'], dtype='datetime64[ns]', freq=None) I see that the pull request with the fix has been sitting since December of last year. Is there some way to somehow get someone to look at that pull request who can merge it?

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!

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!

Ah hah! Well, I don't know why this is working for you @kmuehlbauer, but I can see why it is not working for me. I've been debugging through the code and it looks like the problem is the _decode_datetime_with_pandas function. For me, it's converting a float NaN into an integer, which results in a zero value.

It all starts in the open_zarr function, which by default sets the use_cftime parameter to None by default: https://github.com/pydata/xarray/blob/25d9a28e12141b9b5e4a79454eb76ddd2ee2bc4d/xarray/backends/zarr.py#L701-L817

There's a bunch of stuff that gets called, but eventually we get to the function decode_cf_datetime, which ironically (given the name) also takes this use_cftime parameter, which is still None. Because use_cftime is None, the function calls _decode_datetime_with_pandas:

https://github.com/pydata/xarray/blob/25d9a28e12141b9b5e4a79454eb76ddd2ee2bc4d/xarray/coding/times.py#L265-L289

and then, in _decode_datetime_with_pandas, the code casts a float NaN value to zero:

https://github.com/pydata/xarray/blob/979b99831f5d34d33120312a15dad3e6a0830f32/xarray/coding/times.py#L216-L262

In line 254, flat_num_dates is array([ nan, 1.6726176e+18]). After line 254, flat_nuM-dates_ns_int is array([ 0, 1672617600000000000]).

@kmuehlbauer - I ran https://github.com/pydata/xarray/issues/7790#issuecomment-1529894939 and I get an incorrect fill value:

```

Created with fill value 1900-01-01 <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02

Read back out of the zarr store with xarray <xarray.DataArray 'time' (time: 2)> array(['1970-01-01T00:00:00.000000000', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 1970-01-01 2023-01-02 {} {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -2208988800000000000, 'units': 'nanoseconds since 1970-01-01', 'calendar': 'proleptic_gregorian', 'dtype': dtype('int64')}

Read back out of the zarr store with zarr <zarr.core.Array '/time' (2,) int64 read-only> <zarr.attrs.Attributes object at 0x132802a50> [-2208988800000000000 1672617600000000000]

and here is my show_versions, since it may have changed because I've added some new libraries. It looks like my ipython version is slightly different, but I can't see how that would affect things. INSTALLED VERSIONS

commit: None python: 3.11.3 | packaged by conda-forge | (main, Apr 6 2023, 08:58:31) [Clang 14.0.6 ] python-bits: 64 OS: Darwin OS-release: 22.4.0 machine: arm64 processor: arm byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: ('en_US', 'UTF-8') libhdf5: None libnetcdf: None

xarray: 2023.4.2 pandas: 2.0.1 numpy: 1.24.3 scipy: None netCDF4: None pydap: None h5netcdf: None h5py: None Nio: None zarr: 2.14.2 cftime: None nc_time_axis: None PseudoNetCDF: None iris: None bottleneck: None dask: None distributed: None matplotlib: None cartopy: None seaborn: None numbagg: None fsspec: None cupy: None pint: None sparse: None flox: None numpy_groupies: None setuptools: 67.7.2 pip: 23.1.2 conda: None pytest: None mypy: None IPython: 8.13.1 sphinx: None ```

Oops! Yes. You are right. I had some cross-wording on the variable names. So I started a new notebook. Unfortunately, I think you may have also gotten some wires crossed? You set the time fill value to 1900-01-01, but then use NaT in the actual array?

Here is a fresh notebook with a stand-alone cell with everything that I think you were doing, but I'm not 100%. The fill value is still wrong when it gets read out, but it is at least different? The fill value is now set to the units for some reason. This seems like progress?

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

Create a time array with one fill value, NaT

time = np.array([np.datetime64("NaT", "ns"), '2023-01-02 00:00:00.00000000'], dtype='M8[ns]')

Create xarray with this fill value

xr_time_array = xr.DataArray(data=time,dims=['time'],name='time') xr_ds = xr.Dataset(dict(time=xr_time_array)) print("****") print("xarray created with NaT fill value") print("----------------------") print(xr_ds["time"])

Save as zarr

location_with_units = "xarray_and_units.zarr" encoding = { "time":{"_FillValue":np.datetime64("NaT","ns"),"dtype":np.int64,"units":"nanoseconds since 1970-01-01"} } xr_ds.to_zarr(location_with_units,mode="w",encoding=encoding)

Read it back out again

xr_read = xr.open_zarr(location_with_units) print("****") print("xarray created read with NaT fill value") print("----------------------") print(xr_read["time"]) print(xr_read["time"].attrs) print(xr_read["time"].encoding)

```

```

xarray created with NaT fill value

<xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02

xarray created read with NaT fill value

<xarray.DataArray 'time' (time: 2)> array(['1970-01-01T00:00:00.000000000', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 1970-01-01 2023-01-02 {} {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -9223372036854775808, 'units': 'nanoseconds since 1970-01-01', 'calendar': 'proleptic_gregorian', 'dtype': dtype('int64')}

```

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

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.

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.

I'm also relatively sure that if you are willing to put a floor on the version of Matplotlib you support get_window_extent will use it's internally cached renderer (and when we make it uniformly optional we also fixed the cache invalidation logic).

Where in the code is the time array being decoded? That seems to be where a lot of the issue is?

The zarr store does indeed use an integer in this case according to the .zmetadata file: { "metadata": { ".zattrs": {}, ".zgroup": { "zarr_format": 2 }, "time/.zarray": { "chunks": [ 2 ], "compressor": { "blocksize": 0, "clevel": 5, "cname": "lz4", "id": "blosc", "shuffle": 1 }, "dtype": "<i8", "fill_value": -25567, "filters": null, "order": "C", "shape": [ 2 ], "zarr_format": 2 }, "time/.zattrs": { "_ARRAY_DIMENSIONS": [ "time" ], "calendar": "proleptic_gregorian", "units": "days since 1900-01-01 00:00:00" } }, "zarr_consolidated_format": 1 } Once again the values in the zarr store are correct given the units, but xarray misreads the fill value for some reason: python z_new_fill = zarr.open('from_xarray_new_fill.zarr','r') z_new_fill["time"][:] array([ 0, 44926])

@kmuehlbauer - I think I'm not understanding what you are suggesting because the zarr store is still not being read correctly when I switch the fill value to a different date:

```python

Create a numpy array of type np.datetime64 with one fill value and one date

time_fill_value = np.datetime64("1900-01-01") time = np.array([time_fill_value,'2023-01-02'],dtype='M8[ns]')

Create a dataset with this one array

xr_time_array = xr.DataArray(data=time,dims=['time'],name='time') xr_ds = xr.Dataset(dict(time=xr_time_array))

print("******") print("Created with fill value 1900-01-01") print(xr_ds["time"])

Save the dataset to zarr

location_new_fill = "from_xarray_new_fill.zarr" encoding = { "time":{"_FillValue":time_fill_value,"dtype":np.int64} } xr_ds.to_zarr(location_new_fill,encoding=encoding)

xr_read = xr.open_zarr(location) print("******") print("Read back out of the zarr store with xarray") print(xr_read["time"]) print(xr_read["time"].encoding)

Created with fill value 1900-01-01 <xarray.DataArray 'time' (time: 2)> array(['1900-01-01T00:00:00.000000000', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 1900-01-01 2023-01-02

<xarray.DataArray 'time' (time: 2)> array(['2023-01-02T00:00:00.000000000', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 2023-01-02 2023-01-02 {'chunks': (2,), 'preferred_chunks': {'time': 2}, 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0), 'filters': None, '_FillValue': -9.223372036854776e+18, 'units': 'days since 2023-01-02 00:00:00', 'calendar': 'proleptic_gregorian', 'dtype': dtype('float64')} ```

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!

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

@jdldeauna how did you resolve this issue? I am seeing similar issues, but only encounter this when writing to zarr with a dask cluster.

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!

@kmuehlbauer It looks like a bug in the code if indeed tuples are meant to be treated the same as any sequence of data. These lines https://github.com/pydata/xarray/blob/0f4e99d036b0d6d76a3271e6191eacbc9922662f/xarray/core/variable.py#L259-L260 suggest that when a tuple is passed in, it is converted to a 0-dimension array of type object via https://github.com/pydata/xarray/blob/0f4e99d036b0d6d76a3271e6191eacbc9922662f/xarray/core/utils.py#L344-L348

Maybe removing the tuple type check and relying on this line https://github.com/pydata/xarray/blob/0f4e99d036b0d6d76a3271e6191eacbc9922662f/xarray/core/variable.py#L287-L288 is better?

Interestingly, xarray is also perfectly happy to read a numpy.datetime64 array out of a zarr store as long as the xarray metadata is present. xarray even helpfully creates an '_FillValue" attribute for the array so there is no confusion:

```

Create a zarr store directly with numpy.datetime64 type

location_zarr_direct = "from_zarr.zarr" root = zarr.open(location_zarr_direct,mode='w') z_time_array = root.create_dataset( "time",data=time,shape=time.shape,chunks=time.shape,dtype=time.dtype, fill_value=time_fill_value )

Add xarray metadata

z_time_array.attrs["_ARRAY_DIMENSIONS"] = ["time"] zarr.convenience.consolidate_metadata(location_zarr_direct)

Use xarray to read this data out

xr_read_from_zarr = xr.open_zarr(location_zarr_direct) print(xr_read_from_zarr["time"]) <xarray.DataArray 'time' (time: 2)> array([ 'NaT', '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] NaT 2023-01-02 Attributes: _FillValue: NaT ```

So I am extremely confused as to why xarray encodes time arrays so strangely when it creates the zarr store itself! (Hence https://github.com/pydata/xarray/discussions/7776)

Ah! Okay. I did not know about the .encoding option, which does indeed have the fill value. Thank you.

Interestingly, -9.223372036854776e+18 is just the float equivalent of numpy.datetime64('NaT'):

python float(np.datetime64('NaT').view('i8')) -9.223372036854776e+18

And I know this isn't an issue with zarr and NaT because I can create the zarr store directly with the zarr library and it's perfectly happy: ```python

Create a zarr store directly with numpy.datetime64 type

location_zarr_direct = "from_zarr.zarr" root = zarr.open(location_zarr_direct,mode='w') z_time_array = root.create_dataset( "time",data=time,shape=time.shape,chunks=time.shape,dtype=time.dtype, fill_value=time_fill_value ) zarr.convenience.consolidate_metadata(location_zarr_direct)

Read it back out again

read_zarr = zarr.open(location_zarr_direct,mode='r') print(read_zarr["time"][:]) [ 'NaT' '2023-01-02T00:00:00.000000000'] ```

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!

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.

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

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.

@dcherian Interesting! There should ideally be a way to set that because 32-64 seconds is way to long to wait before timing out in my opinion.

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!

Thank you @dcherian . I cannot reproduced this on main.

Do you need them in 4.9.1 then, or is updating to 4.9.2 an option?

The first point that you mentioned does not seem to be correct. Please see the below code (we took the sparse matrix ) and output: ``` import xarray as xa import numpy as np

def get_data(): lat_dim = 7210 lon_dim = 7440

lat = [0] * lat_dim
lon = [0] * lon_dim 
time = [0] * 5

nlats = lat_dim; nlons = lon_dim; ntimes = 5

var_1 = np.empty((ntimes,nlats,nlons))
var_2 = np.empty((ntimes,nlats,nlons))
var_3 = np.empty((ntimes,nlats,nlons))
var_4 = np.empty((ntimes,nlats,nlons))

data_arr = np.random.uniform(low=0,high=0,size=(ntimes,nlats,nlons))
data_arr[:,0,:] = 1
data_arr[:,:,1] = 1

var_1[:,:,:] = data_arr 
var_2[:,:,:] = data_arr 
var_3[:,:,:] = data_arr 
var_4[:,:,:] = data_arr

dataset = xa.Dataset( 
        data_vars = {
            'var_1': (('time','lat','lon'), var_1),
            'var_2': (('time','lat','lon'), var_2),
            'var_3': (('time','lat','lon'), var_3),
            'var_4': (('time','lat','lon'), var_4)},
        coords = {
            'lat': lat,
            'lon': lon,
            'time':time})

print(sum(v.size * v.dtype.itemsize for v in dataset.variables.values()))
print(dataset.nbytes)

if name == "main": get_data() ```

8582901240 8582901240 As we can observe here both nbytes and self.size * self.dtype.itemsize gives the same size.

And for the 2nd point can you share any solution for the nbytes for the netCDF or grib file as it takes too much memory and killed the process?

I open a new issue in netcdf github - https://github.com/Unidata/netcdf4-python/issues/1246

Sure @dcherian!

I add a simple test as you recommend, and all ubuntu tests crashed. ```python

def test_only_netCDF4_access(): dataset = Dataset( "https://ihthredds.ihcantabria.com/thredds/dodsC/Bathymetry/Global/Gebco_2020.nc" ) assert isinstance(dataset, Dataset)

````

And also I run this code in Googlecolab. Same problem... seems that something happend with netCDF4....

Related issue: https://github.com/pydata/xarray/issues/4122

So really, my question is: how do we support python scalars for libraries that only implement __array_namespace__, given that stopping to do so would be a major breaking change?

I was considering this question for SciPy (xref scipy#18286) this week, and I think I'm happy with this strategy: 1. Cast all "array-like" inputs like Python scalars, lists/sequences, and generators, to numpy.ndarray. 2. Require "same array type" input, forbid mixing numpy-cupy, numpy-pytorch, cupy-pytorch, etc. - this will raise an exception 3. As a result, cupy-pyscalar and pytorch-pyscalar will also raise an exception.

What that results in is an API that's backwards-compatible for numpy and array-like usage, and much stricter when using other array libraries. That strictness to me is a good thing, because: - that's what CuPy, PyTorch & co themselves do, and it works well there - it avoids the complexity raised by arbitrary mixing, which results in questions like the one raised in this issue. - in case you do need to use a scalar from within a function inside your own library, just convert it explicitly to the desired array type with xp.asarray(a_scalar) giving you a 0-D array of the correct type (add dtype=x.dtype to make sure dtypes match if that matters)

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!

Just confirmed that the fix works for the actual application that I was using (where N-NW etc were used as labels in binning wind directions), so all good in that respect.

Thanks, Justus, for expanding on this. It sounds to me the question is "how do we cast dtypes when multiple array libraries are participating in the same computation?" and I am not sure I am knowledgable enough to make any comment.

From the array API point of view, long long ago we decided that this is UB (undefined behavior), meaning it's completely up to each library to decide what to do. You can raise or come up with a special rule that you can make sense of.

It sounds like Xarray has some machinery to deal with this situation, but you'd rather prefer to not keep special-casing for a certain array library? Am I understanding it right?

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!

Thanks for the replies

So while xr.where(cond, x, y) is semantically, "where condition is true, x, else y", da.where(cond, x) is "where condition is true da, else x".

The latter feels quite unintuitive to me. Is the reason they're different only for the mask example you provided? Where NaN is returned as the default 'x' value?

Pandas seems to also look at the policies of its dependencies like Numpy: https://numpy.org/neps/nep-0029-deprecation_policy.html#support-table

See: https://github.com/pandas-dev/pandas/issues/52513