import numpy as np from py_wake.site.shear import PowerShear from py_wake.site.xrsite import XRSite import xarray as xr from py_wake.tests import npt import pytest import matplotlib.pyplot as plt from numpy import newaxis as na from py_wake.tests.test_files import tfp from py_wake.examples.data.hornsrev1 import Hornsrev1Site, V80, wt9_x, wt9_y from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site from py_wake.wind_turbines import WindTurbines from py_wake.deficit_models.gaussian import BastankhahGaussian, BastankhahGaussianDeficit from py_wake.wind_farm_models.engineering_models import PropagateDownwind from py_wake.superposition_models import LinearSum from py_wake.tests.check_speed import timeit from py_wake.site._site import LocalWind 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 @pytest.fixture def uniform_site(): ti = 0.1 ds = xr.Dataset( data_vars={'WS': 10, 'P': ('wd', f), 'TI': ti}, coords={'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def uniform_weibull_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def complex_fixed_pos_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Speedup': ('i', np.arange(.8, 1.3, .1)), 'Turning': ('i', np.arange(-2, 3)), 'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'i': np.arange(5), 'wd': np.linspace(0, 360, len(f), endpoint=False)}) x_i = np.arange(5) return XRSite(ds, initial_position=np.array([x_i, x_i + 1]).T, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def complex_grid_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Speedup': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))), 'Turning': (['x', 'y'], np.arange(-2, 4, 1).reshape((3, 2))), 'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'x': [0, 5, 10], 'y': [0, 5], 'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2), interp_method='linear') def test_uniform_local_wind(uniform_site): site = uniform_site x_i = y_i = np.arange(5) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, 10) npt.assert_array_equal(lw.WD, wdir_lst) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_equal(lw.P, np.array([0.035972, 0.070002, 0.086432, 0.147379]) / 30) npt.assert_array_equal(site.elevation(x_i, y_i), 0) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100) npt.assert_array_equal(lw.WD_ilk.shape, (1, 360, 1)) z = np.arange(1, 200) zero = [0] * len(z) ws100_2 = site.local_wind(x_i=zero, y_i=zero, h_i=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0] site.shear = PowerShear(70, alpha=.3) ws70_3 = site.local_wind(x_i=zero, y_i=zero, h_i=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0] if 0: plt.plot(ws100_2, z) plt.plot(ws70_3, z) plt.show() npt.assert_array_equal(10 * (z / 100)**.2, ws100_2) npt.assert_array_equal(10 * (z / 70)**.3, ws70_3) def test_uniform_weibull_local_wind(uniform_weibull_site): site = uniform_weibull_site x_i = y_i = np.arange(5) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, [3, 4, 5]) npt.assert_array_equal(lw.WD, wdir_lst) npt.assert_array_equal(lw.TI, 0.1) # ref_site = UniformWeibullSite(p_wd=f, a=A, k=k, ti=0.1, shear=site.shear) # lw_ref = ref_site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) # print(lw_ref.P) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_complex_fixed_pos_local_wind(complex_fixed_pos_site): site = complex_fixed_pos_site x_i, y_i = site.initial_position.T npt.assert_array_equal(x_i, np.arange(5)) npt.assert_array_equal(y_i, np.arange(5) + 1) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, [3, 4, 5] * np.arange(.8, 1.3, .1)[:, na]) npt.assert_array_equal(lw.WD, (wdir_lst + np.arange(-2, 3)[:, na]) % 360) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_complex_grid_local_wind(complex_grid_site): site = complex_grid_site y = np.arange(5) x = y * 2 X, Y = np.meshgrid(x, y) x_i, y_i = X.flatten(), Y.flatten() wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) if 0: c = plt.contourf(X, Y, lw.WS.sel(ws=5).data.reshape(X.shape)) plt.colorbar(c) plt.figure() c = plt.contourf(X, Y, lw.WD.sel(wd=90).data.reshape(X.shape)) plt.colorbar(c) plt.show() lw = site.local_wind(x_i=[2.5, 7.5], y_i=[2.5, 2.5], h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_almost_equal(lw.WS, [3, 4, 5] * np.array([0.95, 1.15])[:, na]) npt.assert_array_equal(lw.WD, (wdir_lst + np.array([-.5, 1.5])[:, na]) % 360) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_wd_independent_site(): ti = 0.1 ds = xr.Dataset( data_vars={ 'WS': 10, 'Sector_frequency': 1, 'Weibull_A': 4, 'Weibull_k': 2, 'TI': ti}, coords={}) site = XRSite(ds, shear=None) npt.assert_equal(site.ds.sector_width, 360) def test_load_save(complex_grid_site): complex_grid_site.save(tfp + "tmp.nc") site = XRSite.load(tfp + "tmp.nc", interp_method='linear') test_complex_grid_local_wind(site) def test_elevation(): ti = 0.1 ds = xr.Dataset( data_vars={'Elevation': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))), 'P': 1, 'TI': ti}, coords={'x': [0, 5, 10], 'y': [0, 5]}) site = XRSite(ds) npt.assert_array_almost_equal(site.elevation([2.5, 7.5], [2.5, 2.5]), [0.95, 1.15]) def test_cyclic_interpolation(uniform_site): site = uniform_site lw = site.local_wind([0], [0], 100, wd=np.arange(360), ws=10) if 0: plt.plot(np.linspace(0, 360, len(f) + 1), np.r_[f, f[0]], '.') plt.plot(lw.wd, lw.P) plt.show() npt.assert_array_almost_equal(lw.P[::30], np.array(f) / 30) def test_from_flow_box_2wt(): site = Hornsrev1Site() windTurbines = V80() # simulate current and neighbour wt wfm = BastankhahGaussian(site, windTurbines) wd = np.arange(30) sim_res = wfm([0, 0], [0, 500], wd=wd) ref_aep = sim_res.aep().sel(wt=0) wt_x, wt_y = [0], [0] neighbour_x, neighbour_y = [0], [500] # make site with effects of neighbour wt sim_res = wfm(neighbour_x, neighbour_y, wd=wd) e = 100 box = sim_res.flow_box(x=np.linspace(min(wt_x) - e, max(wt_x) + e, 21), y=np.linspace(min(wt_y) - e, max(wt_y) + e, 21), h=windTurbines.hub_height(windTurbines.types())) site = XRSite.from_flow_box(box) # Simujlate current wt and compare aep wfm = BastankhahGaussian(site, windTurbines) sim_res = wfm(wt_x, wt_y, wd=wd) aep = sim_res.aep() if 0: site.ds.WS.sel(ws=10, wd=3).plot() windTurbines.plot(wt_x, wt_y) windTurbines.plot(neighbour_x, neighbour_y) plt.show() npt.assert_array_almost_equal(ref_aep, aep.sel(wt=0)) def test_neighbour_farm_speed(): # import and setup site and windTurbines site = IEA37Site(16) # setup current, neighbour and all positions wt_x, wt_y = site.initial_position.T neighbour_x, neighbour_y = wt_x - 4000, wt_y all_x, all_y = np.r_[wt_x, neighbour_x], np.r_[wt_y, neighbour_y] windTurbines = WindTurbines.from_WindTurbines([IEA37_WindTurbines(), IEA37_WindTurbines()]) windTurbines._names = ["Current wind farm", "Neighbour wind farm"] types = [0] * len(wt_x) + [1] * len(neighbour_x) wf_model = PropagateDownwind(site, windTurbines, wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True), superpositionModel=LinearSum()) # Consider wd=270 +/- 60 deg only wd_lst = np.arange(210, 331) sim_res, t = timeit(wf_model, verbose=True)(all_x, all_y, type=types, ws=9.8, wd=wd_lst) if 1: ext = 100 flow_box = wf_model(neighbour_x, neighbour_y, wd=wd_lst).flow_box( x=np.linspace(min(wt_x) - ext, max(wt_x) + ext, 51), y=np.linspace(min(wt_y) - ext, max(wt_y) + ext, 51), h=[100, 110, 120]) wake_site = XRSite.from_flow_box(flow_box) wake_site.save('tmp.nc') else: wake_site = XRSite.load('tmp.nc') wf_model_wake_site = PropagateDownwind(wake_site, windTurbines, wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True), superpositionModel=LinearSum()) if 0: from line_profiler import LineProfiler lp = LineProfiler() lp.timer_unit = 1e-6 lp.add_function(XRSite.interp) lp_wrapper = lp(wf_model_wake_site.__call__) sim_res_wake_site = lp_wrapper(wt_x, wt_y, ws=9.8, wd=wd_lst) lp.print_stats() else: sim_res_wake_site, _ = timeit(wf_model_wake_site, verbose=True)(wt_x, wt_y, ws=9.8, wd=wd_lst) npt.assert_allclose(sim_res.aep().sel(wt=np.arange(len(wt_x))).sum(), sim_res_wake_site.aep().sum(), rtol=0.0005) npt.assert_array_almost_equal(sim_res.aep().sel(wt=np.arange(len(wt_x))), sim_res_wake_site.aep(), 2) @pytest.mark.parametrize('h,wd,ws,h_i,wd_l,ws_k', [ ([100, 110, 120], range(360), [9.8], 110, range(360), [9.8]), ([100, 110, 120], range(5, 25), [9.8], 110, range(10, 20), [9.8]), ([100, 110, 120], range(5, 25), [9, 10, 11], 110, range(10, 20), [10]), ([100, 110, 120], range(5, 25), [8, 10, 11], 110, range(10, 20), [9.8]), ]) def test_interp(h, wd, ws, h_i, wd_l, ws_k): ds = xr.Dataset({ 'TI': (['x', 'y', 'h', 'wd', 'ws'], np.random.rand(10, 20, len(h), len(wd), len(ws))), 'P': 1 }, coords={'x': np.linspace(0, 100, 10), 'y': np.linspace(200, 300, 20), 'h': h, 'wd': wd, 'ws': ws} ) site = XRSite(ds) lw = LocalWind(x_i=[25, 50], y_i=[225, 250], h_i=h_i, wd=wd_l, ws=ws_k, wd_bin_size=1) ip1 = site.interp(site.ds.TI, lw.coords) ip2 = ds.TI.sel_interp_all(lw.coords) npt.assert_array_equal(ip1.shape, ip2.shape) npt.assert_array_almost_equal(ip1.data, ip2.data) def test_interp2(): wd = np.arange(5) ws = np.arange(10) ds = xr.Dataset({ 'TI': (['i', 'wd', 'ws'], np.random.rand(10, len(wd), len(ws))), 'P': 1 }, coords={'i': np.arange(10), 'wd': wd, 'ws': ws} ) site = XRSite(ds) with pytest.raises(ValueError, match=r"Number of points, i\(=10\), in site data variable, TI, must match "): lw = LocalWind(x_i=[25, 50], y_i=[225, 250], h_i=110, wd=wd, ws=ws, wd_bin_size=1) site.interp(site.ds.TI, lw.coords) x = y = np.arange(10) lw = LocalWind(x_i=x, y_i=y, h_i=110, wd=wd, ws=ws, wd_bin_size=1) ip1 = site.interp(site.ds.TI, lw.coords) ip2 = ds.TI.sel_interp_all(lw.coords) npt.assert_array_equal(ip1.shape, ip2.shape) npt.assert_array_almost_equal(ip1.data, ip2.data)