xCDAT support for Quadmesh

[gkb999]

Hi, I'm working with a data that it concatenated using mf_Dataset and exported.

My data looks like:

<xarray.Dataset>
Dimensions:  (x: 1264, y: 1328, time: 3287)
Coordinates:
  * x        (x) float64 -3.947e+06 -3.941e+06 ... 3.941e+06 3.947e+06
  * y        (y) float64 -3.947e+06 -3.941e+06 ... 4.341e+06 4.347e+06
  * time     (time) int32 0 1 2 3 4 5 6 7 ... 3280 3281 3282 3283 3284 3285 3286
Data variables:
    z        (time, y, x) float32 ...

After modifying the time variable

<xarray.Dataset>
Dimensions:    (x: 1264, y: 1328, time: 3287, bnds: 2)
Coordinates:
  * x          (x) float64 -3.947e+06 -3.941e+06 ... 3.941e+06 3.947e+06
  * y          (y) float64 -3.947e+06 -3.941e+06 ... 4.341e+06 4.347e+06
  * time       (time) object 2013-01-01 00:00:00 ... 2021-12-31 00:00:00
Dimensions without coordinates: bnds
Data variables:
    z          (time, y, x) float32 ...
    x_bnds     (x, bnds) float64 -3.95e+06 -3.944e+06 ... 3.944e+06 3.95e+06
    time_bnds  (time, bnds) object 2012-12-31 12:00:00 ... 2021-12-31 12:00:00'

Now,

1. The formats of x and y (are in meters), i want to convert them into lat-lon
2. As the data is really high resolution, I request any memory handling tips (As it is taking much time)
3. If you the the data after adding time variable, the data is very much in 'datetime64' format. Yet, I have issues computing anomalies and departures.
   For instance:

d_clm = sic.temporal.climatology("z", freq="day")

gives the follwoing error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[15], line 1
----> 1 d_clm = sic.temporal.climatology("z", freq="day")

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:496, in TemporalAccessor.climatology(self, data_var, freq, weighted, keep_weights, reference_period, season_config)
    358 """Returns a Dataset with the climatology of a data variable.
    359 
    360 Parameters
   (...)
    492 }
    493 """
    494 self._set_data_var_attrs(data_var)
--> 496 return self._averager(
    497     data_var,
    498     "climatology",
    499     freq,
    500     weighted,
    501     keep_weights,
    502     reference_period,
    503     season_config,
    504 )

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:772, in TemporalAccessor._averager(self, data_var, mode, freq, weighted, keep_weights, reference_period, season_config)
    770 # Get the data variable and the required time axis metadata.
    771 dv = _get_data_var(ds, data_var)
--> 772 time_bounds = ds.bounds.get_bounds("T", var_key=dv.name)
    774 if self._mode == "average":
    775     dv = self._average(dv, time_bounds)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\bounds.py:218, in BoundsAccessor.get_bounds(self, axis, var_key)
    215         bounds_keys = []
    217 if len(bounds_keys) == 0:
--> 218     raise KeyError(
    219         f"No bounds data variables were found for the '{axis}' axis. Make sure "
    220         "the dataset has bound data vars and their names match the 'bounds' "
    221         "attributes found on their related time coordinate variables. "
    222         "Alternatively, you can add bounds with `xcdat.add_missing_bounds` "
    223         "or `xcdat.add_bounds`."
    224     )
    226 bounds: Union[xr.Dataset, xr.DataArray] = self._dataset[
    227     bounds_keys if len(bounds_keys) > 1 else bounds_keys[0]
    228 ].copy()
    230 return bounds

KeyError: "No bounds data variables were found for the 'T' axis. Make sure the dataset has bound data vars and their names match the 'bounds' attributes found on their related time coordinate variables. Alternatively, you can add bounds with `xcdat.add_missing_bounds` or `xcdat.add_bounds`."

What is missing in the bounds here?
sic.bounds.get_bounds (where sic is xcdat dataset)

<bound method BoundsAccessor.get_bounds of <xcdat.bounds.BoundsAccessor object at 0x00000233DA29F890>>

I as well added

sic_T = sic.bounds.add_bounds('T')

now looks like:

xarray.Dataset
Dimensions:
x: 1264y: 1328time: 3287bnds: 2
Coordinates:
x
(x)
float64
-3.947e+06 -3.941e+06 ... 3.947e+06
y
(y)
float64
-3.947e+06 -3.941e+06 ... 4.347e+06
time
(time)
datetime64[ns]
2013-01-01 ... 2021-12-31
Data variables:
z
(time, y, x)
float32
...
x_bnds
(x, bnds)
float64
-3.95e+06 -3.944e+06 ... 3.95e+06
time_bnds
(time, bnds)
datetime64[ns]
2012-12-31T12:00:00 ... 2021-12-...
Indexes: (3)
Attributes: (0)

Yet, the departures and climatology doesn't seems to work.
Error:

---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In[22], line 1
----> 1 d_clm = sic_T.temporal.climatology("z", freq="day")

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:496, in TemporalAccessor.climatology(self, data_var, freq, weighted, keep_weights, reference_period, season_config)
    358 """Returns a Dataset with the climatology of a data variable.
    359 
    360 Parameters
   (...)
    492 }
    493 """
    494 self._set_data_var_attrs(data_var)
--> 496 return self._averager(
    497     data_var,
    498     "climatology",
    499     freq,
    500     weighted,
    501     keep_weights,
    502     reference_period,
    503     season_config,
    504 )

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:777, in TemporalAccessor._averager(self, data_var, mode, freq, weighted, keep_weights, reference_period, season_config)
    775     dv = self._average(dv, time_bounds)
    776 elif self._mode in ["group_average", "climatology", "departures"]:
--> 777     dv = self._group_average(dv, time_bounds)
    779 # The original time dimension is dropped from the dataset because
    780 # it becomes obsolete after the data variable is averaged. When the
    781 # averaged data variable is added to the dataset, the new time dimension
    782 # and its associated coordinates are also added.
    783 ds = ds.drop_dims(self.dim)  # type: ignore

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:1150, in TemporalAccessor._group_average(self, data_var, time_bounds)
   1145 # Ensure missing data (`np.nan`) receives no weight (zero). To
   1146 # achieve this, first broadcast the one-dimensional (temporal
   1147 # dimension) shape of the `weights` DataArray to the
   1148 # multi-dimensional shape of its corresponding data variable.
   1149 weights, _ = xr.broadcast(self._weights, dv)
-> 1150 weights = xr.where(np.isnan(dv), 0.0, weights)
   1152 # Perform weighted average using the formula
   1153 # WA = sum(data*weights) / sum(weights). The denominator must be
   1154 # included to take into account zero weight for missing data.
   1155 with xr.set_options(keep_attrs=True):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:1861, in where(cond, x, y, keep_attrs)
   1858     keep_attrs = _get_keep_attrs(default=False)
   1860 # alignment for three arguments is complicated, so don't support it yet
-> 1861 result = apply_ufunc(
   1862     duck_array_ops.where,
   1863     cond,
   1864     x,
   1865     y,
   1866     join="exact",
   1867     dataset_join="exact",
   1868     dask="allowed",
   1869     keep_attrs=keep_attrs,
   1870 )
   1872 # keep the attributes of x, the second parameter, by default to
   1873 # be consistent with the `where` method of `DataArray` and `Dataset`
   1874 # rebuild the attrs from x at each level of the output, which could be
   1875 # Dataset, DataArray, or Variable, and also handle coords
   1876 if keep_attrs is True and hasattr(result, "attrs"):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:1196, in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args)
   1194 # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
   1195 elif any(isinstance(a, DataArray) for a in args):
-> 1196     return apply_dataarray_vfunc(
   1197         variables_vfunc,
   1198         *args,
   1199         signature=signature,
   1200         join=join,
   1201         exclude_dims=exclude_dims,
   1202         keep_attrs=keep_attrs,
   1203     )
   1204 # feed Variables directly through apply_variable_ufunc
   1205 elif any(isinstance(a, Variable) for a in args):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:303, in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    298 result_coords, result_indexes = build_output_coords_and_indexes(
    299     args, signature, exclude_dims, combine_attrs=keep_attrs
    300 )
    302 data_vars = [getattr(a, "variable", a) for a in args]
--> 303 result_var = func(*data_vars)
    305 out: tuple[DataArray, ...] | DataArray
    306 if signature.num_outputs > 1:

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:760, in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    755     if vectorize:
    756         func = _vectorize(
    757             func, signature, output_dtypes=output_dtypes, exclude_dims=exclude_dims
    758         )
--> 760 result_data = func(*input_data)
    762 if signature.num_outputs == 1:
    763     result_data = (result_data,)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\duck_array_ops.py:300, in where(condition, x, y)
    298 """Three argument where() with better dtype promotion rules."""
    299 xp = get_array_namespace(condition)
--> 300 return xp.where(condition, *as_shared_dtype([x, y], xp=xp))

File <__array_function__ internals>:200, in where(*args, **kwargs)

MemoryError: Unable to allocate 41.1 GiB for an array with shape (3285, 1328, 1264) and data type float64

Whereas, import as "xr" (xarray) rather than "xc" (xcdat) allows the functionalities I'm talking about.
Climatology Fig from regular method which I always do:

But, the xcdat one fails.

With xCDAT memory issues are consistent (maybe I'm not handling the data well)

MemoryError                               Traceback (most recent call last)
Cell In[62], line 1
----> 1 mon_anm = xc_sic.temporal.departures("z", "month")
      2 mon_anm.z

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:669, in TemporalAccessor.departures(self, data_var, freq, weighted, keep_weights, reference_period, season_config)
    667 inferred_freq = self._infer_freq()
    668 if inferred_freq != freq:
--> 669     ds_obs = ds_obs.temporal.group_average(
    670         data_var,
    671         freq,
    672         weighted,
    673         keep_weights,
    674         season_config,
    675     )
    677 ds_climo = ds.temporal.climatology(
    678     data_var,
    679     freq,
   (...)
    683     season_config,
    684 )
    686 # 4. Group the averaged data variable values by the time `freq`.
    687 # ----------------------------------------------------------------------
    688 # This step allows us to perform xarray's grouped arithmetic to
    689 # calculate departures.

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:340, in TemporalAccessor.group_average(self, data_var, freq, weighted, keep_weights, season_config)
    215 """Returns a Dataset with average of a data variable by time group.
    216 
    217 Parameters
   (...)
    336 }
    337 """
    338 self._set_data_var_attrs(data_var)
--> 340 return self._averager(
    341     data_var,
    342     "group_average",
    343     freq,
    344     weighted=weighted,
    345     keep_weights=keep_weights,
    346     season_config=season_config,
    347 )

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:777, in TemporalAccessor._averager(self, data_var, mode, freq, weighted, keep_weights, reference_period, season_config)
    775     dv = self._average(dv, time_bounds)
    776 elif self._mode in ["group_average", "climatology", "departures"]:
--> 777     dv = self._group_average(dv, time_bounds)
    779 # The original time dimension is dropped from the dataset because
    780 # it becomes obsolete after the data variable is averaged. When the
    781 # averaged data variable is added to the dataset, the new time dimension
    782 # and its associated coordinates are also added.
    783 ds = ds.drop_dims(self.dim)  # type: ignore

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xcdat\temporal.py:1150, in TemporalAccessor._group_average(self, data_var, time_bounds)
   1145 # Ensure missing data (`np.nan`) receives no weight (zero). To
   1146 # achieve this, first broadcast the one-dimensional (temporal
   1147 # dimension) shape of the `weights` DataArray to the
   1148 # multi-dimensional shape of its corresponding data variable.
   1149 weights, _ = xr.broadcast(self._weights, dv)
-> 1150 weights = xr.where(np.isnan(dv), 0.0, weights)
   1152 # Perform weighted average using the formula
   1153 # WA = sum(data*weights) / sum(weights). The denominator must be
   1154 # included to take into account zero weight for missing data.
   1155 with xr.set_options(keep_attrs=True):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:1861, in where(cond, x, y, keep_attrs)
   1858     keep_attrs = _get_keep_attrs(default=False)
   1860 # alignment for three arguments is complicated, so don't support it yet
-> 1861 result = apply_ufunc(
   1862     duck_array_ops.where,
   1863     cond,
   1864     x,
   1865     y,
   1866     join="exact",
   1867     dataset_join="exact",
   1868     dask="allowed",
   1869     keep_attrs=keep_attrs,
   1870 )
   1872 # keep the attributes of x, the second parameter, by default to
   1873 # be consistent with the `where` method of `DataArray` and `Dataset`
   1874 # rebuild the attrs from x at each level of the output, which could be
   1875 # Dataset, DataArray, or Variable, and also handle coords
   1876 if keep_attrs is True and hasattr(result, "attrs"):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:1196, in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args)
   1194 # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
   1195 elif any(isinstance(a, DataArray) for a in args):
-> 1196     return apply_dataarray_vfunc(
   1197         variables_vfunc,
   1198         *args,
   1199         signature=signature,
   1200         join=join,
   1201         exclude_dims=exclude_dims,
   1202         keep_attrs=keep_attrs,
   1203     )
   1204 # feed Variables directly through apply_variable_ufunc
   1205 elif any(isinstance(a, Variable) for a in args):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:303, in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    298 result_coords, result_indexes = build_output_coords_and_indexes(
    299     args, signature, exclude_dims, combine_attrs=keep_attrs
    300 )
    302 data_vars = [getattr(a, "variable", a) for a in args]
--> 303 result_var = func(*data_vars)
    305 out: tuple[DataArray, ...] | DataArray
    306 if signature.num_outputs > 1:

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\computation.py:760, in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    755     if vectorize:
    756         func = _vectorize(
    757             func, signature, output_dtypes=output_dtypes, exclude_dims=exclude_dims
    758         )
--> 760 result_data = func(*input_data)
    762 if signature.num_outputs == 1:
    763     result_data = (result_data,)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\xarray\core\duck_array_ops.py:300, in where(condition, x, y)
    298 """Three argument where() with better dtype promotion rules."""
    299 xp = get_array_namespace(condition)
--> 300 return xp.where(condition, *as_shared_dtype([x, y], xp=xp))

File <__array_function__ internals>:200, in where(*args, **kwargs)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\array\core.py:1761, in Array.__array_function__(self, func, types, args, kwargs)
   1758 if has_keyword(da_func, "like"):
   1759     kwargs["like"] = self
-> 1761 return da_func(*args, **kwargs)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\array\routines.py:2107, in where(condition, x, y)
   2105     return broadcast_to(out, shape).astype(dtype)
   2106 else:
-> 2107     return elemwise(np.where, condition, x, y)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\array\core.py:4818, in elemwise(op, out, where, dtype, name, *args, **kwargs)
   4813     need_enforce_dtype = any(
   4814         not is_scalar_for_elemwise(a) and a.ndim == 0 for a in args
   4815     )
   4817 if not name:
-> 4818     name = f"{funcname(op)}-{tokenize(op, dtype, *args, where)}"
   4820 blockwise_kwargs = dict(dtype=dtype, name=name, token=funcname(op).strip("_"))
   4822 if where is not True:

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\base.py:930, in tokenize(*args, **kwargs)
    921 def tokenize(*args, **kwargs):
    922     """Deterministic token
    923 
    924     >>> tokenize([1, 2, '3'])
   (...)
    928     True
    929     """
--> 930     hasher = _md5(str(tuple(map(normalize_token, args))).encode())
    931     if kwargs:
    932         hasher.update(str(normalize_token(kwargs)).encode())

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\utils.py:642, in Dispatch.__call__(self, arg, *args, **kwargs)
    638 """
    639 Call the corresponding method based on type of argument.
    640 """
    641 meth = self.dispatch(type(arg))
--> 642 return meth(arg, *args, **kwargs)

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\dask\base.py:1219, in register_numpy.<locals>.normalize_array(x)
   1217 else:
   1218     try:
-> 1219         data = hash_buffer_hex(x.ravel(order="K").view("i1"))
   1220     except (BufferError, AttributeError, ValueError):
   1221         data = hash_buffer_hex(x.copy().ravel(order="K").view("i1"))

MemoryError: Unable to allocate 41.1 GiB for an array with shape (5517531904,) and data type float64

Just those point I mentioned, I'm bothered.
Aim: I need anomalies for the data mentioned uisng XCDAT's temporal departures.

[tomvothecoder]

Hi @gkb999, please follow this new template for asking questions in this forum. It really helps get things moving if questions and comments are succinct. Open-source projects (including ours) usually require users provide a minimum, reproducible code example with the relevant stack trace in proper markdown formatting.

If you want, you can transfer the content relevant from this post over to the new post.

I also encourage you to consider how you frame your questions and comments because it might come off as demanding to others in the open-source community (just an FYI!).

[lee1043]

@gkb999 I made some minor edit to your post to add formatting which helps improving readability. I second @tomvothecoder, following the template and providing a minimal code would be extremely helpful to trace down your issue.

[gkb999]

Well, is there a way I can post errors in those colors?

[lee1043]

@gkb999 you can click "edit" to your original post and see how I changed it. In front of your code part I added "python", and at the end I added "", which is going to look like below.

a = b + c

You can use "Preview" tab above your typing box, to check how it will look like once you post.

You can also google something like "Markdown syntax" to learn more. Here is one example resource: https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet

[gkb999]

Thanks a ton for enlightening me.
Now I know. ```python should be used !!!
Thanks.
Would follow from next posts.
Now should I transfer this question with these modifications?

[tomvothecoder]

@lee1043 Thank you for going out of your way to edit the posts. The new questions template will streamline posts and help a lot with readability.

@gkb999 Yes, please transfer the relevant content over to this template and fill out whatever is needed. The form fields for "Minimal Complete Verifiable Example (MVCE)" and `Relevant log output" will automatically format raw code/output into Python Markdown. There are descriptions under each form section. Thanks!

[tomvothecoder]

Closing this discussion to transfer over to a new discussion.