issue_comments

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.

That's right. I just did a quick 0.1 release of xncml, most likely rough around the edges. Give it a spin. PRs most welcome.

@rabernat If you're happy with it, this issue can probably be closed.

@andersy005 Sounds good !

Ok, another option would be to add that to xncml

@andersy005 What do you think ?

I've got a first draft that parses an NcML document and spits out an xarray.Dataset. It does not cover all the NcML syntax, but the essential elements are there.

It uses xsdata to parse the XML, using a datamodel automatically generated from the NcML 2-2 schema. I've scrapped test files from the netcdf-java repo to create a test suite.

Wondering what's the best place to host the code, tests and test data so others may give it a spin ?

The alignment of data and weights is not done automatically. So I agree this would be ideal, but I'll need some guidance to make it happen without the align call.

This is the error I get python ValueError: operand to apply_ufunc has required core dimensions ['dim_0', 'dim_1'], but some of these dimensions are absent on an input variable: ['dim_1']

@mathause I'm not 100% confident the methods other than linear work as expected, so I suggest we do not expose method until a robust test suite can confirm results are valid.

We have local changes that expose method in the public API, and compare results against numpy for the case with uniform weights to confirm results are identical.

Do you want us to push those changes here ?

Correct, it's not exposed yet because I don't have the bandwidth to create tests for all the different methods.

The equal weight case could be tested against numpy fairly easily though.

@mathause Should this PR be amended to account for #6108 ?

I had the same interrogation. My guess is that the DataArray.quantile and DataArrayWeighted.quantile should have similar signatures.

The way I see it is that another PR would bring the DataArray.quantile signature up to date with numpy 1.22, which includes the 9 h types, and harmonize DataArrayWeighted.quantile with it.

Is there an appetite to add a weighted quantile method as well ?

Thanks, I'll just pin pandas and wait for the next xarray release. Sorry about the noise.

I'd like to revive this issue. We're increasingly using NcML aggregations within our THREDDS server to create "logical" datasets. This allows us to fix some non-CF-conforming metadata fields without changing files on disk (which would break syncing with ESGF nodes). More importantly, by aggregating multiple time periods, variables and realizations, we're able to create catalog entries for simulations instead of files, which we expect will greatly facilitate parsing catalog search results. We'd like to offer the same aggregation functionality outside of the THREDDS server. Ideally, this would be supported right from the netcdf-c library (see https://github.com/Unidata/netcdf-c/issues/1478), but an xarray NcML backend is the second best option. I also imagine that NcML files could be use as a clean mechanism to create Zarr/NCZarr objects ie: *.nc -> open_ncml -> xr.Dataset -> to_zarr -> Zarr store

@andersy005 In terms of API, I think the need is not so much to create or modify NcML files, but rather to return an xarray.Dataset from an NcML description. My understanding is that open_ncml would be a wrapper around open_mfdataset. My hope is that NcML-based xarray.Dataset objects would behave similarly whether they are created from files on disk through xarray.open_ncml('sim.ncml') or xarray.open_dataset('https://.../thredds/sim.ncml').

The THREDDS repo contains a number of unit tests that could be emulated to steer the Python implementation. My understanding is that getting this done could involve a fair amount of work, so I'd like to see who's interested in collaborating on this and maybe schedule a meeting to plan work for this year or the next.

@max-sixty Thanks for taking the time. This PR and #3733 are needed for quantile mapping methods @aulemahal and I are working on.

@shoyer I'm having trouble wrapping my head around this. The example above is essentially a 1D interpolation over y, a coordinate selection over x, and a broadcast over z. I don't think it should interpolate over x, right ?

Just got bit by this as well. Computing monthly quantile correction factors, so I have an array with dimensions (month, quantile, lon, lat). I then want to apply these correction factors to a time series (time, lon, lat), so I compute the month and quantile of my time series, and want to interp into the quantile correction factors. This doesn't work because both the factors and the time series have (lat, lon) dimensions.

Thanks @spencerkclark for shepherding this to completion.

Ah ! No, I tried numpy.testing's version. Works now. Thanks !

It seems like assert_allclose doesn't like nans. Other suggestion ?

I think so. If there's a branch fixing the assert_allclose failures, I can merge it here.

Merged with #3631.

Note that at the moment, if we pass np.datetime64 objects that exceed the allowed time span, the function yields garbage without failing. Is this something we want to fix as well ?

One option is to convert array and offset to microseconds first, then compute the delta, but this may break people's code.

How about replacing array = np.asarray(pd.Series(array.ravel())).reshape(array.shape) by array = array.astype("timedelta64") ? with numpy 1.17 your example works and the test suite only fails on unrelated netcdf string errors.

Got it, thanks !

Another issue with datetime_to_numeric happens with: import xarray as xr import cftime i = xr.CFTimeIndex(xr.cftime_range('2000-01-01', periods=2)) xr.core.duck_array_ops.datetime_to_numeric(i, cftime.DatetimeGregorian(2, 1, 1), datetime_unit='D')

```python

TypeError Traceback (most recent call last) pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.array_to_timedelta64()

pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.parse_timedelta_string()

TypeError: object of type 'datetime.timedelta' has no len()

During handling of the above exception, another exception occurred:

OverflowError Traceback (most recent call last) <ipython-input-50-b03d9c4f220d> in <module> ----> 1 xr.core.duck_array_ops.datetime_to_numeric(i, cftime.DatetimeGregorian(2, 1, 1), datetime_unit='D')

~/src/xarray/xarray/core/duck_array_ops.py in datetime_to_numeric(array, offset, datetime_unit, dtype) 395 else: 396 offset = min(array) --> 397 array = array - offset 398 399 if not hasattr(array, "dtype"): # scalar is converted to 0d-array

~/src/xarray/xarray/coding/cftimeindex.py in sub(self, other) 431 432 if isinstance(other, (CFTimeIndex, cftime.datetime)): --> 433 return pd.TimedeltaIndex(np.array(self) - np.array(other)) 434 elif isinstance(other, pd.TimedeltaIndex): 435 return CFTimeIndex(np.array(self) - other.to_pytimedelta())

~/.conda/envs/xclim3/lib/python3.6/site-packages/pandas/core/indexes/timedeltas.py in new(cls, data, unit, freq, start, end, periods, closed, dtype, copy, name, verify_integrity) 256 257 tdarr = TimedeltaArray._from_sequence( --> 258 data, freq=freq, unit=unit, dtype=dtype, copy=copy 259 ) 260 return cls._simple_new(tdarr._data, freq=tdarr.freq, name=name)

~/.conda/envs/xclim3/lib/python3.6/site-packages/pandas/core/arrays/timedeltas.py in _from_sequence(cls, data, dtype, copy, freq, unit) 270 freq, freq_infer = dtl.maybe_infer_freq(freq) 271 --> 272 data, inferred_freq = sequence_to_td64ns(data, copy=copy, unit=unit) 273 freq, freq_infer = dtl.validate_inferred_freq(freq, inferred_freq, freq_infer) 274

~/.conda/envs/xclim3/lib/python3.6/site-packages/pandas/core/arrays/timedeltas.py in sequence_to_td64ns(data, copy, unit, errors) 971 if is_object_dtype(data.dtype) or is_string_dtype(data.dtype): 972 # no need to make a copy, need to convert if string-dtyped --> 973 data = objects_to_td64ns(data, unit=unit, errors=errors) 974 copy = False 975

~/.conda/envs/xclim3/lib/python3.6/site-packages/pandas/core/arrays/timedeltas.py in objects_to_td64ns(data, unit, errors) 1096 values = np.array(data, dtype=np.object_, copy=False) 1097 -> 1098 result = array_to_timedelta64(values, unit=unit, errors=errors) 1099 return result.view("timedelta64[ns]") 1100

pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.array_to_timedelta64()

pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.convert_to_timedelta64()

pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.delta_to_nanoseconds()

OverflowError: Python int too large to convert to C long ```

@clyne Let me rephrase my question: how do you feel about xarray providing a polyfit/polyval implementation essentially duplicating GeoCat's implementation ?

@maboualidev Is your objective to integrate the GeoCat implementation into xarray or keep it standalone ?

On my end, I'll submit a PR to add support for non-standard calendars to xarray.core.missing.get_clean_interp, which you'd then be able to use to get x values from coordinates.

@maboualidev Nice ! I see you're storing the residuals in the DataArray attributes. From my perspective, it would be useful to have those directly as DataArrays. Thoughts ?

So it looks like there are multiple inspirations to draw from. Here is what I could gather.

• xscale.signal.fitting.polyfit(obj, deg=1, dim=None, coord=None) supports chunking along the fitting dimension using dask.array.linalg.lstsq. No explicit missing data handling.
• xyzpy.signal.xr_polyfit(obj, dim, ix=None, deg=0.5, poly='hermite') applies np.polynomial.polynomial.polyfit using xr.apply_ufunc along dim with the help of numba. Also supports other types of polynomial (legendre, chebyshev, ...). Missing values are masked out 1D wise.
• geocat.comp.ndpolyfit(x: Iterable, y: Iterable, deg: int, axis: int = 0, **kwargs) -> (xr.DataArray, da.Array) reorders the array to apply np.polyfit along dim, returns the full outputs (residuals, rank, etc) as DataArray attributes. Missing values are masked out in bulk if possible, 1D-wise otherwise.

There does not seem to be matching polyval implementations for any of those nor support for indexing along a time dimension with a non-standard calendar.

My current implementation is pretty naive. It's just calling numpy.polyfit using dask.array.apply_along_axis. Happy to put that in a PR as a starting point, but there are a couple of questions I had: * How to return the full output (residuals, rank, singular_values, rcond) ? A tuple of dataarrays or a dataset ? * Do we want to use the dask least square functionality to allow for chunking within the x dimension ? Then it's not just a simple wrapper around polyfit. * Should we use np.polyfit or np.polynomial.polynomial.polyfit ?

Thanks, it seems to work !

Started to work on this and facing some issues with the x-coordinate when its a datetime. For standard calendars, I can use pd.to_numeric(da.time), but for non-standard calendars, it's not clear how to go ahead. If I use xr.coding.times.encode_cf_datetime(coord), the coefficients we'll find will only make sense in the polyval function if we use the same time encoding.

Ok, I'll submit a PR shortly.

Looking at the code, I'm confused. The DataArray.quantile method creates a temporary dataset, copies the variable over, calls the Variable.quantile method, then assigns the attributes from the dataset to this new variable. At no point however are attributes assigned to this temporary dataset. My understanding is that Variable.quantile should have a keep_attrs argument, correct ?

The PR only covered DataArrays, not Datasets, so yes, this is expected.

Ok thanks, will submit to pydap.

Made the change. I had to drop the quantile dimension for scalar q after the apply to maintain the same behavior as DataArray.quantile, but it got me wondering if keeping the quantile dimension even for scalar q should not be the default.

I just realized that I could do da.groupby('time.month').apply(xr.DataArray.quantile, q=.5, dim='time') which seems to yield the same results.

Got it, thanks for the workaround !

Do you think there would be a benefit to implementing a TimeGrouper class based on panda's ?

I'm trying to wrap my head around what is needed to get the resample method to work but I must say I'm confused. Would it be possible/practical to create a branch with stubs in the code for the methods that need to be written (with a #2191 comment) so newbies can help fill-in the gaps?