diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 898941a3f9ddc89dbec6ff8f0a3092394af08901..0aab1c4918ac300c8371f145d715670a3a2950cd 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -11,7 +11,7 @@ test_PyWake:  # name the job what we like
   stage:  # build, test, deploy defined by default [2]
     test
   script:
-  - pip install -e .[test]
+  - pip install -e .[test] --timeout 60
   - pytest
   tags:  # only runners with this tag can do the job [3]
   - ci-ubuntu
diff --git a/py_wake/site/_site.py b/py_wake/site/_site.py
index a51ea3e76702002dd3c5d62f5b857f79cc459751..db71b46dd3b9c34bca0306d44b7f4cd6c41ff4f9 100644
--- a/py_wake/site/_site.py
+++ b/py_wake/site/_site.py
@@ -4,6 +4,7 @@ from py_wake.site.shear import PowerShear
 import py_wake.utils.xarray_utils  # register ilk function @UnusedImport
 import xarray as xr
 from abc import ABC, abstractmethod
+from py_wake.utils.xarray_utils import da2py, coords2py
 
 """
 suffixs:
@@ -48,9 +49,9 @@ class LocalWind(xr.Dataset):
         for k, v in [('x', x_i), ('y', y_i), ('h', h_i)]:
             if v is not None:
                 coords[k] = ('i', np.zeros(n_i) + v)
-        xr.Dataset.__init__(self, data_vars={k: v for k, v in [('WD', WD), ('WS', WS),
-                                                               ('TI', TI), ('P', P)] if v is not None},
-                            coords=coords)
+        xr.Dataset.__init__(self, data_vars={k: da2py(v) for k, v in [('WD', WD), ('WS', WS),
+                                                                      ('TI', TI), ('P', P)] if v is not None},
+                            coords={k: coords2py(v) for k, v in coords.items()})
         self.attrs['wd_bin_size'] = wd_bin_size
 
         # set localWind.WS_ilk etc.
diff --git a/py_wake/utils/xarray_utils.py b/py_wake/utils/xarray_utils.py
index 03aa671a1cd2b2d8d7c0e0ba63fedb0f89b74d35..33a5ced075c82b133419f3fe3f831d3222f868b1 100644
--- a/py_wake/utils/xarray_utils.py
+++ b/py_wake/utils/xarray_utils.py
@@ -4,6 +4,7 @@ import numpy as np
 import xarray as xr
 from xarray.plot.plot import _PlotMethods
 import warnings
+from xarray.core.dataarray import DataArray
 
 
 class ilk():
@@ -52,7 +53,7 @@ class add_ilk():
                     break
         while len(np.shape(value)) > len(d) and np.shape(value)[-1] == 1:
             value = value[..., 0]
-        self.dataset[name] = (d, value)
+        self.dataset[name] = (d, da2py(value, include_dims=False))
 
 
 class add_ijlk():
@@ -127,3 +128,18 @@ if not hasattr(xr.DataArray(None), 'ilk'):
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
         xr.register_dataarray_accessor('plot')(plot_xy_map)
+
+
+def da2py(v, include_dims=True):
+    if isinstance(v, DataArray):
+        if include_dims:
+            return (v.dims, v.values)
+        else:
+            return v.values
+    return v
+
+
+def coords2py(v):
+    if isinstance(v, tuple):
+        return tuple([da2py(v, False) for v in v])
+    return v
diff --git a/py_wake/wind_farm_models/wind_farm_model.py b/py_wake/wind_farm_models/wind_farm_model.py
index 26963921b2436437fc3fc4d4a3ac1b06d2d2c2f9..2ea5016304410dada2e805156113838abe6a0afc 100644
--- a/py_wake/wind_farm_models/wind_farm_model.py
+++ b/py_wake/wind_farm_models/wind_farm_model.py
@@ -7,6 +7,7 @@ import xarray as xr
 from py_wake.utils import xarray_utils, weibull  # register ilk function @UnusedImport
 from numpy import newaxis as na
 from py_wake.utils.model_utils import check_model, fix_shape
+from py_wake.utils.xarray_utils import da2py
 
 
 class WindFarmModel(ABC):
@@ -59,7 +60,8 @@ class WindFarmModel(ABC):
 
         if len(x) == 0:
             lw = UniformSite([1], 0.1).local_wind(x_i=[], y_i=[], h_i=[], wd=wd, ws=ws)
-            z = xr.DataArray(np.zeros((0, len(lw.wd), len(lw.ws))), coords=[('wt', []), ('wd', lw.wd), ('ws', lw.ws)])
+            z = xr.DataArray(np.zeros((0, len(lw.wd), len(lw.ws))), coords=[('wt', []), ('wd', da2py(lw.wd, False)),
+                                                                            ('ws', da2py(lw.ws, False))])
             return SimulationResult(self, lw, [], yaw, tilt, z, z, z, z, kwargs)
         res = self.calc_wt_interaction(x_i=np.asarray(x), y_i=np.asarray(y), h_i=h, type_i=type,
                                        yaw_ilk=yaw_ilk, tilt_ilk=tilt_ilk,
@@ -168,13 +170,14 @@ class SimulationResult(xr.Dataset):
 
         ilk_dims = (['wt', 'wd', 'ws'], ['wt', 'time'])['time' in lw]
         xr.Dataset.__init__(self,
-                            data_vars={k: (ilk_dims, (v, v[:, :, 0])['time' in lw], {'Description': d}) for k, v, d in [
-                                ('WS_eff', WS_eff_ilk, 'Effective local wind speed [m/s]'),
-                                ('TI_eff', np.zeros_like(WS_eff_ilk) + TI_eff_ilk,
-                                 'Effective local turbulence intensity'),
-                                ('Power', power_ilk, 'Power [W]'),
-                                ('CT', ct_ilk, 'Thrust coefficient'),
-                            ]},
+                            data_vars={k: (ilk_dims, da2py((v, v[:, :, 0])['time' in lw], include_dims=False),
+                                           {'Description': d})
+                                       for k, v, d in [('WS_eff', WS_eff_ilk, 'Effective local wind speed [m/s]'),
+                                                       ('TI_eff', np.zeros_like(WS_eff_ilk) + TI_eff_ilk,
+                                                        'Effective local turbulence intensity'),
+                                                       ('Power', power_ilk, 'Power [W]'),
+                                                       ('CT', ct_ilk, 'Thrust coefficient'),
+                                                       ]},
                             coords=coords)
         for n in localWind:
             self[n] = localWind[n]