diff --git a/py_wake/site/_site.py b/py_wake/site/_site.py
index db71b46dd3b9c34bca0306d44b7f4cd6c41ff4f9..97e6ac635878916163fb6753ad52785e84457376 100644
--- a/py_wake/site/_site.py
+++ b/py_wake/site/_site.py
@@ -4,7 +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
+from py_wake.utils.xarray_utils import da2py
 
 """
 suffixs:
@@ -49,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: 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()})
+        xr.Dataset.__init__(self, data_vars={k: da2py(v, include_dims=True) for k, v in [('WD', WD), ('WS', WS),
+                                                                                         ('TI', TI), ('P', P)] if v is not None},
+                            coords={k: da2py(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 33a5ced075c82b133419f3fe3f831d3222f868b1..e11a422e1231ea7da9f4aefc181f4970d33aea0c 100644
--- a/py_wake/utils/xarray_utils.py
+++ b/py_wake/utils/xarray_utils.py
@@ -130,16 +130,12 @@ if not hasattr(xr.DataArray(None), 'ilk'):
         xr.register_dataarray_accessor('plot')(plot_xy_map)
 
 
-def da2py(v, include_dims=True):
+def da2py(v, include_dims=False):
+    if isinstance(v, tuple):
+        return tuple([da2py(v, include_dims) for v in v])
     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 2ea5016304410dada2e805156113838abe6a0afc..dcd0e62d487d2fc013803c585f718b61fd789470 100644
--- a/py_wake/wind_farm_models/wind_farm_model.py
+++ b/py_wake/wind_farm_models/wind_farm_model.py
@@ -60,8 +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', da2py(lw.wd, False)),
-                                                                            ('ws', da2py(lw.ws, False))])
+            z = xr.DataArray(np.zeros((0, len(lw.wd), len(lw.ws))), coords=[('wt', []), ('wd', da2py(lw.wd)),
+                                                                            ('ws', da2py(lw.ws))])
             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,
@@ -170,7 +170,7 @@ class SimulationResult(xr.Dataset):
 
         ilk_dims = (['wt', 'wd', 'ws'], ['wt', 'time'])['time' in lw]
         xr.Dataset.__init__(self,
-                            data_vars={k: (ilk_dims, da2py((v, v[:, :, 0])['time' in lw], include_dims=False),
+                            data_vars={k: (ilk_dims, da2py((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,