fix and speedup flowmap

py_wake/flow_map.py

@@ -20,12 +20,12 @@ class FlowBox(xr.Dataset):
             coords = localWind_j.coords
             X = localWind_j.i
         else:
-            coords = {'x': X[0, :, 0], 'y': Y[:, 0, 0], 'h': H[0, 0, :],
-                      **{k: (dep, v, {'Description': d}) for k, dep, v, d in [
-                          ('wd', ('wd', 'time')[time], wd, 'Ambient reference wind direction [deg]'),
-                          ('ws', ('ws', 'time')[time], ws, 'Ambient reference wind speed [m/s]')]}}
+            coords = {'x': X[0, :, 0], 'y': Y[:, 0, 0], 'h': H[0, 0, :]}
             if time:
                 coords['time'] = lw_j.coords['time']
+        coords.update({k: (dep, v, {'Description': d}) for k, dep, v, d in [
+                          ('wd', ('wd', 'time')[time], wd, 'Ambient reference wind direction [deg]'),
+                          ('ws', ('ws', 'time')[time], ws, 'Ambient reference wind speed [m/s]')]})
 
         def get_da(arr_jlk):
             if len(X.shape) == 1:
@@ -161,7 +161,6 @@ class FlowMap(FlowBox):
             if True (default), lines/circles showing the wind turbine rotors are plotted
         ax : pyplot or matplotlib axes object, default None
         """
-        import matplotlib.pyplot as plt
         if cmap is None:
             cmap = 'Blues_r'
         if ax is None:
@@ -186,12 +185,16 @@ class FlowMap(FlowBox):
             x = np.zeros_like(y) + self.plane[1]
             z = self.simulationResult.windFarmModel.site.elevation(x, y)
             ax.plot(y / n, z / n, 'k')
-        else:
+        elif self.plane[0] == 'XY':
+
             # xarray gives strange levels
             # c = data.isel(h=0).plot(levels=levels, cmap=cmap, ax=ax, add_colorbar=plot_colorbar)
             c = ax.contourf(self.X / n, self.Y / n, data.isel(h=0).data, levels=levels, cmap=cmap)
             if plot_colorbar:
                 plt.colorbar(c, label=clabel, ax=ax)
+        else:
+            raise NotImplementedError(
+                f"Plot not supported for FlowMaps based on Points. Use XYGrid, YZGrid or XZGrid instead")
 
         if plot_windturbines:
             self.plot_windturbines(normalize_with=normalize_with, ax=ax)
@@ -245,10 +248,9 @@ class FlowMap(FlowBox):
             if True (default), lines/circles showing the wind turbine rotors are plotted
         ax : pyplot or matplotlib axes object, default None
         """
-        if 'time' in self:
-            WS_eff = (self.WS_eff * self.P / self.P.sum(['time'])).sum(['time'])
-        else:
-            WS_eff = (self.WS_eff * self.P / self.P.sum(['wd', 'ws'])).sum(['wd', 'ws'])
+        sum_dims = [d for d in ['wd', 'time', 'ws'] if d in self.P.dims]
+        WS_eff = (self.WS_eff * self.P / self.P.sum(sum_dims)).sum(sum_dims)
+
         return self.plot(WS_eff, clabel='wind speed [m/s]',
                          levels=levels, cmap=cmap, plot_colorbar=plot_colorbar,
                          plot_windturbines=plot_windturbines, normalize_with=normalize_with, ax=ax)

py_wake/tests/test_flow_map.py

@@ -369,9 +369,10 @@ def test_wd_dependent_flow_map():
     sim_res = wfm(x=[0], y=[0], wd=[0, 90, 180])
     for wd in [[0], [0, 90], None]:
         fm = sim_res.flow_map(wd=wd)
+        fm.plot_wake_map()
         if 0:
-            fm.plot_wake_map()
             plt.show()
+        plt.close('all')
 
 
 def test_ws_dependent_flow_map():
@@ -379,9 +380,10 @@ def test_ws_dependent_flow_map():
     sim_res = wfm(x=[0], y=[0], ws=[8, 9, 10], wd=270)
     for ws in [[8], [8, 9], None]:
         fm = sim_res.flow_map(ws=ws)
+        fm.plot_wake_map()
         if 0:
-            fm.plot_wake_map()
             plt.show()
+        plt.close('all')
 
 
 def test_time_dependent_flow_map():
@@ -389,6 +391,16 @@ def test_time_dependent_flow_map():
     sim_res = wfm(x=[0], y=[0], wd=[0, 90, 180], ws=[8, 9, 10], time=True)
     for t in [[0], [0, 1], None]:
         fm = sim_res.flow_map(time=t)
+        fm.plot_wake_map()
         if 0:
-            fm.plot_wake_map()
             plt.show()
+        plt.close('all')
+
+
+def test_i_dependent_flow_map():
+    wfm = IEA37CaseStudy1(16)
+    sim_res = wfm(x=[0], y=[0], wd=[0, 90, 180])
+    X, Y = np.meshgrid(np.linspace(-2000, 2000, 50), np.linspace(-2000, 2000, 50))
+    fm = sim_res.flow_map(Points(x=X.flatten(), y=Y.flatten(), h=X.flatten() * 0 + 110))
+    with pytest.raises(NotImplementedError, match="Plot not supported for FlowMaps based on Points. Use XYGrid, YZGrid or XZGrid instead"):
+        fm.plot_wake_map()

py_wake/wind_farm_models/engineering_models.py

@@ -325,7 +325,7 @@ class EngineeringWindFarmModel(WindFarmModel):
                 'D_dst_ijl': lambda l: np.zeros((I, J, 1)),
                 'h_il': lambda l: wt_h_ilk[:, :, 0],
                 'ct_ilk': get_ilk('CT'),
-                'IJLK': lambda l=wd: (I, J, len(np.atleast_1d(l)), K)}, lw_j, wd, WD_il
+                'IJLK': lambda l=slice(None), I=I, J=J, L=L, K=K: (I, J, len(np.arange(L)[l]), K)}, lw_j, wd, WD_il
 
     def _get_flow_l(self, model_kwargs, l, wt_x_ilk, wt_y_ilk, wt_h_ilk, lw_j, wd, WD_ilk):
 
@@ -399,7 +399,10 @@ class EngineeringWindFarmModel(WindFarmModel):
             return (lw_j, np.broadcast_to(lw_j.WS_ilk, (len(x_j), L, K)).astype(float),
                     np.broadcast_to(lw_j.TI_ilk, (len(x_j), L, K)).astype(float))
 
-        l_iter = tqdm([slice(l, l + 1) for l in range(L)], disable=L <= 1 or not self.verbose,
+        size_gb = I * J * L * K * 8 / 1024**3
+        wd_chunks = np.minimum(np.maximum(int(size_gb // 1), 1), L)
+        wd_i = np.round(np.linspace(0, L, wd_chunks + 1)).astype(int)
+        l_iter = tqdm([slice(i0, i1) for i0, i1 in zip(wd_i[:-1], wd_i[1:])], disable=L <= 1 or not self.verbose,
                       desc='Calculate flow map', unit='wd')
         wt_x_ilk, wt_y_ilk, wt_h_ilk = [sim_res_data[k].ilk() for k in ['x', 'y', 'h']]
         WS_eff_jlk, TI_eff_jlk = zip(*[self._get_flow_l(