test_wasp_grid_site.py

import numpy as np
from numpy import newaxis as na
from py_wake.tests import npt
import pytest
from py_wake.examples.data.ParqueFicticio import ParqueFicticio_path, ParqueFicticioSite
from py_wake.site.wasp_grid_site import WaspGridSite
import os
import time
from py_wake.tests.test_files.wasp_grid_site import one_layer
from py_wake.site.distance import TerrainFollowingDistance, StraightDistance, TerrainFollowingDistance2
import math
from py_wake import NOJ
from py_wake.wind_turbines import OneTypeWindTurbines
from py_wake.flow_map import HorizontalGrid
from py_wake.tests.check_speed import timeit
import matplotlib.pyplot as plt


@pytest.fixture
def site():
    return ParqueFicticioSite()


@pytest.fixture
def site2():
    site = ParqueFicticioSite(distance=TerrainFollowingDistance2())
    x, y = site.initial_position.T
    return site, x, y


def test_WaspGridSiteDistanceClass(site):
    wgs = WaspGridSite(site._ds, distance=TerrainFollowingDistance(distance_resolution=2000))
    assert wgs.distance.distance_resolution == 2000
    assert wgs.distance.__call__.__func__ == TerrainFollowingDistance().__call__.__func__
    wgs = WaspGridSite(site._ds, distance=StraightDistance())
    assert wgs.distance.__call__.__func__ == StraightDistance().__call__.__func__


def test_local_wind(site):
    x_i, y_i = site.initial_position.T
    h_i = x_i * 0 + 70
    wdir_lst = np.arange(0, 360, 90)
    wsp_lst = np.arange(3, 6)
    WS_ilk = site.local_wind(x_i=x_i, y_i=y_i, h_i=h_i, wd=wdir_lst, ws=wsp_lst).WS_ilk
    npt.assert_array_equal(WS_ilk.shape, (8, 4, 3))

    WS_ilk = site.local_wind(x_i=x_i, y_i=y_i, h_i=h_i).WS_ilk
    npt.assert_array_equal(WS_ilk.shape, (8, 360, 23))

    # check probability local_wind()[-1]
    npt.assert_almost_equal(site.local_wind(x_i=x_i, y_i=y_i, h_i=h_i, wd=[0], ws=[10], wd_bin_size=1).P_ilk,
                            site.local_wind(x_i=x_i, y_i=y_i, h_i=h_i, wd=[0], ws=[10], wd_bin_size=2).P_ilk / 2, 6)


def test_shear(site):
    npt.assert_array_almost_equal(site._ds['spd'].sel(x=262878, y=6504714, sec=1), [.6240589, .8932919])
    x = [262878.0001] * 3
    y = [6504714.0001] * 3
    z = [30, 115, 200]
    ws = site.local_wind(x_i=x, y_i=y, h_i=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0]

    if 0:
        plt.plot(ws, z, '.-')
        plt.show()
    plt.close()

    # linear interpolation
    npt.assert_array_almost_equal(ws, [6.240589, np.mean([6.240589, 8.932919]), 8.932919])


def test_wasp_resources_grid_point(site):
    #     x = np.array([l.split() for l in """0.6010665    -10.02692    32.71442    -6.746912
    # 0.5007213    -4.591617    37.10247    -11.0699
    # 0.3104101    -1.821247    59.18301    -12.56743
    # 0.4674515    16.14293    44.84665    -9.693183
    # 0.8710347    5.291974    26.01634    -6.154611
    # 0.9998786    -2.777032    15.72486    1.029988
    # 0.9079611    -7.882853    16.33216    6.42329
    # 0.759553    -5.082487    17.23354    10.18187
    # 0.7221162    4.606324    17.96417    11.45838
    # 0.8088576    8.196074    16.16308    9.277925
    # 0.8800673    3.932325    14.82337    5.380589
    # 0.8726974    -3.199536    19.99724    -1.433086""".split("\n")], dtype=np.float)
    #     for x_ in x.T:
    #         print(list(x_))
    x = [262978]
    y = [6504814]
    npt.assert_almost_equal(site.elevation(x, y), 227.8, 1)

    # Data from WAsP:
    # - add turbine (262878,6504814,30)
    # - Turbine (right click) - reports - Turbine Site Report full precision
    wasp_A = [2.197305, 1.664085, 1.353185, 2.651781, 5.28438, 5.038289,
              4.174325, 4.604496, 5.043066, 6.108261, 6.082033, 3.659798]
    wasp_k = [1.771484, 2.103516, 2.642578, 2.400391, 2.357422, 2.306641,
              2.232422, 2.357422, 2.400391, 2.177734, 1.845703, 1.513672]
    wasp_f = [5.188083, 2.509297, 2.869334, 4.966141, 13.16969, 9.514355,
              4.80275, 6.038354, 9.828702, 14.44174, 16.60567, 10.0659]
    wasp_spd = [0.6010665, 0.5007213, 0.3104101, 0.4674515, 0.8710347, 0.9998786,
                0.9079611, 0.759553, 0.7221162, 0.8088576, 0.8800673, 0.8726974]
    wasp_trn = [-10.02692, -4.591617, -1.821247, 16.14293, 5.291974, -
                2.777032, -7.882853, -5.082487, 4.606324, 8.196074, 3.932325, -3.199536]
    wasp_inc = [-6.746912, -11.0699, -12.56743, -9.693183, -6.154611,
                1.029988, 6.42329, 10.18187, 11.45838, 9.277925, 5.380589, -1.433086]
    wasp_ti = [32.71442, 37.10247, 59.18301, 44.84665, 26.01634, 15.72486,
               16.33216, 17.23354, 17.96417, 16.16308, 14.82337, 19.99724]
    rho = 1.179558

    wasp_u_mean = [1.955629, 1.473854, 1.202513, 2.350761, 4.683075,
                   4.463644, 3.697135, 4.080554, 4.470596, 5.409509, 5.402648, 3.300305]
    wasp_p_air = [9.615095, 3.434769, 1.556282, 12.45899, 99.90289,
                  88.03519, 51.41135, 66.09097, 85.69466, 164.5592, 193.3779, 56.86945]
#     wasp_aep = np.array([3725293.0, 33722.71, 0.3093564, 3577990.0, 302099600.0, 188784100.0,
#                          48915640.0, 84636210.0, 189009800.0, 549195100.0, 691258600.0, 120013000.0]) / 1000
    wasp_aep_no_density_correction = np.array([3937022.0, 36046.93, 0.33592, 3796496.0, 314595600.0,
                                               196765700.0, 51195440.0, 88451200.0, 197132700.0, 568584400.0, 712938400.0, 124804600.0]) / 1000
#     wasp_aep_total = 2.181249024
    wasp_aep_no_density_correction_total = 2.26224
    wt_u = np.array([3.99, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25])
    wt_p = np.array([0, 55., 185., 369., 619., 941., 1326., 1741., 2133., 2436., 2617., 2702., 2734.,
                     2744., 2747., 2748., 2748., 2750., 2750., 2750., 2750., 2750., 2750.])
    wt_ct = np.array([0, 0.871, 0.853, 0.841, 0.841, 0.833, 0.797, 0.743, 0.635, 0.543, 0.424,
                      0.324, 0.258, 0.21, 0.175, 0.147, 0.126, 0.109, 0.095, 0.083, 0.074, 0.065, 0.059])
    wt = OneTypeWindTurbines.from_tabular(name="NEG-Micon 2750/92 (2750 kW)", diameter=92, hub_height=70,
                                          ws=wt_u, ct=wt_ct, power=wt_p, power_unit='kw')

    A_lst, k_lst, f_lst, spd_lst, orog_trn_lst, flow_inc_lst, tke_lst = [site.interp_funcs[n](
        (x, y, 30, range(0, 360, 30))) for n in ['A', 'k', 'f', 'spd', 'orog_trn', 'flow_inc', 'tke']]
    f_lst = f_lst * 360 / 12
    pdf_lst = [lambda x, A=A, k=k: k / A * (x / A)**(k - 1) * np.exp(-(x / A)**k) * (x[1] - x[0])
               for A, k in zip(A_lst, k_lst)]
#     cdf_lst = [lambda x, A=A, k=k: 1 - np.exp(-(x / A) ** k) for A, k in zip(A_lst, k_lst)]
    dx = .1
    ws = np.arange(dx / 2, 35, dx)

    # compare to wasp data
    npt.assert_array_equal(A_lst, wasp_A)
    npt.assert_array_equal(k_lst, wasp_k)
    npt.assert_array_almost_equal(f_lst, np.array(wasp_f) / 100)
    npt.assert_array_almost_equal(spd_lst, wasp_spd)
    npt.assert_array_almost_equal(orog_trn_lst, wasp_trn)
    npt.assert_array_almost_equal(flow_inc_lst, wasp_inc)
    npt.assert_array_almost_equal(tke_lst, np.array(wasp_ti) / 100)

    # compare pdf, u_mean and aep to wasp
    lw = site.local_wind(x, np.array(y) + 1e-6, 30, wd=np.arange(0, 360, 30), ws=ws)
    P = lw.P_ilk / lw.P_ilk.sum(2)[:, :, na]  # only wind speed probablity (not wdir)

    # pdf
    for l in range(12):
        npt.assert_array_almost_equal(np.interp(ws, lw.WS_ilk[0, l], np.cumsum(P[0, l])),
                                      np.cumsum(pdf_lst[l](ws)), 1)

    # u_mean
    npt.assert_almost_equal([A * math.gamma(1 + 1 / k) for A, k in zip(A_lst, k_lst)], wasp_u_mean, 5)
    npt.assert_almost_equal([(pdf(ws) * ws).sum() for pdf in pdf_lst], wasp_u_mean, 5)
    npt.assert_almost_equal((P * lw.WS_ilk).sum((0, 2)), wasp_u_mean, 5)

    # air power
    p_air = [(pdf(ws) * 1 / 2 * rho * ws**3).sum() for pdf in pdf_lst]
    npt.assert_array_almost_equal(p_air, wasp_p_air, 3)
    npt.assert_array_almost_equal((P * 1 / 2 * rho * lw.WS_ilk**3).sum((0, 2)), wasp_p_air, 2)

    # AEP
    AEP_ilk = NOJ(site, wt)(x, y, h=30, wd=np.arange(0, 360, 30), ws=ws).aep_ilk(
        with_wake_loss=False, normalize_probabilities=True)
    if 0:
        plt.plot(wasp_aep_no_density_correction / 1000, '.-', label='WAsP')
        plt.plot(AEP_ilk.sum((0, 2)) * 1e3, label='PyWake')
        plt.xlabel('Sector')
        plt.ylabel('AEP [MWh]')
        plt.legend()
        plt.show()
    npt.assert_array_less(np.abs(wasp_aep_no_density_correction - AEP_ilk.sum((0, 2)) * 1e6), 300)
    npt.assert_almost_equal(AEP_ilk.sum(), wasp_aep_no_density_correction_total, 3)


@pytest.mark.parametrize('site,dw_ref', [
    (ParqueFicticioSite(distance=TerrainFollowingDistance2()),
     [0., 207.3842238, 484.3998264, 726.7130743, 1039.148129, 1263.1335982, 1490.3841602, 1840.6508086]),
    (ParqueFicticioSite(distance=TerrainFollowingDistance()),
     [0, 209.803579, 480.8335365, 715.6003233, 1026.9476322, 1249.5510034, 1475.1467251, 1824.1317343]),
    (ParqueFicticioSite(distance=StraightDistance()),
     [-0, 207, 477, 710, 1016, 1236, 1456, 1799])])
def test_distances(site, dw_ref):
    x, y = site.initial_position.T
    dw_ijl, cw_ijl, dh_ijl, _ = site.distances(src_x_i=x, src_y_i=y, src_h_i=np.array([70]),
                                               dst_x_j=x, dst_y_j=y, dst_h_j=np.array([70]),
                                               wd_il=np.array([[0]]))
    npt.assert_almost_equal(dw_ijl[0, :, 0], dw_ref)

    cw_ref = [236.1, 0., -131.1, -167.8, -204.5, -131.1, -131.1, -45.4]
    npt.assert_almost_equal(cw_ijl[:, 1, 0], cw_ref)
    npt.assert_almost_equal(dh_ijl, np.zeros_like(dh_ijl))


def test_distances_different_points(site2):
    site, x, y = site2
    with pytest.raises(NotImplementedError):
        site.distances(src_x_i=x, src_y_i=y, src_h_i=np.array([70]),
                       dst_x_j=x[1:], dst_y_j=y[1:], dst_h_j=np.array([70]),
                       wd_il=np.array([[0]]))


# def test_distances_wd_shape():
#     site = ParqueFicticioSite(distance=TerrainFollowingDistance2())
#     x, y = site.initial_position.T
#     dw_ijl, cw_ijl, dh_ijl, dwo = site.distances(src_x_i=x, src_y_i=y, src_h_i=np.array([70]),
#                                                  dst_x_j=x, dst_y_j=y, dst_h_j=np.array([70]),
#                                                  wd_il=np.ones((len(x), 1)) * 180)
#     npt.assert_almost_equal(dw_ijl[0, :, 0], np.array([0., -207., -477., -710., -1016., -1236., -1456., -1799.]))
#     npt.assert_almost_equal(cw_ijl[:, 1, 0], np.array([-236.1, 0., 131.1, 167.8, 204.5, 131.1, 131.1, 45.4]))
#     npt.assert_almost_equal(dh_ijl, np.zeros_like(dh_ijl))


def test_speed_up_using_pickle():
    pkl_fn = ParqueFicticio_path + "ParqueFicticio.pkl"
    if os.path.exists(pkl_fn):
        os.remove(pkl_fn)
    start = time.time()
    site = WaspGridSite.from_wasp_grd(ParqueFicticio_path, speedup_using_pickle=False)
    time_wo_pkl = time.time() - start
    site = WaspGridSite.from_wasp_grd(ParqueFicticio_path, speedup_using_pickle=True)
    assert os.path.exists(pkl_fn)
    start = time.time()
    site = WaspGridSite.from_wasp_grd(ParqueFicticio_path, speedup_using_pickle=True)
    time_w_pkl = time.time() - start
    npt.assert_array_less(time_w_pkl * 10, time_wo_pkl)


def test_interp_funcs_initialization_missing_key(site):
    site = ParqueFicticioSite(distance=TerrainFollowingDistance2())
    site.interp_funcs_initialization(['missing'])


def test_one_layer():
    site = WaspGridSite.from_wasp_grd(os.path.dirname(one_layer.__file__) + "/", speedup_using_pickle=False)


def test_missing_path():
    with pytest.raises(NotImplementedError):
        WaspGridSite.from_wasp_grd("missing_path/", speedup_using_pickle=True)

    with pytest.raises(Exception, match='Path was not a directory'):
        WaspGridSite.from_wasp_grd("missing_path/", speedup_using_pickle=False)


def test_elevation(site):
    x_i, y_i = site.initial_position.T
    npt.assert_array_less(np.abs(site.elevation(x_i=x_i, y_i=y_i) -
                                 [519.4, 567.7, 583.6, 600, 574.8, 559.9, 517.7, 474.5]  # ref from wasp
                                 ), 5)


def test_plot_map(site):
    with pytest.raises(AttributeError, match="missing not found in dataset. Available data variables are:\nflow_inc,"):
        site.plot_map('missing')

    with pytest.raises(AttributeError, match=r"Sector None not found. Available sectors are: \[ 1"):
        site.plot_map('ws_mean')
    with pytest.raises(AttributeError, match="Height missing for 'ws_mean'"):
        site.plot_map('ws_mean', sector=1)

    site.plot_map('elev')
    plt.figure()
    site.plot_map('ws_mean', 80, sector=1)
    if 0:
        plt.show()
    plt.close()


def test_elevation_outside_map(site):
    site.plot_map('elev')
    x = np.linspace(262500, 265500, 500)
    y = x * 0 + 6505450
    plt.plot(x, y, '--', label='Terrain line')
    plt.plot(x, y + site.elevation(x, y), label='Elevation')
    npt.assert_array_equal(np.round(site.elevation(x, y)[::50]),
                           [np.nan, np.nan, 303, 390, 491, 566, 486, 524, np.nan, np.nan])
    if 0:
        plt.legend()
        plt.show()
    plt.close()


def test_plot_ws_distribution(site):
    x, y = site.initial_position[0]
    p1 = site.plot_ws_distribution(x=x, y=y, h=70, wd=[0, 90, 180, 270])
    p2 = site.plot_ws_distribution(x=x, y=y, h=70, wd=[0, 90, 180, 270], include_wd_distribution=True)
    if 0:
        print(np.round(p1[-1, ::30].data, 4).tolist())
        print(np.round(p2[-1, ::30].data, 4).tolist())
        plt.show()

    npt.assert_array_almost_equal(p1[-1, ::30], [0.0001, 0.0079, 0.011, 0.0082, 0.0039,
                                                 0.0013, 0.0003, 0.0, 0.0, 0.0], 4)
    npt.assert_array_almost_equal(p2[-1, ::30], [0.0001, 0.0036, 0.0047, 0.0033, 0.0014,
                                                 0.0004, 0.0001, 0.0, 0.0, 0.0], 4)
    plt.close()


def test_plot_wd_distribution(site):
    x, y = site.initial_position[0]
    p = site.plot_wd_distribution(x=x, y=y, h=70, n_wd=12, ax=plt)
    # print(np.round(p, 3).tolist())

    if 0:
        plt.show()
    plt.close()
    npt.assert_array_almost_equal(p, [0.052, 0.043, 0.058, 0.085, 0.089, 0.061,
                                      0.047, 0.083, 0.153, 0.152, 0.108, 0.068], 3)


def test_plot_wd_distribution_with_ws_levels(site):
    x, y = site.initial_position[0]
    p = site.plot_wd_distribution(x=x, y=y, n_wd=12, ws_bins=[0, 5, 10, 15, 20, 25])
    if 0:
        print(np.round(p, 3).data.tolist())
        plt.show()
    npt.assert_array_almost_equal(p, [[0.039, 0.013, 0.001, 0.0, 0.0],
                                      [0.034, 0.009, 0.0, 0.0, 0.0],
                                      [0.035, 0.022, 0.0, 0.0, 0.0],
                                      [0.034, 0.048, 0.003, 0.0, 0.0],
                                      [0.031, 0.052, 0.007, 0.0, 0.0],
                                      [0.03, 0.03, 0.002, 0.0, 0.0],
                                      [0.027, 0.019, 0.001, 0.0, 0.0],
                                      [0.034, 0.043, 0.006, 0.0, 0.0],
                                      [0.047, 0.082, 0.023, 0.001, 0.0],
                                      [0.048, 0.074, 0.026, 0.003, 0.0],
                                      [0.044, 0.046, 0.015, 0.002, 0.0],
                                      [0.041, 0.023, 0.003, 0.0, 0.0]], 3)
    plt.close()


def test_additional_input():
    site = ParqueFicticioSite()
    wgs = WaspGridSite(site._ds, distance=TerrainFollowingDistance(distance_resolution=2000))
    wgs.interp_funcs_initialization(['ws_mean'])
    x, y = site.initial_position.T
    h = 70 * np.ones_like(x)
    ws_mean, = wgs.interpolate(['ws_mean'], x, y, h)
    npt.assert_array_almost_equal(ws_mean[0, :50],
                                  np.array([4.77080802, 4.77216214, 4.77351626, 4.77487037, 4.77622449,
                                            4.77757861, 4.77893273, 4.78028685, 4.78164097, 4.78299508,
                                            4.7843492, 4.78570332, 4.78705744, 4.78841156, 4.78976567,
                                            4.79111979, 4.79247391, 4.79382803, 4.79518215, 4.79653626,
                                            4.79789038, 4.7992445, 4.80059862, 4.80195274, 4.80330685,
                                            4.80466097, 4.80601509, 4.80736921, 4.80872333, 4.81007744,
                                            4.81143156, 4.87114138, 4.9308512, 4.99056102, 5.05027084,
                                            5.10998066, 5.16969047, 5.22940029, 5.28911011, 5.34881993,
                                            5.40852975, 5.46823957, 5.52794939, 5.5876592, 5.64736902,
                                            5.70707884, 5.76678866, 5.82649848, 5.8862083, 5.94591812]))


def test_interpolation_speed():
    import xarray as xr
    da = xr.DataArray(np.sin(0.3 * np.arange(20).reshape(5, 4)),
                      [('x', np.arange(5)),
                       ('y', [0.1, 0.2, 0.3, 0.4])])
    x = xr.DataArray([0.5, 1.5, 2.5], dims='z')
    y = xr.DataArray([0.15, 0.25, 0.35], dims='z')
    da.interp(x=x, y=y)

    site = ParqueFicticioSite()
    x, y = site.initial_position.T
    X, Y, x_j, y_j, h_j = HorizontalGrid()(x, y, 70)
    wd = [270]  # site.default_wd
    ws = site.default_ws
    res1, t_lst = timeit(site.interp_funcs['A'])((x_j, y_j, h_j, x_j * 0 + 270))
    print(res1.shape)
    res2, t_lst = timeit(lambda x, y, z, sec:
                         site._ds.A.interp(x=xr.DataArray(x, dims='z'),
                                           y=xr.DataArray(y, dims='z'),
                                           z=xr.DataArray(z, dims='z'),
                                           sec=xr.DataArray(sec, dims='z')).data)(x_j, y_j, h_j, x_j * 0 + 10)
    npt.assert_array_almost_equal(res1, res2)
    if 0:
        c = plt.contourf(X, Y, res1.reshape(X.shape))
        plt.colorbar(c)
        plt.figure()
        c = plt.contourf(X, Y, res2.reshape(X.shape))
        plt.colorbar(c)
        plt.show()


if __name__ == '__main__':
    test_wasp_resources_grid_point(ParqueFicticioSite())