new site based on global wind atlas data

py_wake/site/xrsite.py

-import xarray as xr
-from py_wake.site.distance import StraightDistance
 import numpy as np
-from py_wake.utils.grid_interpolator import GridInterpolator, EqDistRegGrid2DInterpolator
+from numpy import newaxis as na
+import xarray as xr
+from requests import get
 from py_wake.site._site import Site
+from py_wake.site.distance import StraightDistance
 from py_wake.utils import weibull
+from py_wake.utils.grid_interpolator import GridInterpolator, EqDistRegGrid2DInterpolator
 
 
 class XRSite(Site):
@@ -309,3 +311,47 @@ class UniformWeibullSite(XRSite):
         if ti is not None:
             ds['TI'] = ti
         XRSite.__init__(self, ds, interp_method=interp_method, shear=shear)
+
+
+class GlobalWindAtlasSite(XRSite):
+    """Site with Global Wind Climate (GWC) from the Global Wind Atlas based on
+    lat and long which is interpolated at specific roughness and height.
+    """
+
+    def __init__(self, lat, long, height, roughness, ti=None, **kwargs):
+        """
+        Parameters
+        ----------
+        lat: float
+            Latitude of the location
+        long: float
+            Longitude of the location
+        height: float
+            Height of the location
+        roughness: float
+            roughness length at the location
+        """
+        self.gwc_ds = self._read_gwc(lat, long)
+        if ti is not None:
+            self.gwc_ds['TI'] = ti
+        XRSite.__init__(self, ds=self.gwc_ds.interp(height=height, roughness=roughness), **kwargs)
+
+    def _read_gwc(self, lat, long):
+        url_str = f'https://globalwindatlas.info/api/gwa/custom/Lib/?lat={lat}&long={long}'
+        lines = get(url_str).text.strip().split('\r\n')
+
+        # Read header information one line at a time
+        # desc = txt[0].strip()  # File Description
+        nrough, nhgt, nsec = map(int, lines[1].split())  # dimensions
+        roughnesses = np.array(lines[2].split(), dtype=float)  # Roughness classes
+        heights = np.array(lines[3].split(), dtype=float)  # heights
+        data = np.array([l.split() for l in lines[4:]], dtype=float).reshape((nrough, nhgt * 2 + 1, nsec))
+        freq = data[:, 0] / data[:, 0].sum(1)[:, na]
+        A = data[:, 1::2]
+        k = data[:, 2::2]
+        ds = xr.Dataset({'Weibull_A': (["roughness", "height", "wd"], A),
+                         'Weibull_k': (["roughness", "height", "wd"], k),
+                         "Sector_frequency": (["roughness", "wd"], freq)},
+                        coords={"height": heights, "roughness": roughnesses,
+                                "wd": np.linspace(0, 360, nsec, endpoint=False)})
+        return ds

py_wake/tests/test_sites/test_xrsite.py

 import numpy as np
 from py_wake.site.shear import PowerShear
-from py_wake.site.xrsite import XRSite
+from py_wake.site.xrsite import XRSite, GlobalWindAtlasSite
 import xarray as xr
 from py_wake.tests import npt
 import pytest
@@ -190,6 +190,31 @@ def test_complex_grid_local_wind(complex_grid_site):
                                          [0.01079997, 0.01656828, 0.02257487]])
 
 
+def test_GlobalWindAtlasSite():
+    ref = Hornsrev1Site()
+    lat, long = 55.52972, 7.906111  # hornsrev
+
+    site = GlobalWindAtlasSite(lat, long, height=70, roughness=0.001, ti=0.075)
+    ref_mean = weibull.mean(ref.ds.Weibull_A, ref.ds.Weibull_k)
+    gwa_mean = weibull.mean(site.ds.Weibull_A, site.ds.Weibull_k)
+
+    if 0:
+        plt.figure()
+        plt.plot(ref.ds.wd, ref_mean, label='HornsrevSite')
+        plt.plot(site.ds.wd, gwa_mean, label='HornsrevSite')
+        for r in [0, 1.5]:
+            for h in [10, 200]:
+                A, k = [site.gwc_ds[v].sel(roughness=r, height=h) for v in ['Weibull_A', 'Weibull_k']]
+                plt.plot(site.gwc_ds.wd, weibull.mean(A, k), label=f'{h}, {r}')
+        plt.legend()
+
+        plt.show()
+
+    npt.assert_allclose(gwa_mean, ref_mean, atol=1.4)
+    for var, atol in [('Sector_frequency', 0.03), ('Weibull_A', 1.6), ('Weibull_k', 0.4)]:
+        npt.assert_allclose(site.ds[var], ref.ds[var], atol=atol)
+
+
 def test_wrong_height():
     ti = 0.1
     ds = xr.Dataset(

py_wake/utils/weibull.py

 import numpy as np
+from scipy.special import gamma
+
+
+def mean(A, k):
+    return A * gamma(1 + 1 / k)
 
 
 def cdf(ws, A, k):