add LogShear + test and test of custom shear function

py_wake/site/_site.py

@@ -117,14 +117,6 @@ class Site(ABC):
             ws = np.atleast_1d(ws)
         return wd, ws
 
-    def wref(self, wd, ws, ws_bins=None):
-        wd, ws = self.get_defaults(wd, ws)
-        WS = xr.DataArray(ws, [('ws', ws)])
-        ds = self.ws_bins(WS, ws_bins)
-        ds['WS'] = WS
-        ds['WD'] = xr.DataArray(wd, [('wd', wd)])
-        return ds
-
     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
 

py_wake/site/shear.py

@@ -23,7 +23,7 @@ class Shear(ABC):
         """
 
 
-class PowerShear():
+class PowerShear(Shear):
     def __init__(self, h_ref, alpha, interp_method='nearest'):
         self.h_ref = h_ref
         from py_wake.site._site import get_sector_xr
@@ -37,6 +37,20 @@ class PowerShear():
         return (h / self.h_ref) ** alpha * WS
 
 
+class LogShear(Shear):
+    def __init__(self, h_ref, z0, interp_method='nearest'):
+        self.h_ref = h_ref
+        from py_wake.site._site import get_sector_xr
+        self.z0 = get_sector_xr(z0, "Roughness length")
+        self.interp_method = interp_method
+
+    def __call__(self, WS, WD, h):
+        z0 = self.z0.interp_all(WD, method=self.interp_method)
+        if z0.shape == ():
+            z0 = z0.data
+        return np.log(h / z0) / np.log(self.h_ref / z0) * WS
+
+
 # ======================================================================================================================
 # Potentially the code below can be used to implement power/log shear interpolation between grid layers
 # ======================================================================================================================

py_wake/tests/test_sites/test_shear_models.py

-from py_wake.site.shear import PowerShear
+from py_wake.site.shear import PowerShear, LogShear
 import numpy as np
 from py_wake.tests import npt
 from py_wake.site._site import UniformSite
@@ -7,22 +7,65 @@ import matplotlib.pyplot as plt
 
 
 def test_power_shear():
-    shear = PowerShear(70, alpha=[.1, .2])
     h_lst = np.arange(10, 100, 10)
-    site = UniformSite([1], .1)
-    wref = site.wref([0, 180], [10])
-    h = xr.DataArray(np.arange(10, 100, 10), [('i', np.arange(9))])
-    u = shear(wref.WS, wref.WD, h)
+    site = UniformSite([1], .1, shear=PowerShear(70, alpha=[.1, .2]))
+    WS = site.local_wind(x_i=h_lst * 0, y_i=h_lst * 0, h_i=h_lst, wd=[0, 180], ws=[10, 12]).WS
 
     if 0:
-        plt.plot(u.sel(wd=0), h, label='alpha=0.1')
+        plt.plot(WS.sel(wd=0, ws=10), WS.h, label='alpha=0.1')
         plt.plot((h_lst / 70)**0.1 * 10, h_lst, ':')
-        plt.plot(u.sel(wd=180), h, label='alpha=0.2')
-        plt.plot((h_lst / 70)**0.2 * 10, h_lst, ':')
+        plt.plot(WS.sel(wd=180, ws=12), WS.h, label='alpha=0.2')
+        plt.plot((h_lst / 70)**0.2 * 12, h_lst, ':')
         plt.legend()
         plt.show()
-    npt.assert_array_almost_equal(u.sel(wd=0, ws=10), [8.23, 8.82, 9.19, 9.46, 9.67, 9.85, 10., 10.13, 10.25], 2)
-    npt.assert_array_almost_equal(u.sel(wd=180, ws=10), [6.78, 7.78, 8.44, 8.94, 9.35, 9.7, 10., 10.27, 10.52], 2)
+    npt.assert_array_equal(WS.sel(wd=0, ws=10), (h_lst / 70)**0.1 * 10)
+    npt.assert_array_equal(WS.sel(wd=180, ws=12), (h_lst / 70)**0.2 * 12)
+
+
+def test_log_shear():
+
+    h_lst = np.arange(10, 100, 10)
+    site = UniformSite([1], .1, shear=LogShear(70, z0=[.02, 2]))
+    WS = site.local_wind(x_i=h_lst * 0, y_i=h_lst * 0, h_i=h_lst, wd=[0, 180], ws=[10, 12]).WS
+
+    if 0:
+        plt.plot(WS.sel(wd=0, ws=10), WS.h, label='z0=0.02')
+        plt.plot(np.log(h_lst / 0.02) / np.log(70 / 0.02) * 10, h_lst, ':')
+        plt.plot(WS.sel(wd=180, ws=12), WS.h, label='z0=2')
+        plt.plot(np.log(h_lst / 2) / np.log(70 / 2) * 12, h_lst, ':')
+        plt.legend()
+        plt.show()
+    npt.assert_array_equal(WS.sel(wd=0, ws=10), np.log(h_lst / 0.02) / np.log(70 / 0.02) * 10)
+    npt.assert_array_equal(WS.sel(wd=180, ws=12), np.log(h_lst / 2) / np.log(70 / 2) * 12)
+
+
+def test_log_shear_constant_z0():
+    h_lst = np.arange(10, 100, 10)
+    site = UniformSite([1], .1, shear=LogShear(70, z0=.02))
+    WS = site.local_wind(x_i=h_lst * 0, y_i=h_lst * 0, h_i=h_lst, wd=[0, 180], ws=[10, 12]).WS
+
+    if 0:
+        plt.plot(WS.sel(ws=10), WS.h, label='z0=0.02')
+        plt.plot(np.log(h_lst / 0.02) / np.log(70 / 0.02) * 10, h_lst, ':')
+        plt.legend()
+        plt.show()
+    npt.assert_array_equal(WS.sel(ws=10), np.log(h_lst / 0.02) / np.log(70 / 0.02) * 10)
+
+
+def test_custom_shear():
+    def my_shear(WS, WD, h):
+        return WS * (0.02 * (h - 70) + 1) * (WD * 0 + 1)
+    h_lst = np.arange(10, 100, 10)
+
+    site = UniformSite([1], .1, shear=my_shear)
+    WS = site.local_wind(x_i=h_lst * 0, y_i=h_lst * 0, h_i=h_lst, wd=[0, 180], ws=[10, 12]).WS
+
+    if 0:
+        plt.plot(WS.sel(wd=0, ws=10), WS.h, label='2z-2')
+        plt.plot((h_lst - 70) * 0.2 + 10, h_lst, ':')
+        plt.legend()
+        plt.show()
+    npt.assert_array_almost_equal(WS.sel(wd=0, ws=10), (h_lst - 70) * 0.2 + 10)
 
 
 if __name__ == '__main__':