From a17d0193c19e9935ec54664a51da7501a7a9ea8c Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Thu, 4 Dec 2025 16:26:04 +0100
Subject: [PATCH] add an xarray wrapper

---
 pyproject.toml               |  3 ++-
 rqadeforestation/__init__.py | 22 +++++++++++++++++++---
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 7cf7f12..603290d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,7 +15,8 @@ readme = "README.md"
 requires-python = ">=3.8"
 dependencies = [
     "numpy",
-    "dask"
+    "dask", 
+    "xarray"
 ]
 classifiers = [
     "Programming Language :: Python :: 3",
diff --git a/rqadeforestation/__init__.py b/rqadeforestation/__init__.py
index fd202eb..5d08531 100644
--- a/rqadeforestation/__init__.py
+++ b/rqadeforestation/__init__.py
@@ -2,6 +2,7 @@
 import ctypes as ct
 import numpy as np
 import dask.array as da
+import xarray as xr
 
 class MallocVector(ct.Structure):
     _fields_ = [("pointer", ct.c_void_p),
@@ -38,7 +39,7 @@ def rqatrend(y: np.ndarray, threshold: float, border: int = 10, theiler: int = 1
     :param theiler: Theiler window size for the RQA calculation.
     :return: The RQA trend value.
     """
-    py = mvptr(y)
+    py = mvptr(y.astype(np.float64))
     lib.rqatrend.argtypes = (ct.POINTER(MallocVector), ct.c_double, ct.c_int64, ct.c_int64)
     lib.rqatrend.restype = ct.c_double
     result_single = lib.rqatrend(py, threshold, border, theiler)
@@ -56,16 +57,31 @@ def rqatrend_matrix(matrix: np.ndarray, threshold: float, border: int = 10, thei
     :return: Numpy array of all RQA trend values of size n_timeseries.
     """
 
+    if not len(matrix.shape) == 2:
+        raise Exception("Input to rqatrend_matrix must be 2d") 
+
     n = matrix.shape[0]
     result_several = np.ones(n)
     p_result_several = mvptr(result_several)
-    p_matrix = mmptr(matrix)
-    
+    p_matrix = mmptr(matrix.astype(np.float64))
     # arguments: result_vector, data, threshhold, border, theiler
     lib.rqatrend_inplace.argtypes = (ct.POINTER(MallocVector), ct.POINTER(MallocMatrix), ct.c_double, ct.c_int64, ct.c_int64)
     return_value = lib.rqatrend_inplace(p_result_several, p_matrix, threshold, border, theiler)
     return result_several
 
+def rqatrend_xarray(x, threshold:float, border: int = 10, 
+                    theiler: int=1, out_dtype = np.float64,
+                    timeaxis_name = "time"):
+    return xr.apply_ufunc(
+        rqatrend,
+        x.chunk({timeaxis_name: -1}),
+        kwargs = {'threshold': threshold,'border':border,'theiler':theiler},
+        input_core_dims = [[timeaxis_name]],
+        output_core_dims = [[]],
+        dask = "parallelized",
+        vectorize=True,
+    )
+
 
 def rqatrend_dask(x: da.Array, timeseries_axis: int, threshold: float, border: int = 10, theiler: int = 1, out_dtype: type = np.float64) -> da.Array:
     """