Slope of linear fit: Using xskillscore.linslope

[HARSHBARDHANK94]

import glob
import numpy as np
import xarray as xr
import xskillscore as xs

# fetching files from foldersb and opening using mfdataset()............................................
filename1= glob.glob("E:\\ACI_data\\AOD\\*.nc")
data1=xr.open_mfdataset(filename1)

# Region of Interest Subsetting...........................................................
AOI1=data1.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))

# Variable selection from the dataset ................................................
AOD=AOI1.MYD08_M3_6_1_AOD_550_Dark_Target_Deep_Blue_Combined_Mean_Mean

The Dimension of the Variables are:
lat=28
lon=62
time=216
AOD=(216,28,62)

In xskillscore.linslope (), parameters requires two labeled array(s) (e.g., a and b) as detailed in documentation: https://xskillscore.readthedocs.io/en/stable/api/xskillscore.linslope.html

But I want to compute the linear slope against time in AOD. As my AOD variable has 216 month of time dimension, I want to compute linear fit between time and AOD, that should produce slope of trend. May be I getting the documentation wrong or something different.........................

slope=xs.linslope (AOD.chunk(dict(time=-1)),...............b................., dim='time')

What I should declare in the above line at position 'b'

Thanks for constructive comments!!

[aaronspring]

Try b=AOD.time or np.arange(AOD.time.size)

[HARSHBARDHANK94]

np.arange(AOD.time.size) didn't work, bcz numpy array has not 'time' dim and in parameter call both data arrays should be labeled and having dim='time'
So again I supplied:

slope=xs.linslope(AOD.chunk(dict(time=-1)), AOD.time, dim='time')
slope.values

produced error..........

UFuncTypeError                            Traceback (most recent call last)
C:\Users\HARSHB~1\AppData\Local\Temp/ipykernel_11224/158551320.py in <module>
      1 slope=xs.linslope(AOD.chunk(dict(time=-1)),AOI1.time, dim='time')
----> 2 slope.values

~\anaconda3\lib\site-packages\xarray\core\dataarray.py in values(self)
    649         type does not support coercion like this (e.g. cupy).
    650         """
--> 651         return self.variable.values
    652 
    653     @values.setter

~\anaconda3\lib\site-packages\xarray\core\variable.py in values(self)
    515     def values(self):
    516         """The variable's data as a numpy.ndarray"""
--> 517         return _as_array_or_item(self._data)
    518 
    519     @values.setter

~\anaconda3\lib\site-packages\xarray\core\variable.py in _as_array_or_item(data)
    257     TODO: remove this (replace with np.asarray) once these issues are fixed
    258     """
--> 259     data = np.asarray(data)
    260     if data.ndim == 0:
    261         if data.dtype.kind == "M":

~\anaconda3\lib\site-packages\numpy\core\_asarray.py in asarray(a, dtype, order, like)
    100         return _asarray_with_like(a, dtype=dtype, order=order, like=like)
    101 
--> 102     return array(a, dtype, copy=False, order=order)
    103 
    104 

~\anaconda3\lib\site-packages\dask\array\core.py in __array__(self, dtype, **kwargs)
   1532 
   1533     def __array__(self, dtype=None, **kwargs):
-> 1534         x = self.compute()
   1535         if dtype and x.dtype != dtype:
   1536             x = x.astype(dtype)

~\anaconda3\lib\site-packages\dask\base.py in compute(self, **kwargs)
    286         dask.base.compute
    287         """
--> 288         (result,) = compute(self, traverse=False, **kwargs)
    289         return result
    290 

~\anaconda3\lib\site-packages\dask\base.py in compute(*args, **kwargs)
    568         postcomputes.append(x.__dask_postcompute__())
    569 
--> 570     results = schedule(dsk, keys, **kwargs)
    571     return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])
    572 

~\anaconda3\lib\site-packages\dask\threaded.py in get(dsk, result, cache, num_workers, pool, **kwargs)
     77             pool = MultiprocessingPoolExecutor(pool)
     78 
---> 79     results = get_async(
     80         pool.submit,
     81         pool._max_workers,

~\anaconda3\lib\site-packages\dask\local.py in get_async(submit, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, chunksize, **kwargs)
    505                             _execute_task(task, data)  # Re-execute locally
    506                         else:
--> 507                             raise_exception(exc, tb)
    508                     res, worker_id = loads(res_info)
    509                     state["cache"][key] = res

~\anaconda3\lib\site-packages\dask\local.py in reraise(exc, tb)
    313     if exc.__traceback__ is not tb:
    314         raise exc.with_traceback(tb)
--> 315     raise exc
    316 
    317 

~\anaconda3\lib\site-packages\dask\local.py in execute_task(key, task_info, dumps, loads, get_id, pack_exception)
    218     try:
    219         task, data = loads(task_info)
--> 220         result = _execute_task(task, data)
    221         id = get_id()
    222         result = dumps((result, id))

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\core.py in <genexpr>(.0)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\optimization.py in __call__(self, *args)
    967         if not len(args) == len(self.inkeys):
    968             raise ValueError("Expected %d args, got %d" % (len(self.inkeys), len(args)))
--> 969         return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))
    970 
    971     def __reduce__(self):

~\anaconda3\lib\site-packages\dask\core.py in get(dsk, out, cache)
    147     for key in toposort(dsk):
    148         task = dsk[key]
--> 149         result = _execute_task(task, cache)
    150         cache[key] = result
    151     result = _execute_task(out, cache)

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\xskillscore\core\np_deterministic.py in _linslope(a, b, weights, axis, skipna)
    189         weights = np.rollaxis(weights, axis)
    190 
--> 191     am, bm = __compute_anomalies(a, b, weights=weights, axis=0, skipna=skipna)
    192 
    193     if weights is not None:

~\anaconda3\lib\site-packages\xskillscore\core\np_deterministic.py in __compute_anomalies(a, b, weights, axis, skipna)
     97         with suppress_warnings("Mean of empty slice"):
     98             ma = meanfunc(a, axis=axis)
---> 99             mb = meanfunc(b, axis=axis)
    100     am, bm = a - ma, b - mb
    101     return am, bm

<__array_function__ internals> in mean(*args, **kwargs)

~\anaconda3\lib\site-packages\numpy\core\fromnumeric.py in mean(a, axis, dtype, out, keepdims, where)
   3417             return mean(axis=axis, dtype=dtype, out=out, **kwargs)
   3418 
-> 3419     return _methods._mean(a, axis=axis, dtype=dtype,
   3420                           out=out, **kwargs)
   3421 

~\anaconda3\lib\site-packages\numpy\core\_methods.py in _mean(a, axis, dtype, out, keepdims, where)
    176             is_float16_result = True
    177 
--> 178     ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)
    179     if isinstance(ret, mu.ndarray):
    180         ret = um.true_divide(

UFuncTypeError: ufunc 'add' cannot use operands with types dtype('<M8[ns]') and dtype('<M8[ns]')

[HARSHBARDHANK94]

[HARSHBARDHANK94]

When I added skipna=True

slope=xs.linslope(AOD.chunk(dict(time=-1)),AOD.time, dim='time', skipna=True)
slope.values

ValueError                                Traceback (most recent call last)
C:\Users\HARSHB~1\AppData\Local\Temp/ipykernel_11224/3043353494.py in <module>
      1 slope=xs.linslope(AOD.chunk(dict(time=-1)),AOD.time, dim='time', skipna=True)
----> 2 slope.values

~\anaconda3\lib\site-packages\xarray\core\dataarray.py in values(self)
    649         type does not support coercion like this (e.g. cupy).
    650         """
--> 651         return self.variable.values
    652 
    653     @values.setter

~\anaconda3\lib\site-packages\xarray\core\variable.py in values(self)
    515     def values(self):
    516         """The variable's data as a numpy.ndarray"""
--> 517         return _as_array_or_item(self._data)
    518 
    519     @values.setter

~\anaconda3\lib\site-packages\xarray\core\variable.py in _as_array_or_item(data)
    257     TODO: remove this (replace with np.asarray) once these issues are fixed
    258     """
--> 259     data = np.asarray(data)
    260     if data.ndim == 0:
    261         if data.dtype.kind == "M":

~\anaconda3\lib\site-packages\numpy\core\_asarray.py in asarray(a, dtype, order, like)
    100         return _asarray_with_like(a, dtype=dtype, order=order, like=like)
    101 
--> 102     return array(a, dtype, copy=False, order=order)
    103 
    104 

~\anaconda3\lib\site-packages\dask\array\core.py in __array__(self, dtype, **kwargs)
   1532 
   1533     def __array__(self, dtype=None, **kwargs):
-> 1534         x = self.compute()
   1535         if dtype and x.dtype != dtype:
   1536             x = x.astype(dtype)

~\anaconda3\lib\site-packages\dask\base.py in compute(self, **kwargs)
    286         dask.base.compute
    287         """
--> 288         (result,) = compute(self, traverse=False, **kwargs)
    289         return result
    290 

~\anaconda3\lib\site-packages\dask\base.py in compute(*args, **kwargs)
    568         postcomputes.append(x.__dask_postcompute__())
    569 
--> 570     results = schedule(dsk, keys, **kwargs)
    571     return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])
    572 

~\anaconda3\lib\site-packages\dask\threaded.py in get(dsk, result, cache, num_workers, pool, **kwargs)
     77             pool = MultiprocessingPoolExecutor(pool)
     78 
---> 79     results = get_async(
     80         pool.submit,
     81         pool._max_workers,

~\anaconda3\lib\site-packages\dask\local.py in get_async(submit, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, chunksize, **kwargs)
    505                             _execute_task(task, data)  # Re-execute locally
    506                         else:
--> 507                             raise_exception(exc, tb)
    508                     res, worker_id = loads(res_info)
    509                     state["cache"][key] = res

~\anaconda3\lib\site-packages\dask\local.py in reraise(exc, tb)
    313     if exc.__traceback__ is not tb:
    314         raise exc.with_traceback(tb)
--> 315     raise exc
    316 
    317 

~\anaconda3\lib\site-packages\dask\local.py in execute_task(key, task_info, dumps, loads, get_id, pack_exception)
    218     try:
    219         task, data = loads(task_info)
--> 220         result = _execute_task(task, data)
    221         id = get_id()
    222         result = dumps((result, id))

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\core.py in <genexpr>(.0)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\dask\optimization.py in __call__(self, *args)
    967         if not len(args) == len(self.inkeys):
    968             raise ValueError("Expected %d args, got %d" % (len(self.inkeys), len(args)))
--> 969         return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))
    970 
    971     def __reduce__(self):

~\anaconda3\lib\site-packages\dask\core.py in get(dsk, out, cache)
    147     for key in toposort(dsk):
    148         task = dsk[key]
--> 149         result = _execute_task(task, cache)
    150         cache[key] = result
    151     result = _execute_task(out, cache)

~\anaconda3\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    117         # temporaries by their reference count and can execute certain
    118         # operations in-place.
--> 119         return func(*(_execute_task(a, cache) for a in args))
    120     elif not ishashable(arg):
    121         return arg

~\anaconda3\lib\site-packages\xskillscore\core\np_deterministic.py in _linslope(a, b, weights, axis, skipna)
    182     sumfunc, meanfunc = _get_numpy_funcs(skipna)
    183     if skipna:
--> 184         a, b, weights = _match_nans(a, b, weights)
    185     weights = _check_weights(weights)
    186     a = np.rollaxis(a, axis)

~\anaconda3\lib\site-packages\xskillscore\core\np_deterministic.py in _match_nans(a, b, weights)
     44         # Find pairwise indices in a and b that have nans.
     45         idx = np.logical_or(np.isnan(a), np.isnan(b))
---> 46         a[idx], b[idx] = np.nan, np.nan
     47         # https://github.com/xarray-contrib/xskillscore/issues/168
     48         if isinstance(weights, np.ndarray):

ValueError: Could not convert object to NumPy datetime

[aaronspring]

please try without using dask, so the traceback can be better debugged.

slope=xs.linslope(AOD.chunk(dict(time=-1)), xr.DataArray(np.arange(AOD.time.size),dims='time'), dim='time')

[aaronspring]

the linslope function doesnt work with cftime often. Therefore try setting time to an equally spaced range:

import xskillscore as xs
ds = xr.tutorial.load_dataset('rasm')
ds['time']=np.arange(ds.time.size)
xs.linslope(ds, ds.time, dim='time')

[HARSHBARDHANK94]

the linslope function doesnt work with cftime often. Therefore try setting time to an equally spaced range:

import xskillscore as xs
ds = xr.tutorial.load_dataset('rasm')
ds['time']=np.arange(ds.time.size)
xs.linslope(ds, ds.time, dim='time')

Hello, this worked well........Thanks,

I want to know does the xskillscore.linslope works same as the scipy.stats.linregress to calculate the slope.
Bcz I tried to loop through every lat-lon along the time dimension and tried to calculate the slope using scipy.stats.linregress, I got different result than this xs.linslope(ds, ds.time, dim='time'). Also xs.linslope does not return p_value.

I used the following lines for that:

from scipy.stats import linregress
array1=np.array([])
list1=[]
p_array1= np.array([])
p_list1=[]

for i in range(0,56375):                                                       # 205*275=56375
    slope, intercept, r, p,se = linregress (time [:,0], value [:,i])
    list1.append(slope), p_list1.append (p)

array1 = np.append (array1, list1)
slope = array1.reshape (205,275)                                      # return array (56375) to lat lon (205*275)

p_array1 = np.append (p_array1,p_list1)
p_value = p_array1.reshape(205,275)                                # return array (56375) to lat lon (205*275)