_site.py

import matplotlib.pyplot as plt
import numpy as np
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

"""
suffixs:
- i: Local point (wind turbines)
- j: Local point (downstream turbines or positions)
- l: Wind directions
- k: Wind speeds
- m: Height above ground
"""


class LocalWind(xr.Dataset):
    __slots__ = ('wd_bin_size')

    def __init__(self, x_i, y_i, h_i, wd, ws, time, wd_bin_size, WD=None, WS=None, TI=None, P=None):
        """

        Parameters
        ----------
        WD : array_like
            local free flow wind directions
        WS : array_like
            local free flow wind speeds
        TI : array_like
            local free flow turbulence intensity
        P : array_like
            Probability/weight
        """
        ws = np.atleast_1d(ws)
        if time is not False:
            assert len(wd) == len(ws)
            if time is True:
                time = np.arange(len(wd))
            coords = {'time': time, 'wd': ('time', wd), 'ws': ('time', ws)}
        else:
            coords = {'wd': wd, 'ws': np.atleast_1d(ws)}

        assert len(np.atleast_1d(x_i)) == len(np.atleast_1d(y_i))
        n_i = max(len(np.atleast_1d(x_i)), len(np.atleast_1d(h_i)))
        coords['i'] = np.arange(n_i)

        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, 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.
        for k in ['WD', 'WS', 'TI', 'P', 'TI_std']:
            setattr(self.__class__, "%s_ilk" % k, property(lambda self, k=k: self[k].ilk()))

    def set_data_array(self, data_array, name, description):
        if data_array is not None:
            data_array.attrs.update({'Description': description})
            self[name] = data_array

    def set_W(self, ws, wd, ti, ws_bins, use_WS=False):
        for da, name, desc in [(ws, 'WS', 'Local free-stream wind speed [m/s]'),
                               (wd, 'WD', 'Local free-stream wind direction [deg]'),
                               (ti, 'TI', 'Local free-stream turbulence intensity')]:
            self.set_data_array(da, name, desc)

        # upper and lower bounds of wind speed bins
        WS = [self.ws, self.WS][use_WS]
        lattr = {'Description': 'Lower bound of wind speed bins [m/s]'}
        uattr = {'Description': 'Upper bound of wind speed bins [m/s]'}
        if not hasattr(ws_bins, '__len__') or len(ws_bins) != len(WS) + 1:
            if len(WS.shape) and WS.shape[-1] > 1:
                d = np.diff(WS) / 2
                ws_bins = 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
                ws_bins = WS.data + np.array([-ws_bins / 2, ws_bins / 2])

            self['ws_lower'] = xr.DataArray(ws_bins[..., :-1], dims=WS.dims, attrs=lattr)
            self['ws_upper'] = xr.DataArray(ws_bins[..., 1:], dims=WS.dims, attrs=uattr)
        else:
            self['ws_lower'] = xr.DataArray(ws_bins[:-1], dims=['ws'], attrs=lattr)
            self['ws_upper'] = xr.DataArray(ws_bins[1:], dims=['ws'], attrs=uattr)


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

    @property
    def distance(self):
        return self._distance

    @distance.setter
    def distance(self, distance):
        self._distance = distance
        distance.site = self

    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

    def local_wind(self, x_i, y_i, h_i=None, wd=None, ws=None, time=False, 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
        time : boolean or array_like
            If True or array_like, wd and ws is interpreted as a time series
            If False, full wd x ws matrix is computed
        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
        """
        wd, ws = self.get_defaults(wd, ws)
        wd_bin_size = self.wd_bin_size(wd, wd_bin_size)
        lw = LocalWind(x_i, y_i, h_i, wd, ws, time, wd_bin_size)
        return self._local_wind(lw, ws_bins)

    @abstractmethod
    def _local_wind(self, localWind, ws_bins=None):
        """Local free flow wind conditions

        Parameters
        ----------
        localWind  : LocalWind
            xarray dataset containing coordinates x, y, h, wd, ws
        ws_bin : array_like or None, optional
            Wind speed bin edges

        Returns
        -------
        LocalWind xarray dataset containing:
            WD : DataArray
                local free flow wind directions
            WS : DataArray
                local free flow wind speeds
            TI : DataArray
                local free flow turbulence intensity
            P : DataArray
                Probability/weight
        """

    def wt2wt_distances(self, WD_il):
        return self.distance(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):
        wd = np.atleast_1d(wd)
        if wd_bin_size is not None:
            return wd_bin_size
        elif len(wd) > 1 and len(np.unique(np.diff(wd))) == 1:
            return wd[1] - wd[0]
        else:
            return 360 / len(np.atleast_1d(wd))

    def ws_bins(self, WS, ws_bins=None):
        # TODO: delete function
        if not isinstance(WS, xr.DataArray):
            WS = xr.DataArray(WS, [('ws', np.atleast_1d(WS))])
        if not hasattr(ws_bins, '__len__') or len(ws_bins) != len(WS) + 1:

            if len(WS.shape) and WS.shape[-1] > 1:
                d = np.diff(WS) / 2
                ws_bins = 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
                ws_bins = WS.data + np.array([-ws_bins / 2, ws_bins / 2])
        else:
            ws_bins = np.asarray(ws_bins)
        return xr.Dataset({'ws_lower': (WS.dims, ws_bins[..., :-1]),
                           'ws_upper': (WS.dims, ws_bins[..., 1:])},
                          coords=WS.coords)

    def _sector(self, wd):
        sector = np.zeros(360, dtype=int)

        d_wd = (np.diff(np.r_[wd, wd[0]]) % 360) / 2
        assert np.all(d_wd == d_wd[0]), "Wind directions must be equidistant"
        lower = np.ceil(wd - d_wd).astype(int)
        upper = np.ceil(wd + d_wd).astype(int)
        for i, (lo, up) in enumerate(zip(lower, upper)):
            if lo < 0:
                sector[lo % 360 + 1:] = i
                lo = 0
            if up > 359:
                sector[:up % 360 + 1] = i
                up = 359
            sector[lo + 1:up + 1] = i
        return sector

    def plot_ws_distribution(self, x=0, y=0, h=70, wd=[0], ws=np.arange(0.05, 30.05, .1),
                             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)
        ws : array_like, optional
            Wind speeds to calculate for
        include_wd_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:
            ax = plt

        lbl = "Wind direction: %d deg"
        if include_wd_distribution:

            lw = self.local_wind(x_i=x, y_i=y, h_i=h, wd=np.arange(360), ws=ws, wd_bin_size=1)
            lw.coords['sector'] = ('wd', self._sector(wd))
            P = lw.P.groupby('sector').sum()
            v = 360 / len(wd) / 2
            lbl += r"$\pm$%s deg" % ((int(v), v)[(v % 2) != 0])
        else:
            lw = self.local_wind(x_i=x, y_i=y, h_i=h, wd=wd, ws=ws, wd_bin_size=1)
            P = lw.P
            if 'ws' not in P.dims:
                P = P.broadcast_like(lw.WS).T
            P = P / P.sum('ws')  # exclude wd probability
        if 'i' in P.dims:
            P = P.squeeze('i')
        for wd, p in zip(wd, P):
            ax.plot(ws, p * 10, label=lbl % wd)
            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/wd_bin_size and 10-20 m/wd_bin_size
        ax : pyplot or matplotlib axes object, default None
        """
        if ax is None:
            ax = plt
        assert 360 % n_wd == 0

        wd_bin_size = 360 // n_wd
        wd = np.arange(0, 360, wd_bin_size)
        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)

        if ws_bins is None:
            if any(['ws' in v.dims for v in self.ds.data_vars.values()]):
                lw = self.local_wind(x_i=x, y_i=y, h_i=h, wd=np.arange(360), wd_bin_size=1)
                lw['P'] = lw.P.sum('ws')
            else:
                lw = self.local_wind(x_i=x, y_i=y, h_i=h, wd=np.arange(360),
                                     ws=[100], ws_bins=[0, 200], wd_bin_size=1)
        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)
            lw = self.local_wind(x_i=x, y_i=y, h_i=h, wd=np.arange(360), ws=ws, wd_bin_size=1)

        lw.coords['sector'] = ('wd', self._sector(wd))
        p = lw.P.groupby('sector').sum()
        if 'i' in p.dims:
            p = p.squeeze('i')

        if ws_bins is None or 'ws' not in p.dims:
            if 'ws' in p.dims:
                p = p.squeeze('ws')
            ax.bar(theta, p.data, width=np.deg2rad(wd_bin_size), bottom=0.0)
        else:
            p = p.T
            start_p = np.vstack([np.zeros_like(p[:1]), p.cumsum('ws')[:-1]])
            for ws1, ws2, p_ws0, p_ws in zip(lw.ws_lower.data, lw.ws_upper.data, start_p, p):
                ax.bar(theta, p_ws, width=np.deg2rad(wd_bin_size), bottom=p_ws0,
                       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.T


from py_wake.site import xrsite  # @NoMove # nopep8
UniformSite = xrsite.UniformSite
UniformWeibullSite = xrsite.UniformWeibullSite


def get_sector_xr(v, name):
    if isinstance(v, (int, float)):
        return xr.DataArray(v, coords=[], name=name)
    v = np.r_[v, np.atleast_1d(v)[0]]
    return xr.DataArray(v, coords=[('wd', np.linspace(0, 360, len(v)))], name=name)


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, h_i=h_ref, wd=wdir_lst, ws=wsp_lst)
        print(local_wind.WS_ilk.shape)

        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()