From 852ccceec53c6bd43f68e75c4cded9dde174fc4c Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Tue, 22 Jun 2021 09:54:01 +0000
Subject: [PATCH] new site based on global wind atlas data
---
py_wake/site/xrsite.py | 52 +++++++++++++++++++++++--
py_wake/tests/test_sites/test_xrsite.py | 27 ++++++++++++-
py_wake/utils/weibull.py | 5 +++
3 files changed, 80 insertions(+), 4 deletions(-)
diff --git a/py_wake/site/xrsite.py b/py_wake/site/xrsite.py
index af53d989c..db92dca22 100644
--- a/py_wake/site/xrsite.py
+++ b/py_wake/site/xrsite.py
@@ -1,9 +1,11 @@
-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
diff --git a/py_wake/tests/test_sites/test_xrsite.py b/py_wake/tests/test_sites/test_xrsite.py
index f095949b7..3f7c9c539 100644
--- a/py_wake/tests/test_sites/test_xrsite.py
+++ b/py_wake/tests/test_sites/test_xrsite.py
@@ -1,6 +1,6 @@
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(
diff --git a/py_wake/utils/weibull.py b/py_wake/utils/weibull.py
index d6714e658..0b3998906 100644
--- a/py_wake/utils/weibull.py
+++ b/py_wake/utils/weibull.py
@@ -1,4 +1,9 @@
import numpy as np
+from scipy.special import gamma
+
+
+def mean(A, k):
+ return A * gamma(1 + 1 / k)
def cdf(ws, A, k):
--
GitLab