diff --git a/py_wake/site/xrsite.py b/py_wake/site/xrsite.py
index db92dca227f75fed6874bd5d492c088f4a52bd63..e892fe8f9d23ae47f3ec8eee5a106fa85952b1e0 100644
--- a/py_wake/site/xrsite.py
+++ b/py_wake/site/xrsite.py
@@ -128,14 +128,11 @@ class XRSite(Site):
         else:
             data, l_indices = pre_sel(data, 'wd')
 
-        if 'i' in ip_dims:
-            if 'i' in coords and len(var.i) != len(coords['i']):
-                raise ValueError(
-                    "Number of points, i(=%d), in site data variable, %s, must match number of requested points(=%d)" %
-                    (len(var.i), var.name, len(coords['i'])))
-            # requesting all points(wt positions) in site
-            # ip_dims.remove('i')
-            # ip_data_dims = ['i']
+        if 'i' in ip_dims and 'i' in coords and len(var.i) != len(coords['i']):
+            raise ValueError(
+                "Number of points, i(=%d), in site data variable, %s, must match number of requested points(=%d)" %
+                (len(var.i), var.name, len(coords['i'])))
+        data, i_indices = pre_sel(data, 'i')
 
         if len(ip_dims) > 0:
             grid_interp = GridInterpolator([var.coords[k].data for k in ip_dims], data,
@@ -167,8 +164,9 @@ class XRSite(Site):
             ip_data = data
             ip_data_dims = []
 
-#         if 'i' in var.dims:
-#             ip_data_dims.insert(0, 'i')
+        if i_indices is not None:
+            ip_data_dims.append('i')
+            ip_data = sel(ip_data, ip_data_dims, i_indices, 'i')
         if l_indices is not None:
             if 'time' in coords:
                 ip_data_dims.append('time')
diff --git a/py_wake/tests/test_sites/test_xrsite.py b/py_wake/tests/test_sites/test_xrsite.py
index 3f7c9c539feb5b2f56ee990390278de8f4557ba5..ab3a71e8f34e429d0414fd62aec705832880d894 100644
--- a/py_wake/tests/test_sites/test_xrsite.py
+++ b/py_wake/tests/test_sites/test_xrsite.py
@@ -16,6 +16,7 @@ from py_wake.superposition_models import LinearSum
 from py_wake.tests.check_speed import timeit
 from py_wake.site._site import LocalWind
 from py_wake.utils import weibull
+from py_wake.deficit_models.noj import NOJ
 
 
 f = [0.035972, 0.039487, 0.051674, 0.070002, 0.083645, 0.064348,
@@ -244,6 +245,36 @@ def test_wd_independent_site():
     npt.assert_equal(site.ds.sector_width, 360)
 
 
+def test_i_dependent_WS():
+    ds = xr.Dataset(
+        data_vars={'WS': ('i', [8, 9, 10]), 'P': ('wd', f)},
+        coords={'wd': np.linspace(0, 360, len(f), endpoint=False)})
+    site = XRSite(ds)
+    lw = site.local_wind([0, 200, 400], [0, 0, 0], [70, 70, 70], wd=0, ws=10)
+    npt.assert_array_equal(lw.WS, [8, 9, 10])
+
+    WS = np.arange(6).reshape(3, 2) + 9
+    ds = xr.Dataset(
+        data_vars={'WS': (('i', 'ws'), WS), 'Sector_frequency': ('wd', f),
+                   'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k)},
+        coords={'wd': np.linspace(0, 360, len(f), endpoint=False), 'ws': [9, 10], 'i': [0, 1, 2]})
+    site = XRSite(ds)
+    lw = site.local_wind([0, 200, 400], [0, 0, 0], [70, 70, 70], wd=0, ws=10)
+    npt.assert_array_equal(lw.WS.squeeze(), [10, 12, 14])
+
+
+def test_i_time_dependent_WS():
+    t = np.arange(4)
+    WS_it = t[na] / 10 + np.array([9, 10])[:, na]
+    ds = xr.Dataset(
+        data_vars={'WS': (('i', 'time'), WS_it), 'P': ('wd', f), 'TI': 0.1},
+        coords={'wd': np.linspace(0, 360, len(f), endpoint=False)})
+    site = XRSite(ds)
+    wfm = NOJ(site, V80())
+    sim_res = wfm([0, 200], [0, 0], ws=WS_it.mean(0), wd=np.zeros(4), time=t)
+    npt.assert_array_equal(sim_res.WS, WS_it)
+
+
 def test_load_save(complex_grid_site):
     complex_grid_site.save(tfp + "tmp.nc")
     site = XRSite.load(tfp + "tmp.nc", interp_method='linear')
@@ -334,7 +365,7 @@ def test_neighbour_farm_speed():
     neighbour_x, neighbour_y = wt_x - 4000, wt_y
     all_x, all_y = np.r_[wt_x, neighbour_x], np.r_[wt_y, neighbour_y]
 
-    windTurbines = WindTurbines.from_WindTurbines([IEA37_WindTurbines(), IEA37_WindTurbines()])
+    windTurbines = WindTurbines.from_WindTurbine_lst([IEA37_WindTurbines(), IEA37_WindTurbines()])
     windTurbines._names = ["Current wind farm", "Neighbour wind farm"]
     types = [0] * len(wt_x) + [1] * len(neighbour_x)
 
diff --git a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
index f4ffcc9215ec176ca4038e238f4ba8a53d2f307a..8c6ac2c4b7f2a6e08015bc13544747fbecbe8cd3 100644
--- a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
+++ b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
@@ -352,7 +352,37 @@ def test_time_series_dates():
     npt.assert_array_equal(sim_res.time, t)
 
 
-def test_time_series_override_ti():
+def test_time_series_override_WS():
+    d = np.load(os.path.dirname(examples.__file__) + "/data/time_series.npz")
+    wd, ws = [d[k][:6 * 24] for k in ['wd', 'ws']]
+    t = pd.date_range("2000-01-01", freq="10T", periods=24 * 6)
+    WS_it = (np.arange(80) / 100)[:, na] + ws[na]
+    wt = V80()
+    site = Hornsrev1Site()
+    x, y = site.initial_position.T
+    wfm = NOJ(site, wt)
+    sim_res = wfm(x, y, ws=ws, wd=wd, time=t, WS=WS_it, verbose=False)
+    npt.assert_array_equal(sim_res.WS, WS_it)
+    npt.assert_array_equal(sim_res.WD, wd)
+    npt.assert_array_equal(sim_res.time, t)
+
+
+def test_time_series_override_WD():
+    d = np.load(os.path.dirname(examples.__file__) + "/data/time_series.npz")
+    wd, ws = [d[k][:6 * 24] for k in ['wd', 'ws']]
+    t = pd.date_range("2000-01-01", freq="10T", periods=24 * 6)
+    WD_it = (np.arange(80) / 100)[:, na] + wd[na]
+    wt = V80()
+    site = Hornsrev1Site()
+    x, y = site.initial_position.T
+    wfm = NOJ(site, wt)
+    sim_res = wfm(x, y, ws=ws, wd=wd, time=t, WD=WD_it, verbose=False)
+    npt.assert_array_equal(sim_res.WS, ws)
+    npt.assert_array_equal(sim_res.WD, WD_it)
+    npt.assert_array_equal(sim_res.time, t)
+
+
+def test_time_series_override_TI():
 
     d = np.load(os.path.dirname(examples.__file__) + "/data/time_series.npz")
     wd, ws, ws_std = [d[k][:6 * 24] for k in ['wd', 'ws', 'ws_std']]
diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index 82278dc99767f2c85b2c8a337e222fa5f51d00b7..98e2dd75199cd57ba9be4fe8d15d4afb8d40d054 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -163,8 +163,9 @@ class EngineeringWindFarmModel(WindFarmModel):
         self._validate_input(x_i, y_i)
 
         I, L, K, = len(x_i), len(wd), (1, len(ws))[time is False]
-        if 'TI' in kwargs:
-            lw.add_ilk('TI', kwargs['TI'])
+        for v in ['WS', 'WD', 'TI']:
+            if v in kwargs:
+                lw.add_ilk(v, kwargs[v])
 
         WS_eff_ilk = lw.WS.ilk((I, L, K)).copy()
         TI_eff_ilk = lw.TI.ilk((I, L, K)).copy()
@@ -176,7 +177,7 @@ class EngineeringWindFarmModel(WindFarmModel):
         # cw_iil = np.sqrt(hcw_iil**2 + dh_iil**2 + eps)
         wt_kwargs = kwargs
         ri, oi = self.windTurbines.function_inputs
-        unused_inputs = set(wt_kwargs) - set(ri) - set(oi) - {'TI'}
+        unused_inputs = set(wt_kwargs) - set(ri) - set(oi) - {'WS', 'WD', 'TI'}
         if unused_inputs:
             raise TypeError("""got unexpected keyword argument(s): '%s'
             required arguments: %s