_site.py

import numpy as np
from abc import ABC, abstractmethod

"""
suffixs:
- i: Local point (wind turbines)
- j: Local point (downstream turbines or positions)
- l: Wind directions
- k: Wind speeds
- m: Height above ground
"""
from numpy import newaxis as na
from scipy import interpolate
from py_wake.site.distance import StraightDistance
from py_wake.site.shear import NoShear, PowerShear


class LocalWind():
    def __init__(self, WD_ilk, WS_ilk, TI_ilk, P_ilk):
        """

        Parameters
        ----------
        WD_ilk : array_like
            local free flow wind directions
        WS_ilk : array_like
            local free flow wind speeds
        TI_ilk : array_like
            local free flow turbulence intensity
        P_ilk : array_like
            Probability/weight
        """
        self.WD_ilk = WD_ilk
        self.WS_ilk = WS_ilk
        self.TI_ilk = TI_ilk
        self.P_ilk = P_ilk


class Site(ABC):
    def __init__(self, distance):
        self.distance = distance
        self.default_ws = np.arange(3, 26.)
        self.default_wd = np.arange(360)

    def get_defaults(self, wd=None, ws=None):
        if wd is None:
            wd = self.default_wd
        else:
            wd = np.atleast_1d(wd)
        if ws is None:
            ws = self.default_ws
        else:
            ws = np.atleast_1d(ws)
        return wd, ws

    @abstractmethod
    def local_wind(self, x_i, y_i, h_i=None, wd=None, ws=None, wd_bin_size=None, ws_bins=None):
        """Local free flow wind conditions

        Parameters
        ----------
        x_i  :  array_like
            Local x coordinate
        y_i : array_like
            Local y coordinate
        h_i : array_like, optional
            Local h coordinate, i.e., heights above ground
        wd : float, int or array_like, optional
            Global wind direction(s). Override self.default_wd
        ws : float, int or array_like, optional
            Global wind speed(s). Override self.default_ws
        wd_bin_size : int or float, optional
            Size of wind direction bins. default is size between first and
            second element in default_wd
        ws_bin : array_like or None, optional
            Wind speed bin edges

        Returns
        -------
        LocalWind object containing:
            WD_ilk : array_like
                local free flow wind directions
            WS_ilk : array_like
                local free flow wind speeds
            TI_ilk : array_like
                local free flow turbulence intensity
            P_ilk : array_like
                Probability/weight
        """

    @abstractmethod
    def probability(self, x_i, y_i, h_i, WD_ilk, WS_ilk, wd_bin_size, ws_bins):
        """Probability of wind situation (wind speed and direction)

        Parameters
        ----------
        x_i : array_like
            Local x coordinate
        y_i : array_like
            Local y coordinate
        h_i : array_like
            Local height
        WD_lk : array_like
            Wind direction
        WS_lk : array_like
            Wind speed
        wd_bin_size : int or float
            size of wind direction sectors
        ws_bins : array_like
            ws bin edges, size=k+1

        Returns
        -------
        P_ilk : float or array_like
            Probability of wind speed and direction at local positions
        """

    def distances(self, src_x_i, src_y_i, src_h_i, dst_x_j, dst_y_j, dst_h_j, wd_il):
        """Calculate down/crosswind distance between source and destination points

        Parameters
        ----------
        src_x_i : array_like
            Source x position
        src_y_i : array_like
            Source y position
        src_h_i : array_like
            Source height above ground level
        dst_x_j : array_like
            Destination x position
        dst_y_j : array_like
            Destination y position
        dst_h_j : array_like
            Destination height above ground level
        wd_il : array_like, shape (#src, #wd)
            Local wind direction at the source points for all global wind directions


        Returns
        -------
        dw_ijl : array_like
            down wind distances. Positive downstream
        hcw_ijl : array_like
            horizontal cross wind distances. Positive when the wind is in your
            back and  the turbine lies on the left.
        dh_ijl : array_like
            vertical distances
        dw_order_indices_l : array_like
            indices that gives the downwind order of source points
        """
        return self.distance(self, src_x_i, src_y_i, src_h_i, dst_x_j, dst_y_j, dst_h_j, wd_il)

    def wt2wt_distances(self, x_i, y_i, h_i, wd_il):
        return self.distances(x_i, y_i, h_i, x_i, y_i, h_i, wd_il)

    @abstractmethod
    def elevation(self, x_i, y_i):
        """Local terrain elevation (height above mean sea level)

        Parameters
        ----------
        x_i : array_like
            Local x coordinate
        y_i : array_like
            Local y coordinate

        Returns
        -------
        elevation : array_like
        """

    def wd_bin_size(self, wd, wd_bin_size=None):
        if wd_bin_size is not None:
            return wd_bin_size
        else:
            return 1

    def ws_bins(self, ws, ws_bins=None):
        ws = np.asarray(ws)
        if hasattr(ws_bins, '__len__') and len(ws_bins) == len(ws) + 1:
            return ws_bins
        if len(ws.shape) and ws.shape[-1] > 1:
            d = np.diff(ws) / 2
            return np.maximum(np.concatenate(
                [ws[..., :1] - d[..., :1], ws[..., :-1] + d, ws[..., -1:] + d[..., -1:]], -1), 0)
        else:
            # ws is single value
            if ws_bins is None:
                ws_bins = 1
            return ws + np.array([-ws_bins / 2, ws_bins / 2])

    def plot_ws_distribution(self, x=0, y=0, h=70, wd=[0], include_wd_distribution=False, ax=None):
        """Plot wind speed distribution

        Parameters
        ----------
        x : int or float
            Local x coordinate
        y : int or float
            Local y coordinate
        h : int or float
            Local height above ground
        wd : int or array_like
            Wind direction(s) (one curve pr wind direction)
        include_wwd_distributeion : bool, default is False
            If true, the wind speed probability distributions are multiplied by
            the wind direction probability. The sector size is set to 360 / len(wd).
            This only makes sense if the wd array is evenly distributed
        ax : pyplot or matplotlib axes object, default None

        """
        if ax is None:
            import matplotlib.pyplot as plt
            ax = plt
        ws = np.arange(0.05, 30.05, .1)
        x = np.atleast_1d(x)
        y = np.atleast_1d(y)
        h = np.atleast_1d(h)

        wd = np.atleast_1d(wd)
        for wd_ in wd:
            wd_bin_size = 360 / len(wd)

            if include_wd_distribution:
                v = wd_bin_size / 2
                wd_l = np.arange(wd_ - v, wd_ + v) % 360
                WD_lk, WS_lk = np.meshgrid(wd_l, ws, indexing='ij')

                p = self.probability(x, y, h, WD_lk[na], WS_lk[na], wd_bin_size=1,
                                     ws_bins=self.ws_bins(WS_lk[na])).sum((0, 1))
                lbl = r"Wind direction: %d$\pm$%s deg" % (wd_, (int(v), v)[(wd_bin_size % 2) != 0])
            else:
                #                 WD_lk = np.array([wd_])[:, na]
                #                 WS_lk = ws
                WD_lk, WS_lk = np.meshgrid([wd_], ws, indexing='ij')

                p = self.probability(x, y, h, WD_lk[na, :, :], WS_lk[na, :, :],
                                     wd_bin_size=wd_bin_size, ws_bins=self.ws_bins(WS_lk[na, :, :]))[0, 0]
                p /= p.sum()
                lbl = "Wind direction: %d deg" % (wd_)

            ax.plot(ws, p * 10, label=lbl)
            ax.xlabel('Wind speed [m/s]')
            ax.ylabel('Probability')
        ax.legend(loc=1)
        return p

    def plot_wd_distribution(self, x=0, y=0, h=70, n_wd=12, ws_bins=None, ax=None):
        """Plot wind direction (and speed) distribution

        Parameters
        ----------
        x : int or float
            Local x coordinate
        y : int or float
            Local y coordinate
        h : int or float
            Local height above ground
        n_wd : int
            Number of wind direction sectors
        ws_bins : None, int or array_like, default is None
            Splits the wind direction sector pies into different colors to show
            the probability of different wind speeds\n
            If int, number of wind speed bins in the range 0-30\n
            If array_like, limits of the wind speed bins limited by ws_bins,
            e.g. [0,10,20], will show 0-10 m/s and 10-20 m/s
        ax : pyplot or matplotlib axes object, default None
        """
        if ax is None:
            import matplotlib.pyplot as plt
            ax = plt
        x = np.atleast_1d(x)
        y = np.atleast_1d(y)
        h = np.atleast_1d(h)
        wd = np.linspace(0, 360, n_wd, endpoint=False)
        theta = wd / 180 * np.pi
        if not ax.__class__.__name__ == 'PolarAxesSubplot':
            if hasattr(ax, 'subplot'):
                ax.clf()
                ax = ax.subplot(111, projection='polar')
            else:
                ax.figure.clf()
                ax = ax.figure.add_subplot(111, projection='polar')
        ax.set_theta_direction(-1)
        ax.set_theta_offset(np.pi / 2.0)

        s = 360 / n_wd
        if ws_bins is None:

            WD_lk, WS_lk = np.meshgrid(np.arange(-s / 2, s / 2) + 1, [100], indexing='ij')

            p = [self.probability(x_i=x, y_i=y, h_i=h, WD_ilk=(WD_lk[na] + wd_) % 360,
                                  WS_ilk=WS_lk[na],
                                  wd_bin_size=1, ws_bins=[0, 200]).sum() for wd_ in wd]
            ax.bar(theta, p, width=s / 180 * np.pi, bottom=0.0)
        else:
            if not hasattr(ws_bins, '__len__'):
                ws_bins = np.linspace(0, 30, ws_bins)
            else:
                ws_bins = np.asarray(ws_bins)
            ws = ((ws_bins[1:] + ws_bins[:-1]) / 2)
            ws_bins = self.ws_bins(ws)

            WD_lk, WS_lk = np.meshgrid(np.arange(-s / 2, s / 2) + 1, ws, indexing='ij')
            p = [self.probability(x_i=x, y_i=y, h_i=h, WD_ilk=(WD_lk[na] + wd_) % 360, WS_ilk=WS_lk[na],
                                  wd_bin_size=1, ws_bins=ws_bins).sum((0, 1)) for wd_ in wd]
            cum_p = np.cumsum(p, 1).T
            start_p = np.vstack([np.zeros_like(cum_p[:1]), cum_p[:-1]])

            for ws1, ws2, p_ws1, p_ws2 in zip(ws_bins[:-1], ws_bins[1:], start_p, cum_p):
                ax.bar(theta, p_ws2 - p_ws1, width=s / 180 * np.pi, bottom=p_ws1, label="%s-%s m/s" % (ws1, ws2))
            ax.legend(bbox_to_anchor=(1.15, 1.1))

        ax.set_rlabel_position(-22.5)  # Move radial labels away from plotted line
        ax.grid(True)
        return p


class UniformSite(Site):
    """Site with uniform (same wind over all, i.e. flat uniform terrain) and
    constant wind speed probability of 1. Only for one fixed wind speed
    """

    def __init__(self, p_wd, ti, ws=12, interp_method='piecewise', shear=NoShear()):
        super().__init__(StraightDistance())
        self.default_ws = ws
        self.ti = Sector2Subsector(np.atleast_1d(ti), interp_method=interp_method)
        self.p_wd = Sector2Subsector(p_wd / np.sum(p_wd), interp_method=interp_method) / (360 / len(p_wd))
        self.shear = shear

    def probability(self, x_i, y_i, h_i, WD_ilk, WS_ilk, wd_bin_size, ws_bins):
        P_lk = np.ones_like(WS_ilk[0], dtype=np.float) * \
            self.p_wd[np.round(WD_ilk[0]).astype(int) % 360] * wd_bin_size
        return P_lk[na]

    def local_wind(self, x_i=None, y_i=None, h_i=None, wd=None, ws=None, wd_bin_size=None, ws_bins=None):
        if wd is None:
            wd = self.default_wd
        if ws is None:
            ws = self.default_ws

        ws_bins = self.ws_bins(ws, ws_bins)
        wd_bin_size = self.wd_bin_size(wd, wd_bin_size)
        WD_ilk, WS_ilk = [np.tile(W, (len(x_i), 1, 1)).astype(np.float)
                          for W in np.meshgrid(wd, ws, indexing='ij')]
        WD_index_ilk = np.round(WD_ilk).astype(int)
        if h_i is not None:
            WS_ilk = self.shear(WS_ilk, WD_ilk, h_i)

        TI_ilk = self.ti[WD_index_ilk]

        P_ilk = self.probability(0, 0, 0, WD_ilk, WS_ilk, wd_bin_size, ws_bins)
        return LocalWind(WD_ilk, WS_ilk, TI_ilk, P_ilk)

    def elevation(self, x_i, y_i):
        return np.zeros_like(x_i)


class UniformWeibullSite(UniformSite):
    """Site with uniform (same wind over all, i.e. flat uniform terrain) and
    weibull distributed wind speed
    """

    def __init__(self, p_wd, a, k, ti, interp_method='nearest', shear=NoShear()):
        """Initialize UniformWeibullSite

        Parameters
        ----------
        p_wd : array_like
            Probability of wind direction sectors
        a : array_like
            Weilbull scaling parameter of wind direction sectors
        k : array_like
            Weibull shape parameter
        ti : float or array_like
            Turbulence intensity
        interp_method : 'nearest', 'linear' or 'spline'
            p_wd, a, k, ti and alpha are interpolated to 1 deg sectors using this
            method
        shear : Shear object
            Shear object, e.g. NoShear(), PowerShear(h_ref, alpha)

        Notes
        ------
        The wind direction sectors will be: [0 +/- w/2, w +/- w/2, ...]
        where w is 360 / len(p_wd)

        """
        super().__init__(p_wd, ti, interp_method=interp_method, shear=shear)
        self.default_ws = np.arange(3, 26)
        self.a = Sector2Subsector(a, interp_method=interp_method)
        self.k = Sector2Subsector(k, interp_method=interp_method)

    def weibull_weight(self, WS, A, k, ws_bins):
        def cdf(ws, A=A, k=k):
            return 1 - np.exp(-(ws / A) ** k)
        ws_bins = np.asarray(ws_bins)
        return cdf(ws_bins[..., 1:]) - cdf(ws_bins[..., :-1])

    def probability(self, x_i, y_i, h_i, WD_ilk, WS_ilk, wd_bin_size, ws_bins):
        i_ilk = np.round(WD_ilk).astype(int) % 360
        p_wd_ilk = self.p_wd[i_ilk] * wd_bin_size
        if wd_bin_size == 360:
            p_wd_ilk[:] = 1
        P_ilk = self.weibull_weight(WS_ilk, self.a[i_ilk], self.k[i_ilk], ws_bins) * p_wd_ilk
        return P_ilk


def Sector2Subsector(para, axis=-1, wd_binned=None, interp_method='piecewise'):
    """ Expand para on the wind direction dimension, i.e., increase the nubmer
    of sectors (sectors to subsectors), by interpolating between sectors, using
    specified method.

    Parameters
    ----------
    para : array_like
        Parameter to be expand, it can be sector-wise Weibull A, k, frequency.
    axis : integer
        Denotes which dimension of para corresponds to wind direction.
    wd_binned : array_like
        Wind direction of subsectors to be expanded to.
    inter_method : string
        'piecewise'/'linear'/'spline', based on interp1d in scipy.interpolate,
        'spline' means cubic spline.

    --------------------------------------
    Note: the interpolating method for sector-wise Weibull distributions and
    joint distribution of wind speed and wind direction is referred to the
    following paper:
        Feng, J. and Shen, W.Z., 2015. Modelling wind for wind farm layout
        optimization using joint distribution of wind speed and wind direction.
        Energies, 8(4), pp.3075-3092. [https://doi.org/10.3390/en8043075]
    """
    if wd_binned is None:
        wd_binned = np.arange(360)
    para = np.array(para)
    num_sector = para.shape[axis]
    wd_sector = np.linspace(0, 360, num_sector, endpoint=False)

    try:
        interp_index = ['nearest', 'piecewise', 'linear', 'spline'].index(interp_method)
        interp_kind = ['nearest', 'nearest', 'linear', 'cubic'][interp_index]
    except ValueError:
        raise NotImplementedError(
            'interp_method={0} not implemeted yet.'.format(interp_method))
    wd_sector_extended = np.hstack((wd_sector, 360.0))
    para_sector_extended = np.concatenate((para, para.take([0], axis=axis)),
                                          axis=axis)
    if interp_kind == 'cubic' and len(wd_sector_extended) < 4:
        interp_kind = 'linear'
    f_interp = interpolate.interp1d(wd_sector_extended, para_sector_extended,
                                    kind=interp_kind, axis=axis)
    para_expanded = f_interp(wd_binned % 360)

    return para_expanded


def main():
    if __name__ == '__main__':
        f = [0.035972, 0.039487, 0.051674, 0.070002, 0.083645, 0.064348,
             0.086432, 0.117705, 0.151576, 0.147379, 0.10012, 0.05166]
        A = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921,
             9.584007, 10.51499, 11.39895, 11.68746, 11.63732, 10.08803]
        k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703,
             2.583984, 2.548828, 2.470703, 2.607422, 2.626953, 2.326172]
        ti = .1
        h_ref = 100
        alpha = .1
        site = UniformWeibullSite(f, A, k, ti, shear=PowerShear(h_ref=h_ref, alpha=alpha))

        x_i = y_i = np.arange(5)
        wdir_lst = np.arange(0, 360, 90)
        wsp_lst = np.arange(1, 20)
        local_wind = site.local_wind(x_i=x_i, y_i=y_i, wd=wdir_lst, ws=wsp_lst)
        print(local_wind.WS_ilk.shape)
        import matplotlib.pyplot as plt

        site.plot_ws_distribution(0, 0, wdir_lst)

        plt.figure()
        z = np.arange(1, 100)
        u = [site.local_wind(x_i=[0], y_i=[0], h_i=[z_], wd=0, ws=10).WS_ilk[0][0] for z_ in z]
        plt.plot(u, z)
        plt.xlabel('Wind speed [m/s]')
        plt.ylabel('Height [m]')
        plt.show()


main()