flow_map.py

from py_wake import np
import xarray as xr
from numpy import newaxis as na
from scipy.interpolate import InterpolatedUnivariateSpline
import matplotlib.pyplot as plt
from py_wake.utils.xarray_utils import ilk2da


class FlowBox(xr.Dataset):
    __slots__ = ('simulationResult', 'windFarmModel')

    def __init__(self, simulationResult, X, Y, H, localWind_j, WS_eff_jlk, TI_eff_jlk):
        self.simulationResult = simulationResult
        self.windFarmModel = self.simulationResult.windFarmModel
        lw_j = localWind_j
        wd, ws = lw_j.wd, lw_j.ws
        time = 'time' in simulationResult

        if X is None and Y is None and H is None:
            coords = localWind_j.coords
            X = localWind_j.i
        else:
            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:
                return ilk2da(arr_jlk, coords)
            else:
                shape = [X.shape + (len(wd), len(ws)), X.shape + (len(wd), )][time]
                dims = ['y', 'x', 'h'] + [['wd', 'ws'], ['time']][time]
                return xr.DataArray(arr_jlk.reshape(shape), coords, dims)
        JLK = WS_eff_jlk.shape
        xr.Dataset.__init__(self, data_vars={k: get_da(np.broadcast_to(v, JLK)) for k, v in [
            ('WS_eff', WS_eff_jlk), ('TI_eff', TI_eff_jlk),
            ('WD', lw_j.WD_ilk), ('WS', lw_j.WS_ilk), ('TI', lw_j.TI_ilk), ('P', lw_j.P_ilk)]})


class FlowMap(FlowBox):
    __slots__ = ('simulationResult', 'windFarmModel', 'X', 'Y', 'plane', 'WS_eff_xylk', 'TI_eff_xylk')

    def __init__(self, simulationResult, X, Y, localWind_j, WS_eff_jlk, TI_eff_jlk, plane):
        self.X = X
        self.Y = Y
        self.plane = plane

        if plane[0] == 'XY':
            X = X[:, :, na]
            Y = Y[:, :, na]
            H = np.reshape(localWind_j.h, X.shape)
        elif plane[0] == 'YZ':
            H = Y.T[:, na, :]
            Y = X.T[:, na, :]
            X = np.reshape(localWind_j.x, Y.shape)
        elif plane[0] == 'XZ':
            H = Y.T[:, na, :]
            X = X.T[na, :, :]
            Y = np.reshape(localWind_j.y, X.shape)
        elif plane[0] == 'xyz':
            X = None
            Y = None
            H = None
        else:
            raise NotImplementedError()
        FlowBox.__init__(self, simulationResult, X, Y, H, localWind_j, WS_eff_jlk, TI_eff_jlk)

        # set flowMap.WS_xylk etc.
        if plane[0] == "XY":
            for k in ['WS_eff', 'TI_eff', 'WS', 'WD', 'TI', 'P']:
                setattr(self.__class__, "%s_xylk" % k, property(lambda self, k=k: self[k].isel(h=0)))
        elif plane[0] == "YZ":
            for k in ['WS_eff', 'TI_eff', 'WS', 'WD', 'TI', 'P']:
                self[k] = self[k].transpose('h', 'y', ...)
                setattr(self.__class__, "%s_xylk" % k,
                        property(lambda self, k=k: self[k].isel(x=0).transpose('y', 'h', ...)))
        elif plane[0] == "XZ":
            for k in ['WS_eff', 'TI_eff', 'WS', 'WD', 'TI', 'P']:
                self[k] = self[k].transpose('h', 'x', ...)
                setattr(self.__class__, "%s_xylk" % k,
                        property(lambda self, k=k: self[k].isel(x=0).transpose('x', 'h', ...)))
        elif plane[0] == "xyz":
            for k in ['WS_eff', 'TI_eff', 'WS', 'WD', 'TI', 'P']:
                setattr(self.__class__, "%s_xylk" % k, property(lambda self, k=k: self[k].ilk()))

    @property
    def XY(self):
        return self.X, self.Y

    def power_xylk(self, with_wake_loss=True, **wt_kwargs):
        if with_wake_loss:
            ws = self.WS_eff_xylk

        else:
            ws = self.WS_xylk

        power_xylk = self.windFarmModel.windTurbines.power(ws, **wt_kwargs)
        if self.plane[0] == "xyz":
            return power_xylk
        else:
            return xr.DataArray(power_xylk[:, :, na], self.coords, dims=['y', 'x', 'h', 'wd', 'ws'])

    def aep_xylk(self, normalize_probabilities=False, with_wake_loss=True, **wt_kwargs):
        """Anual Energy Production of a potential wind turbine at all grid positions (x,y)
        for all wind directions (l) and wind speeds (k) in GWh.

        Parameters
        ----------
        normalize_propabilities : Optional bool, defaults to False
            In case only a subset of all wind speeds and/or wind directions is simulated,
            this parameter determines whether the returned AEP represents the energy produced in the fraction
            of a year where these flow cases occurs or a whole year of northern wind.
            If for example, wd=[0], then
            - False means that the AEP only includes energy from the faction of year\n
            with northern wind (359.5-0.5deg), i.e. no power is produced the rest of the year.
            - True means that the AEP represents a whole year of northen wind.
            default is False
        with_wake_loss : Optional bool, defaults to True
            If True, wake loss is included, i.e. power is calculated using local effective wind speed\n
            If False, wake loss is neglected, i.e. power is calculated using local free flow wind speed
        wt_type : Optional arguments
            Additional required/optional arguments needed by the WindTurbines to computer power, e.g. type, Air_density
         """
        power_xylk = self.power_xylk(with_wake_loss, **wt_kwargs)
        P_xylk = self.P_xylk  # .isel.ilk((1,) + power_xylk.shape[2:])
        if normalize_probabilities:
            P_xylk = P_xylk / P_xylk.sum(['wd', 'ws'])
        return power_xylk * P_xylk * 24 * 365 * 1e-9

    def aep_xy(self, normalize_probabilities=False, with_wake_loss=True, **wt_kwargs):
        """Anual Energy Production of a potential wind turbine at all grid positions (x,y)
        (sum of all wind directions and wind speeds)  in GWh.

        see aep_xylk
        """
        return self.aep_xylk(normalize_probabilities, with_wake_loss, **wt_kwargs).sum(['wd', 'ws'])

    def plot(self, data, clabel, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
             normalize_with=1, ax=None):
        """Plot data as contouf map

        Parameters
        ----------
        data : array_like
            2D data array to plot
        clabel : str
            colorbar label
        levels : int or array-like, default 100
            Determines the number and positions of the contour lines / regions.
            If an int n, use n data intervals; i.e. draw n+1 contour lines. The level heights are automatically chosen.
            If array-like, draw contour lines at the specified levels. The values must be in increasing order.
        cmap : str or Colormap, defaults 'Blues_r'.
            A Colormap instance or registered colormap name.
            The colormap maps the level values to colors.
        plot_colorbar : bool, default True
            if True (default), colorbar is drawn
        plot_windturbines : bool, default True
            if True (default), lines/circles showing the wind turbine rotors are plotted
        ax : pyplot or matplotlib axes object, default None
        """
        if cmap is None:
            cmap = 'Blues_r'
        if ax is None:
            ax = plt.gca()
        n = normalize_with
        if self.plane[0] in ['XZ', 'YZ']:
            if self.plane[0] == "YZ":
                y = self.y.values
                x = np.zeros_like(y) + self.plane[1]
                z = self.simulationResult.windFarmModel.site.elevation(x, y)
                c = ax.contourf(self.X, self.Y + z, data.isel(x=0), levels=levels, cmap=cmap)
            elif self.plane[0] == 'XZ':
                x = self.x.values
                y = np.zeros_like(x) + self.plane[1]
                z = self.simulationResult.windFarmModel.site.elevation(x, y)
                c = ax.contourf(self.X, self.Y + z, data.isel(y=0), levels=levels, cmap=cmap)

            if plot_colorbar:
                plt.colorbar(c, label=clabel, ax=ax)
            # plot terrain
            y = np.arange(y.min(), y.max())
            x = np.zeros_like(y) + self.plane[1]
            z = self.simulationResult.windFarmModel.site.elevation(x, y)
            ax.plot(y / n, z / n, 'k')
        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)

        return c

    def plot_windturbines(self, normalize_with=1, ax=None):
        fm = self.windFarmModel

        def get(k):
            v = self.simulationResult[k]
            if 'wd' in v.dims:
                v = v.sel(wd=self.wd[0])
            if 'ws' in v.dims:
                v = v.mean(['ws'])
            return v.data

        yaw, tilt = get('yaw'), get('tilt')

        x_i, y_i = self.simulationResult.x.values, self.simulationResult.y.values
        type_i = self.simulationResult.type.data
        if self.plane[0] in ['XZ', "YZ"]:
            h_i = self.simulationResult.h.values
            z_i = self.simulationResult.windFarmModel.site.elevation(x_i, y_i)
            if self.plane[0] == 'XZ':
                fm.windTurbines.plot_yz(x_i, z_i, h_i, types=type_i, wd=self.wd - 90, yaw=yaw, tilt=tilt,
                                        normalize_with=normalize_with, ax=ax)
            else:
                fm.windTurbines.plot_yz(y_i, z_i, h_i, types=type_i, wd=self.wd, yaw=yaw, tilt=tilt,
                                        normalize_with=normalize_with, ax=ax)
        else:  # self.plane[0] == "XY":
            fm.windTurbines.plot_xy(x_i, y_i, type_i,
                                    wd=self.wd.values, yaw=yaw, tilt=tilt, normalize_with=normalize_with, ax=ax)

    def plot_wake_map(self, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
                      normalize_with=1, ax=None):
        """Plot effective wind speed contourf map

        Parameters
        ----------
        levels : int or array-like, default 100
            Determines the number and positions of the contour lines / regions.
            If an int n, use n data intervals; i.e. draw n+1 contour lines. The level heights are automatically chosen.
            If array-like, draw contour lines at the specified levels. The values must be in increasing order.
        cmap : str or Colormap, defaults 'Blues_r'.
            A Colormap instance or registered colormap name.
            The colormap maps the level values to colors.
        plot_colorbar : bool, default True
            if True (default), colorbar is drawn
        plot_windturbines : bool, default True
            if True (default), lines/circles showing the wind turbine rotors are plotted
        ax : pyplot or matplotlib axes object, default None
        """
        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)

    def plot_ti_map(self, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True, ax=None):
        """Plot effective turbulence intensity contourf map

        Parameters
        ----------
        levels : int or array-like, default 100
            Determines the number and positions of the contour lines / regions.
            If an int n, use n data intervals; i.e. draw n+1 contour lines. The level heights are automatically chosen.
            If array-like, draw contour lines at the specified levels. The values must be in increasing order.
        cmap : str or Colormap, defaults 'Blues'.
            A Colormap instance or registered colormap name.
            The colormap maps the level values to colors.
        plot_colorbar : bool, default True
            if True (default), colorbar is drawn
        plot_windturbines : bool, default True
            if True (default), lines/circles showing the wind turbine rotors are plotted
        ax : pyplot or matplotlib axes object, default None

        """
        if cmap is None:
            cmap = 'Blues'
        c = self.plot(self.TI_eff.mean(['wd', 'ws']), clabel="Turbulence intensity [-]",
                      levels=levels, cmap=cmap, plot_colorbar=plot_colorbar,
                      plot_windturbines=plot_windturbines, ax=ax)

        return c

    def min_WS_eff(self, x=None, h=None):
        if x is None:
            x = self.x
        if h is None:
            h = self.h[0].item()
        WS_eff = self.WS_eff.sel_interp_all(xr.Dataset(coords={'x': x, 'h': h}))
        y = WS_eff.y.values

        def get_min(y, v):
            i = np.argmin(v)
            s = slice(i - 3, i + 4)
            if len(v[s]) < 7 or len(np.unique(v[s])) == 1:
                return np.nan
#             import matplotlib.pyplot as plt
#             plt.plot(y, v)
#             y_ = np.linspace(y[s][0], y[s][-1], 100)
#             plt.plot(y_, InterpolatedUnivariateSpline(y[s], v[s])(y_))
#             plt.axvline(np.interp(0, InterpolatedUnivariateSpline(y[s], v[s]).derivative()(y[s]), y[s]))
#             plt.axvline(0, color='k')
#             plt.show()
            return np.interp(0, InterpolatedUnivariateSpline(y[s], v[s]).derivative()(y[s]), y[s])

        y_min_ws = [get_min(y, ws) for ws in WS_eff.squeeze(['ws', 'wd']).T.values]
        return xr.DataArray(y_min_ws, coords={'x': x, 'h': h}, dims='x')

    def plot_deflection_grid(self, normalize_with=1, ax=None):
        assert self.windFarmModel.deflectionModel is not None
        assert len(self.simulationResult.wt) == 1
        assert len(self.simulationResult.ws) == 1
        assert len(self.simulationResult.wd) == 1
        x, y = self.x, self.y
        y = y[::len(y) // 10]

        X, Y = np.meshgrid(x, y)

        from py_wake.utils.model_utils import get_model_input
        kwargs = get_model_input(self.windFarmModel, X.flatten(), Y.flatten(), ws=self.ws, wd=self.wd,
                                 yaw=self.simulationResult.yaw.ilk())
        dw, hcw, dh = self.windFarmModel.deflectionModel.calc_deflection(**kwargs)
        Yp = -hcw[0, :, 0, 0].reshape(X.shape)
        ax = ax or plt.gca()
        X, Y, Yp = [v / normalize_with for v in [X, Y, Yp]]
        # ax.plot(X[255, :], Y[255, :], 'grey', lw=3)
        for x, y, yp in zip(X, Y, Yp):
            ax.plot(x, y, 'grey', lw=1, zorder=-32)
            ax.plot(x, yp, 'k', lw=1)


class Grid():
    default_resolution = 500


class HorizontalGrid(Grid):

    def __init__(self, x=None, y=None, h=None, resolution=None, extend=.2):
        """Generate a horizontal grid for a flow map

        Parameters
        ----------
        x : array_like, optional
            x coordinates used for generating meshgrid\n
        y : array_like, optional
            y coordinates used for generating meshgrid
        h : array_like, optional
            height above ground, defaults to mean wind turbine hub height
        resolution : int or None, optional
            grid resolution if x or y is not specified. defaults to self.default_resolution
        extend : float, optional
            defines the oversize of the grid if x or y is not specified

        Notes
        -----
        if x or y is not specified then a grid with <resolution> number of points
        covering the wind turbines + <extend> x range
        """
        self.resolution = resolution or self.default_resolution
        self.x = x
        self.y = y
        self.h = h
        self.extend = extend
        self.plane = "XY", h

    def __call__(self, x_i, y_i, h_i, **_):
        # setup horizontal X,Y grid
        def f(x, N=self.resolution, ext=self.extend):
            ext *= np.max([1000, (np.max(x) - np.min(x))])
            return np.linspace(np.min(x) - ext, np.max(x) + ext, N)
        x, y, h = self.x, self.y, self.h
        if x is None:
            x = f(x_i)
        if y is None:
            y = f(y_i)
        if self.h is None:
            h = np.mean(h_i)
        else:
            h = self.h
        self.plane = "XY", h

        X, Y = np.meshgrid(x, y)
        H = np.broadcast_to(h, X.shape)
        return X, Y, X.flatten(), Y.flatten(), H.flatten()


XYGrid = HorizontalGrid


class YZGrid(Grid):

    def __init__(self, x, y=None, z=None, resolution=None, extend=.2):
        """Generate a vertical grid for a flow map in the yz-plane

        Parameters
        ----------
        x : array_like, optional
            x coordinate for the yz-grid\n
        y : array_like, optional
            y coordinates used for generating meshgrid
        z : array_like, optional
            z coordinates(height above ground) used for generating meshgrid
        resolution : int or None, optional
            grid resolution if x or y is not specified. defaults to self.default_resolution
        extend : float, optional
            defines the oversize of the grid if x or y is not specified

        Notes
        -----
        if y or z is not specified then a grid with <resolution> number of points
        covering the wind turbines + <extend> * range
        """
        self.resolution = resolution or self.default_resolution
        self.x = x
        self.y = y
        self.z = z
        self.extend = extend
        self.plane = "YZ", x

    def __call__(self, x_i, y_i, h_i, d_i):
        # setup horizontal X,Y grid
        def f(x, N=self.resolution, ext=self.extend):
            ext *= max(1000, (max(x) - min(x)))
            return np.linspace(min(x) - ext, max(x) + ext, N)
        x, y, z = self.x, self.y, self.z
        if y is None:
            y = f(y_i)
        if self.z is None:
            z = np.arange(0, (1 + self.extend) * (h_i.max() + d_i.max() / 2), np.diff(y[:2])[0])
        else:
            z = self.z

        Y, Z = np.meshgrid(y, z)
        X = np.zeros_like(Y) + x
        return Y, Z, X.T.flatten(), Y.T.flatten(), Z.T.flatten()


class XZGrid(YZGrid):
    def __init__(self, y, x=None, z=None, resolution=None, extend=.2):
        """Generate a vertical grid for a flow map in the xz-plane

        Parameters
        ----------
        y : array_like, optional
            y coordinate used for generating meshgrid
        x : array_like, optional
            x coordinatex for the yz-grid\n
        z : array_like, optional
            z coordinates(height above ground) used for generating meshgrid
        resolution : int or None, optional
            grid resolution if x or y is not specified. defaults to self.default_resolution
        extend : float, optional
            defines the oversize of the grid if x or y is not specified

        Notes
        -----
        if x or z is not specified then a grid with <resolution> number of points
        covering the wind turbines + <extend> * range
        """
        YZGrid.__init__(self, x, y=y, z=z, resolution=resolution, extend=extend)
        self.plane = "XZ", y

    def __call__(self, x_i, y_i, h_i, d_i):
        # setup horizontal X,Y grid
        def f(x, N=self.resolution, ext=self.extend):
            ext *= max(1000, (max(x) - min(x)))
            return np.linspace(min(x) - ext, max(x) + ext, N)
        x, y, z = self.x, self.y, self.z
        if x is None:
            x = f(x_i)
        if self.z is None:
            z = np.arange(0, (1 + self.extend) * (h_i.max() + d_i.max() / 2), np.diff(x[:2])[0])
        else:
            z = self.z

        X, Z = np.meshgrid(x, z)
        Y = np.zeros_like(X) + y
        return X, Z, X.T.flatten(), Y.T.flatten(), Z.T.flatten()


class Points(Grid):
    def __init__(self, x, y, h):
        assert len(x) == len(y) == len(h)
        self.x = x
        self.y = y
        self.h = h
        self.plane = 'xyz', None

    def __call__(self, **_):
        return None, None, self.x, self.y, self.h