Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • toolbox/WindEnergyToolbox
  • tlbl/WindEnergyToolbox
  • cpav/WindEnergyToolbox
  • frza/WindEnergyToolbox
  • borg/WindEnergyToolbox
  • mmpe/WindEnergyToolbox
  • ozgo/WindEnergyToolbox
  • dave/WindEnergyToolbox
  • mmir/WindEnergyToolbox
  • wluo/WindEnergyToolbox
  • welad/WindEnergyToolbox
  • chpav/WindEnergyToolbox
  • rink/WindEnergyToolbox
  • shfe/WindEnergyToolbox
  • shfe1/WindEnergyToolbox
  • acdi/WindEnergyToolbox
  • angl/WindEnergyToolbox
  • wliang/WindEnergyToolbox
  • mimc/WindEnergyToolbox
  • wtlib/WindEnergyToolbox
  • cmos/WindEnergyToolbox
  • fabpi/WindEnergyToolbox
22 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 487 additions and 236 deletions
wetb/wind/tests/test_Shear.py
View file @ 62f0f2e4
...@@ -3,13 +3,6 @@ Created on 05/06/2012 ...@@ -3,13 +3,6 @@ Created on 05/06/2012
@author: Mads @author: Mads
''' '''
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import zip
from future import standard_library
standard_library.install_aliases()
import os import os
from wetb.wind.utils import xyz2uvw from wetb.wind.utils import xyz2uvw
import wetb.gtsdf import wetb.gtsdf
......
wetb/wind/tests/test_weibull.py
View file @ 62f0f2e4
...@@ -8,8 +8,8 @@ import matplotlib.pyplot as plt ...@@ -8,8 +8,8 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
from wetb.wind import weibull from wetb.wind import weibull
class TestWeibull(unittest.TestCase):
class TestWeibull(unittest.TestCase):
def test_weibull(self): def test_weibull(self):
wb = weibull.pdf(4, 2) wb = weibull.pdf(4, 2)
...@@ -21,37 +21,39 @@ class TestWeibull(unittest.TestCase): ...@@ -21,37 +21,39 @@ class TestWeibull(unittest.TestCase):
def test_random_weibull(self): def test_random_weibull(self):
y = weibull.random(4, 2, 1000000) y = weibull.random(4, 2, 1000000)
pdf, x = np.histogram(y, 100, normed=True) pdf, x = np.histogram(y, 100, density=True)
x = (x[1:] + x[:-1]) / 2 x = (x[1:] + x[:-1]) / 2
self.assertLess(sum((pdf - weibull.pdf(4, 2)(x)) ** 2), 0.0001) self.assertLess(sum((pdf - weibull.pdf(4, 2)(x)) ** 2), 0.0001)
if 0: if 0:
plt.plot(x, pdf) plt.plot(x, pdf)
plt.plot(x, weibull.pdf(4, 2)(x)) plt.plot(x, weibull.pdf(4, 2)(x))
plt.show() plt.show()
def test_fit_weibull(self): def test_fit_weibull(self):
np.random.seed(1)
y = weibull.random(4, 2, 1000000) y = weibull.random(4, 2, 1000000)
A, k = weibull.fit(y) A, k = weibull.fit(y)
self.assertAlmostEqual(A, 4, delta=0.01) self.assertAlmostEqual(A, 4, delta=0.01)
self.assertAlmostEqual(k, 2, delta=0.01) self.assertAlmostEqual(k, 2, delta=0.01)
if 0: if 0:
plt.hist(y, 100, normed=True) plt.hist(y, 100, density=True)
x = np.arange(0, 20, .1) x = np.arange(0, 20, .1)
plt.plot(x, weibull.pdf(4, 2)(x)) plt.plot(x, weibull.pdf(4, 2)(x))
plt.show() plt.show()
def test_weibull_cdf(self): def test_weibull_cdf(self):
wb = weibull.pdf(4, 2) wb = weibull.pdf(4, 2)
x = np.arange(0, 20, .01) x = np.arange(0, 20, .01)
cdf = weibull.cdf(4,2) cdf = weibull.cdf(4, 2)
self.assertEqual(cdf(999),1) self.assertEqual(cdf(999), 1)
np.testing.assert_array_almost_equal(np.cumsum(wb(x))/100, cdf(x), 3) np.testing.assert_array_almost_equal(np.cumsum(wb(x)) / 100, cdf(x), 3)
if 0: if 0:
plt.plot(x, np.cumsum(wb(x))*.01) plt.plot(x, np.cumsum(wb(x)) * .01)
plt.plot(x, cdf(x)) plt.plot(x, cdf(x))
plt.show() plt.show()
if __name__ == "__main__": if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testweibull'] #import sys;sys.argv = ['', 'Test.testweibull']
unittest.main() unittest.main()
wetb/wind/turbulence/__init__.py 0 → 100644
View file @ 62f0f2e4
wetb/wind/turbulence/constraint_file.py
View file @ 62f0f2e4
import numpy as np import numpy as np
import os import os
from wetb.hawc2.htc_file import HTCFile from wetb.hawc2.htc_file import HTCFile
class ConstraintFile(object): class ConstraintFile(object):
def __init__(self, center_gl_xyz, box_transport_speed, no_grid_points=(4096, 32, 32), box_size=(6000, 100, 100)): def __init__(self, center_gl_xyz, box_transport_speed, no_grid_points=(4096, 32, 32), box_size=(6000, 100, 100)):
"""Generate list of constraints for turbulence constraint simulator """Generate list of constraints for turbulence constraint simulator
...@@ -23,22 +25,21 @@ class ConstraintFile(object): ...@@ -23,22 +25,21 @@ class ConstraintFile(object):
self.box_size = box_size self.box_size = box_size
self.dxyz = [d / (n - 1) for d, n in zip(self.box_size, self.no_grid_points)] self.dxyz = [d / (n - 1) for d, n in zip(self.box_size, self.no_grid_points)]
self.constraints = {'u':[], 'v':[], 'w':[]} self.constraints = {'u': [], 'v': [], 'w': []}
def load(self, filename): def load(self, filename):
"""Load constraint from existing constraint file (usable for plotting constraints via time_series()""" """Load constraint from existing constraint file (usable for plotting constraints via time_series()"""
with open(filename) as fid: with open(filename) as fid:
lines = fid.readlines() lines = fid.readlines()
for l in lines: for l in lines:
values = l.split(";") values = l.split(";")
mx,my,mz= map(int, values[:3]) mx, my, mz = map(int, values[:3])
comp = values[3]#{"100":'u','010':'v','001':'w'}["".join(values[3:6])] comp = values[3] # {"100":'u','010':'v','001':'w'}["".join(values[3:6])]
self.constraints[comp].append((mx,my,mz,float(values[-1]))) self.constraints[comp].append((mx, my, mz, float(values[-1])))
def add_constraints(self, glpos, tuvw, subtract_mean=True, fail_outside_box=True, nearest=True): def add_constraints(self, glpos, tuvw, subtract_mean=True, fail_outside_box=True, nearest=True):
"""Add constraints to constraint file """Add constraints to constraint file
parameters parameters
---------- ----------
glpos : (float or array_like, float or array_like, float or array_like) glpos : (float or array_like, float or array_like, float or array_like)
...@@ -54,7 +55,7 @@ class ConstraintFile(object): ...@@ -54,7 +55,7 @@ class ConstraintFile(object):
if False, mann coordinates modulo N is used if False, mann coordinates modulo N is used
nearest : boolean, optional nearest : boolean, optional
if True, the wsp at the position closest to the mann grid point are applied\n if True, the wsp at the position closest to the mann grid point are applied\n
if False, the average of all the wsp, that map to the mann grid point are applied if False, the average of all the wsp, that map to the mann grid point are applied
""" """
nx, ny, nz = self.no_grid_points nx, ny, nz = self.no_grid_points
dx, dy, dz = self.dxyz dx, dy, dz = self.dxyz
...@@ -62,50 +63,53 @@ class ConstraintFile(object): ...@@ -62,50 +63,53 @@ class ConstraintFile(object):
time, u, v, w = (list(tuvw.T) + [None, None])[:4] time, u, v, w = (list(tuvw.T) + [None, None])[:4]
x, y, z = glpos x, y, z = glpos
mxs = np.round((time * self.box_transport_speed + (center_y - y)) / dx).astype(np.int) mxs = np.round((time * self.box_transport_speed + (center_y - y)) / dx).astype(int)
if not np.array(y).shape==mxs.shape: if not np.array(y).shape == mxs.shape:
y = np.zeros_like(mxs)+y y = np.zeros_like(mxs) + y
if not np.array(x).shape==mxs.shape: if not np.array(x).shape == mxs.shape:
x = np.zeros_like(mxs)+x x = np.zeros_like(mxs) + x
if not np.array(z).shape==mxs.shape: if not np.array(z).shape == mxs.shape:
z = np.zeros_like(mxs)+z z = np.zeros_like(mxs) + z
x_err = y-(-(mxs*dx-time*self.box_transport_speed-center_y)) x_err = y - (-(mxs * dx - time * self.box_transport_speed - center_y))
mys = np.round(((ny - 1) / 2. + (-x+center_x) / dy)).astype(np.int) mys = np.round(((ny - 1) / 2. + (-x + center_x) / dy)).astype(int)
mzs = np.round(((nz - 1) / 2. + (center_z - z) / dz)).astype(np.int) mzs = np.round(((nz - 1) / 2. + (center_z - z) / dz)).astype(int)
y_err = x-(-((mys-(ny-1)/2)*dy-center_x)) y_err = x - (-((mys - (ny - 1) / 2) * dy - center_x))
z_err = z-(-(mzs-(nz-1)/2)*dz+center_z) z_err = z - (-(mzs - (nz - 1) / 2) * dz + center_z)
pos_err = np.sqrt(x_err**2+y_err**2+z_err**2) pos_err = np.sqrt(x_err**2 + y_err**2 + z_err**2)
if fail_outside_box: if fail_outside_box:
for ms, n in zip([mxs,mys,mzs],self.no_grid_points): for ms, n in zip([mxs, mys, mzs], self.no_grid_points):
if ms.min() + 1 < 1: if ms.min() + 1 < 1:
i = np.argmin(ms) i = np.argmin(ms)
raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz)) raise ValueError(
"At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" %
(time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx, ny, nz))
if ms.max() + 1 > n: if ms.max() + 1 > n:
i = np.argmax(ms) i = np.argmax(ms)
raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz)) raise ValueError(
"At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" %
(time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx, ny, nz))
else: else:
mxs%=nx mxs %= nx
mys%=ny mys %= ny
mzs%=nz mzs %= nz
mann_pos = ["%06d%02d%02d"%(mx,my,mz) for mx,my,mz in zip(mxs,mys,mzs)] mann_pos = ["%06d%02d%02d" % (mx, my, mz) for mx, my, mz in zip(mxs, mys, mzs)]
i = 0 i = 0
N = len(mann_pos) N = len(mann_pos)
uvw_lst = [uvw for uvw in [u,v,w] if uvw is not None] uvw_lst = [uvw for uvw in [u, v, w] if uvw is not None]
if subtract_mean: if subtract_mean:
uvw_lst = [uvw - uvw.mean() for uvw in uvw_lst] uvw_lst = [uvw - uvw.mean() for uvw in uvw_lst]
while i<N: while i < N:
for j in range(i+1,N+1): for j in range(i + 1, N + 1):
if j==N or mann_pos[i]!=mann_pos[j]: if j == N or mann_pos[i] != mann_pos[j]:
break break
for comp, wsp in zip(['u', 'v', 'w'], uvw_lst): for comp, wsp in zip(['u', 'v', 'w'], uvw_lst):
if nearest: if nearest:
value = wsp[i:j][np.argmin(pos_err[i:j])] value = wsp[i:j][np.argmin(pos_err[i:j])]
else: #mean else: # mean
value = wsp[i:j].mean() value = wsp[i:j].mean()
self.constraints[comp].append((mxs[i] + 1, mys[i] + 1, mzs[i] + 1, value)) self.constraints[comp].append((mxs[i] + 1, mys[i] + 1, mzs[i] + 1, value))
i=j i = j
def __str__(self): def __str__(self):
return "\n".join(["%d;%d;%d;%s;%.10f" % (mx, my, mz, comp, value) for comp in ['u', 'v', 'w'] return "\n".join(["%d;%d;%d;%s;%.10f" % (mx, my, mz, comp, value) for comp in ['u', 'v', 'w']
...@@ -117,15 +121,14 @@ class ConstraintFile(object): ...@@ -117,15 +121,14 @@ class ConstraintFile(object):
os.makedirs(path, exist_ok=True) os.makedirs(path, exist_ok=True)
with open(os.path.join(path, "%s.con" % name), 'w') as fid: with open(os.path.join(path, "%s.con" % name), 'w') as fid:
fid.write(str(self)) fid.write(str(self))
def time_series(self, y_position=0): def time_series(self, y_position=0):
time = (np.arange(self.no_grid_points[0])*float(self.dxyz[0]) + y_position)/ self.box_transport_speed time = (np.arange(self.no_grid_points[0]) * float(self.dxyz[0]) + y_position) / self.box_transport_speed
u,v,w = [(np.zeros_like(time)+ np.nan)]*3 u, v, w = [(np.zeros_like(time) + np.nan)] * 3
for uvw, constr in zip([u,v,w], [np.array(self.constraints[uvw]) for uvw in 'uvw']): for uvw, constr in zip([u, v, w], [np.array(self.constraints[uvw]) for uvw in 'uvw']):
if constr.shape[0]>0: if constr.shape[0] > 0:
uvw[constr[:,0].astype(np.int)-1] =constr[:,3] uvw[constr[:, 0].astype(int) - 1] = constr[:, 3]
return time, np.array([u,v,w]).T return time, np.array([u, v, w]).T
def simulation_cmd(self, ae23, L, Gamma, seed, name, folder="./constraints/"): def simulation_cmd(self, ae23, L, Gamma, seed, name, folder="./constraints/"):
assert isinstance(seed, int) and seed > 0, "seed must be a positive integer" assert isinstance(seed, int) and seed > 0, "seed must be a positive integer"
...@@ -133,17 +136,13 @@ class ConstraintFile(object): ...@@ -133,17 +136,13 @@ class ConstraintFile(object):
cmd += "%.6f %.2f %.2f %d " % (ae23, L, Gamma, seed) cmd += "%.6f %.2f %.2f %d " % (ae23, L, Gamma, seed)
cmd += "./turb/%s_s%s %s/%s.con" % (name, seed, folder, name) cmd += "./turb/%s_s%s %s/%s.con" % (name, seed, folder, name)
return cmd return cmd
def hawc2_mann_section(self, name, seed): def hawc2_mann_section(self, name, seed):
htc = HTCFile() htc = HTCFile()
mann = htc.wind.add_section("mann") mann = htc.wind.add_section("mann")
for uvw in 'uvw': for uvw in 'uvw':
setattr(mann, 'filename_%s'%uvw, "./turb/%s_s%04d%s.bin" % ( name, int(seed),uvw)) setattr(mann, 'filename_%s' % uvw, "./turb/%s_s%04d%s.bin" % (name, int(seed), uvw))
for uvw, l,n in zip('uvw', self.box_size, self.no_grid_points): for uvw, l, n in zip('uvw', self.box_size, self.no_grid_points):
setattr(mann, 'box_dim_%s'%uvw, [n,l/n]) setattr(mann, 'box_dim_%s' % uvw, [n, l / n])
mann.dont_scale = 1 mann.dont_scale = 1
return mann return mann
wetb/wind/turbulence/mann_parameters.py
View file @ 62f0f2e4
...@@ -14,8 +14,6 @@ from wetb.wind.turbulence.spectra import spectra, logbin_spectra, plot_spectra,\ ...@@ -14,8 +14,6 @@ from wetb.wind.turbulence.spectra import spectra, logbin_spectra, plot_spectra,\
detrend_wsp detrend_wsp
sp1, sp2, sp3, sp4 = np.load(os.path.dirname(__file__).replace("library.zip", '') + "/mann_spectra_data.npy") sp1, sp2, sp3, sp4 = np.load(os.path.dirname(__file__).replace("library.zip", '') + "/mann_spectra_data.npy")
yp = np.arange(-3, 3.1, 0.1) yp = np.arange(-3, 3.1, 0.1)
...@@ -26,8 +24,7 @@ RBS3 = RectBivariateSpline(xp, yp, sp3) ...@@ -26,8 +24,7 @@ RBS3 = RectBivariateSpline(xp, yp, sp3)
RBS4 = RectBivariateSpline(xp, yp, sp4) RBS4 = RectBivariateSpline(xp, yp, sp4)
# def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
#def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
# if isinstance(fs, (int, float)): # if isinstance(fs, (int, float)):
# fs = [fs] * len(u_lst) # fs = [fs] * len(u_lst)
# if v_lst is None: # if v_lst is None:
...@@ -41,10 +38,6 @@ RBS4 = RectBivariateSpline(xp, yp, sp4) ...@@ -41,10 +38,6 @@ RBS4 = RectBivariateSpline(xp, yp, sp4)
# return [np.mean(np.array([kuvw[i] for kuvw in uvw_lst]), 0) for i in range(5)] # return [np.mean(np.array([kuvw[i] for kuvw in uvw_lst]), 0) for i in range(5)]
def get_mann_model_spectra(ae, L, G, k1): def get_mann_model_spectra(ae, L, G, k1):
"""Mann model spectra """Mann model spectra
...@@ -73,15 +66,11 @@ def get_mann_model_spectra(ae, L, G, k1): ...@@ -73,15 +66,11 @@ def get_mann_model_spectra(ae, L, G, k1):
xq = np.log10(L * k1) xq = np.log10(L * k1)
yq = (np.zeros_like(xq) + G) yq = (np.zeros_like(xq) + G)
f = L ** (5 / 3) * ae f = L ** (5 / 3) * ae
uu = f * RBS1.ev(yq, xq) m = (yq >= 0) & (yq <= 5) & (xq >= -3) & (xq <= 3)
vv = f * RBS2.ev(yq, xq) uu, vv, ww, uw = [f * np.where(m, RBS.ev(yq, xq), np.nan) for RBS in [RBS1, RBS2, RBS3, RBS4]]
ww = f * RBS3.ev(yq, xq)
uw = f * RBS4.ev(yq, xq)
return uu, vv, ww, uw return uu, vv, ww, uw
def _local_error(x, k1, uu, vv, ww=None, uw=None): def _local_error(x, k1, uu, vv, ww=None, uw=None):
ae, L, G = x ae, L, G = x
...@@ -95,7 +84,9 @@ def _local_error(x, k1, uu, vv, ww=None, uw=None): ...@@ -95,7 +84,9 @@ def _local_error(x, k1, uu, vv, ww=None, uw=None):
val += np.sum((k1 * ww - k1 * tmpww) ** 2) + np.sum((k1 * uw - k1 * tmpuw) ** 2) val += np.sum((k1 * ww - k1 * tmpww) ** 2) + np.sum((k1 * uw - k1 * tmpuw) ** 2)
return val return val
def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
"""Fit a mann model to the spectra """Fit a mann model to the spectra
Bins the spectra, into logarithmic sized bins and find the mann model parameters, Bins the spectra, into logarithmic sized bins and find the mann model parameters,
...@@ -139,7 +130,8 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, ...@@ -139,7 +130,8 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
>>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v)) >>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v))
""" """
from scipy.optimize import fmin from scipy.optimize import fmin
x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False) x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(
k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False)
if plt: if plt:
if not hasattr(plt, 'plot'): if not hasattr(plt, 'plot'):
...@@ -147,7 +139,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, ...@@ -147,7 +139,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
# plot_spectra(k1, uu, vv, ww, uw, plt=plt) # plot_spectra(k1, uu, vv, ww, uw, plt=plt)
# plot_mann_spectra(*x, plt=plt) # plot_mann_spectra(*x, plt=plt)
ae, L, G = x ae, L, G = x
plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt) plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt)
plt.title('ae:%.3f, L:%.1f, G:%.2f' % tuple(x)) plt.title('ae:%.3f, L:%.1f, G:%.2f' % tuple(x))
plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]') plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$') plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
...@@ -155,6 +147,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, ...@@ -155,6 +147,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
plt.show() plt.show()
return x return x
def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2): def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
"""Fit a mann model to the spectra """Fit a mann model to the spectra
...@@ -198,40 +191,59 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2): ...@@ -198,40 +191,59 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1)) return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
def var2ae(variance, spatial_resolution, N, L, G): def var2ae(variance, L, G, U, T=600, sample_frq=10, plt=False):
"""Fit alpha-epsilon to match variance of time series """Fit alpha-epsilon to match variance of time series
Parameters Parameters
---------- ----------
variance : array-like variance : array-like
variance of u vind component variance of u vind component
spatial_resolution : int, float or array_like
Distance between samples in meters
- For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U
N : int
L : int or float L : int or float
Length scale of Mann model Length scale of Mann model
G : int or float G : int or float
Gamma of Mann model Gamma of Mann model
U : int or float
Mean wind speed
T: int or float
Length [s] of signal, from which the variance is calculated
sample_frq: int or float
Sample frequency [Hz] of signal from which the variance is calculated
Returns Returns
------- -------
ae : float ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u Alpha epsilon^(2/3) of Mann model that makes the energy of the model in the
frequency range [1/length, sample_frq] equal to the variance of u
""" """
k1 = np.logspace(1,10,1000)/100000000 k_low, k_high = 2 * np.pi / (U * np.array([T, 2 / sample_frq]))
k1 = 10 ** (np.linspace(np.log10(k_low), np.log10(k_high), 1000))
def get_var(uu): def get_var(uu):
return np.trapz(2 * uu[:], k1[:]) return np.trapz(2 * uu[:], k1[:])
v1 = get_var(get_mann_model_spectra(0.1, L, G, k1)[0]) v1 = get_var(get_mann_model_spectra(0.1, L, G, k1)[0])
v2 = get_var(get_mann_model_spectra(0.2, L, G, k1)[0]) v2 = get_var(get_mann_model_spectra(0.2, L, G, k1)[0])
ae = (variance - v1) / (v2 - v1) * .1 + .1 ae = (variance - v1) / (v2 - v1) * .1 + .1
if plt is not False:
if not hasattr(plt, 'plot'):
import matplotlib.pyplot as plt
muu = get_mann_model_spectra(ae, L, G, k1)[0]
plt.semilogx(k1, k1 * muu, label='ae:%.3f, L:%.1f, G:%.2f' % (ae, L, G))
plt.legend()
plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
return ae return ae
def ae2ti(ae23, L, G, U, T=600, sample_frq=10):
k_low, k_high = 2 * np.pi / (U * np.array([T, 2 / sample_frq]))
k1 = 10 ** (np.linspace(np.log10(k_low), np.log10(k_high), 1000))
uu = get_mann_model_spectra(ae23, L, G, k1)[0]
var = np.trapz(2 * uu[:], k1[:])
return np.sqrt(var) / U
def fit_ae(spatial_resolution, u, L, G, plt=False): def fit_ae(spatial_resolution, u, L, G, plt=False):
"""Fit alpha-epsilon to match variance of time series """Fit alpha-epsilon to match variance of time series
...@@ -239,8 +251,8 @@ def fit_ae(spatial_resolution, u, L, G, plt=False): ...@@ -239,8 +251,8 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
Parameters Parameters
---------- ----------
spatial_resolution : int, float or array_like spatial_resolution : int, float or array_like
Distance between samples in meters Number of points pr meterDistance between samples in meters
- For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points, - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
Nx is number of points and Lx is box length in meters Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U - For time series: Sample frequency / U
u : array-like u : array-like
...@@ -255,22 +267,22 @@ def fit_ae(spatial_resolution, u, L, G, plt=False): ...@@ -255,22 +267,22 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
ae : float ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
""" """
#if len(u.shape) == 1: # if len(u.shape) == 1:
# u = u.reshape(len(u), 1) # u = u.reshape(len(u), 1)
# if min_bin_count is None: # if min_bin_count is None:
# min_bin_count = max(2, 6 - u.shape[0] / 2) # min_bin_count = max(2, 6 - u.shape[0] / 2)
# min_bin_count = 1 # min_bin_count = 1
def get_var(k1, uu): def get_var(k1, uu):
l = 0 #128 // np.sqrt(u.shape[1]) l = 0 # 128 // np.sqrt(u.shape[1])
return np.trapz(2 * uu[l:], k1[l:]) return np.mean(np.trapz(2 * uu[l:], k1[l:], axis=0))
k1, uu = spectra(spatial_resolution, u)[:2] k1, uu = spectra(spatial_resolution, u)[:2]
v = get_var(k1,uu) v = get_var(k1, uu)
v1 = get_var(k1, get_mann_model_spectra(0.1, L, G, k1)[0]) v1 = get_var(k1, get_mann_model_spectra(0.1, L, G, k1)[0])
v2 = get_var(k1, get_mann_model_spectra(0.2, L, G, k1)[0]) v2 = get_var(k1, get_mann_model_spectra(0.2, L, G, k1)[0])
ae = (v - v1) / (v2 - v1) * .1 + .1 ae = (v - v1) / (v2 - v1) * .1 + .1
# print (ae) # print (ae)
# #
# k1 = spectra(sf, u)[0] # k1 = spectra(sf, u)[0]
# v1 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.1, L, G, k1)[0], min_bin_count=min_bin_count)[:2]) # v1 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.1, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
# v2 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.2, L, G, k1)[0], min_bin_count=min_bin_count)[:2]) # v2 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.2, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
...@@ -283,7 +295,7 @@ def fit_ae(spatial_resolution, u, L, G, plt=False): ...@@ -283,7 +295,7 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
if not hasattr(plt, 'plot'): if not hasattr(plt, 'plot'):
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.semilogx(k1, k1 * uu, 'b-', label='uu') plt.semilogx(k1, k1 * uu, 'b-', label='uu')
k1_lb, uu_lb = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2] k1_lb, uu_lb = logbin_spectra(*spectra(spatial_resolution, u), min_bin_count=1)[:2]
plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin') plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin')
muu = get_mann_model_spectra(ae, L, G, k1)[0] muu = get_mann_model_spectra(ae, L, G, k1)[0]
...@@ -296,24 +308,27 @@ def fit_ae(spatial_resolution, u, L, G, plt=False): ...@@ -296,24 +308,27 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None): def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
# if plt is None: # if plt is None:
# import matplotlib.pyplot as plt # import matplotlib.pyplot as plt
plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt) plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
plot_mann_spectra(ae, L, G, "-", mean_u, plt) plot_mann_spectra(ae, L, G, "-", mean_u, plt)
def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']): def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
if plt is None: if plt is None:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
mf = 10 ** (np.linspace(-4, 1, 1000)) mf = 10 ** (np.linspace(-4, 3, 1000))
mf = 10 ** (np.linspace(-4, 3, 1000000))
muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf) muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G)) plt.title("ae: %.3f, L: %.2f, G:%.2f" % (ae, L, G))
if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style) if 'uu' in spectra:
if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style) plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style) if 'vv' in spectra:
if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style) plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
if 'ww' in spectra:
plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
if 'uw' in spectra:
plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
if __name__ == "__main__": if __name__ == "__main__":
...@@ -323,30 +338,29 @@ if __name__ == "__main__": ...@@ -323,30 +338,29 @@ if __name__ == "__main__":
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
"""Example of fitting Mann parameters to a time series""" """Example of fitting Mann parameters to a time series"""
ds = gtsdf.Dataset(os.path.dirname(wind.__file__)+"/tests/test_files/WspDataset.hdf5")#'unit_test/test_files/wspdataset.hdf5') ds = gtsdf.Dataset(os.path.dirname(wind.__file__) +
"/tests/test_files/WspDataset.hdf5") # 'unit_test/test_files/wspdataset.hdf5')
f = 35 f = 35
u, v = wsp_dir2uv(ds.Vhub_85m, ds.Dir_hub_) u, v = wsp_dir2uv(ds.Vhub_85m, ds.Dir_hub_)
u_ref = np.mean(u) u_ref = np.mean(u)
u -= u_ref u -= u_ref
sf = f / u_ref sf = f / u_ref
ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt = None) ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt=None)
print (ae, L, G) print(ae, L, G)
print (u.shape) print(u.shape)
print (ds.Time[-1]) print(ds.Time[-1])
plt.plot(u) plt.plot(u)
plt.plot(detrend_wsp(u)[0]) plt.plot(detrend_wsp(u)[0])
plt.show() plt.show()
print (fit_ae(sf, detrend_wsp(u)[0], L, G, plt)) print(fit_ae(sf, detrend_wsp(u)[0], L, G, plt))
print (var2ae(detrend_wsp(u)[0].var(), L, G,)) print(var2ae(detrend_wsp(u)[0].var(), L, G,))
print () print()
print (fit_ae(sf, u[:21000], L, G)) print(fit_ae(sf, u[:21000], L, G))
print (var2ae(u[:21000].var(), L, G,)) print(var2ae(u[:21000].var(), L, G,))
# """Example of fitting Mann parameters to a "series" of a turbulence box""" # """Example of fitting Mann parameters to a "series" of a turbulence box"""
# l = 16384 # l = 16384
# nx = 8192 # nx = 8192
...@@ -359,7 +373,6 @@ if __name__ == "__main__": ...@@ -359,7 +373,6 @@ if __name__ == "__main__":
# print (ae, L, G) # print (ae, L, G)
# """Example of fitting Mann parameters to a "series" of a turbulence box""" # """Example of fitting Mann parameters to a "series" of a turbulence box"""
# l = 4778.3936 # l = 4778.3936
# nx = 8192 # nx = 8192
...@@ -373,4 +386,4 @@ if __name__ == "__main__": ...@@ -373,4 +386,4 @@ if __name__ == "__main__":
# #ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt) # #ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt)
# #print (fit_ae(sf, u, 73.0730383576, 2.01636095317)) # #print (fit_ae(sf, u, 73.0730383576, 2.01636095317))
# print (u.std()) # print (u.std())
# #print (ae, L, G) # #print (ae, L, G)
\ No newline at end of file
wetb/wind/turbulence/mann_turbulence.py
View file @ 62f0f2e4
...@@ -8,38 +8,42 @@ import numpy as np ...@@ -8,38 +8,42 @@ import numpy as np
from wetb.wind.turbulence.spectra import spectra, spectra_from_time_series from wetb.wind.turbulence.spectra import spectra, spectra_from_time_series
name_format = "mann_l%.1f_ae%.4f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb" name_format = "mann_l%.1f_ae%.4f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
def load(filename, N=(32,32)):
def load(filename, N=(32, 32)):
"""Load mann turbulence box """Load mann turbulence box
Parameters Parameters
---------- ----------
filename : str filename : str
Filename of turbulence box Filename of turbulence box
N : tuple, (ny,nz) or (nx,ny,nz) N : tuple, (ny,nz) or (nx,ny,nz)
Number of grid points Number of grid points
Returns Returns
------- -------
turbulence_box : nd_array turbulence_box : nd_array
Examples Examples
-------- --------
>>> u = load('turb_u.dat') >>> u = load('turb_u.dat')
""" """
data = np.fromfile(filename, np.dtype('<f'), -1) data = np.fromfile(filename, np.dtype('<f'), -1)
if len(N)==2: if len(N) == 2:
ny,nz = N ny, nz = N
nx = len(data)/(ny*nz) nx = len(data) / (ny * nz)
assert nx==int(nx), "Size of turbulence box (%d) does not match ny x nz (%d), nx=%.2f" % (len(data),ny*nz, nx) assert nx == int(nx), "Size of turbulence box (%d) does not match ny x nz (%d), nx=%.2f" % (
nx=int(nx) len(data), ny * nz, nx)
nx = int(nx)
else: else:
nx,ny,nz = N nx, ny, nz = N
assert len(data) == nx*ny*nz, "Size of turbulence box (%d) does not match nx x ny x nz (%d)" % (len(data),nx*ny*nz) assert len(data) == nx * ny * \
return data.reshape(nx , ny * nz) nz, "Size of turbulence box (%d) does not match nx x ny x nz (%d)" % (len(data), nx * ny * nz)
return data.reshape(nx, ny * nz)
def load_uvw(filenames, N=(1024,32,32)):
def load_uvw(filenames, N=(1024, 32, 32)):
"""Load u, v and w turbulence boxes """Load u, v and w turbulence boxes
Parameters Parameters
---------- ----------
filenames : list or str filenames : list or str
...@@ -47,34 +51,40 @@ def load_uvw(filenames, N=(1024,32,32)): ...@@ -47,34 +51,40 @@ def load_uvw(filenames, N=(1024,32,32)):
if str: filename pattern where u,v,w are replaced with '%s' if str: filename pattern where u,v,w are replaced with '%s'
N : tuple N : tuple
Number of grid point in the x, y and z direction Number of grid point in the x, y and z direction
Returns Returns
------- -------
u,v,w : list of np_array u,v,w : list of np_array
Examples Examples
-------- --------
>>> u,v,w =load_uvw('turb_%s.dat') >>> u,v,w =load_uvw('turb_%s.dat')
""" """
if isinstance(filenames,str): if isinstance(filenames, str):
return [load(filenames%uvw, N) for uvw in 'uvw'] return [load(filenames % uvw, N) for uvw in 'uvw']
else: else:
return [load(f, N) for f in filenames] return [load(f, N) for f in filenames]
def parameters2name(no_grid_points,box_dimension,ae23,L, Gamma,high_frq_compensation, seed, folder="./turb/"):
def save(turb, filename):
np.array(turb).astype('<f').tofile(filename)
def parameters2name(no_grid_points, box_dimension, ae23, L, Gamma, high_frq_compensation, seed, folder="./turb/"):
dxyz = tuple(np.array(box_dimension) / no_grid_points) dxyz = tuple(np.array(box_dimension) / no_grid_points)
return ["./turb/" + name_format % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']] return ["./turb/" + name_format % ((L, ae23, Gamma, high_frq_compensation) +
no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None): def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
"""Fit mann parameters, ae, L, G """Fit mann parameters, ae, L, G
Parameters Parameters
---------- ----------
spatial_resolution : inf, float or array_like spatial_resolution : inf, float or array_like
Distance between samples in meters Distance between samples in meters
- For turbulence boxes: 1/dx = Lx/Nx where dx is distance between points, - For turbulence boxes: 1/dx = Lx/Nx where dx is distance between points,
Nx is number of points and Lx is box length in meters Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U - For time series: Sample frequency / U
u : array_like u : array_like
...@@ -89,7 +99,7 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None): ...@@ -89,7 +99,7 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
The w-wind component The w-wind component
- if shape is (r,): One time series with *r* observations\n - if shape is (r,): One time series with *r* observations\n
- if shape is (r,c): *c* different time series with *r* observations\n - if shape is (r,c): *c* different time series with *r* observations\n
Returns Returns
------- -------
ae : int or float ae : int or float
...@@ -105,4 +115,4 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None): ...@@ -105,4 +115,4 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
def fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=None): def fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=None):
from wetb.wind.turbulence.mann_parameters import fit_mann_model_spectra from wetb.wind.turbulence.mann_parameters import fit_mann_model_spectra
return fit_mann_model_spectra(*spectra_from_time_series(sample_frq, Uvw_lst), plt=plt) return fit_mann_model_spectra(*spectra_from_time_series(sample_frq, Uvw_lst), plt=plt)
\ No newline at end of file
wetb/wind/turbulence/spectra.py
View file @ 62f0f2e4
...@@ -7,6 +7,8 @@ Created on 27/11/2015 ...@@ -7,6 +7,8 @@ Created on 27/11/2015
import numpy as np import numpy as np
import warnings import warnings
def spectrum(x, y=None, k=1): def spectrum(x, y=None, k=1):
"""PSD or Cross spectrum (only positive half) """PSD or Cross spectrum (only positive half)
...@@ -25,6 +27,7 @@ def spectrum(x, y=None, k=1): ...@@ -25,6 +27,7 @@ def spectrum(x, y=None, k=1):
""" """
if x is None: if x is None:
return None return None
def fft(x): def fft(x):
return np.fft.fft(x.T).T / len(x) return np.fft.fft(x.T).T / len(x)
...@@ -34,22 +37,23 @@ def spectrum(x, y=None, k=1): ...@@ -34,22 +37,23 @@ def spectrum(x, y=None, k=1):
else: else:
fftx = fft(x) * np.conj(fft(y)) # Cross spectra fftx = fft(x) * np.conj(fft(y)) # Cross spectra
fftx = fftx[:len(fftx) // 2 ] * 2 # positive half * 2 fftx = fftx[:len(fftx) // 2] * 2 # positive half * 2
# if len(fftx.shape) == 2: # if len(fftx.shape) == 2:
# fftx = np.mean(fftx, 1) # fftx = np.mean(fftx, 1)
with warnings.catch_warnings(): with warnings.catch_warnings():
warnings.simplefilter("ignore") warnings.simplefilter("ignore")
return np.real(fftx * len(x) / (2 * k))[1:] return np.real(fftx * len(x) / (2 * k))[1:]
def spectra(spatial_resolution, u, v=None, w=None, detrend=True): def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
"""Return the wave number, the uu, vv, ww autospectra and the uw cross spectra """Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
Parameters Parameters
---------- ----------
spatial_resolution : int, float or array_like spatial_resolution : int, float or array_like
Distance between samples in meters Sample points per meter
- For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points, - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
Nx is number of points and Lx is box length in meters Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U - For time series: Sample frequency / U
u : array_like u : array_like
...@@ -81,24 +85,30 @@ def spectra(spatial_resolution, u, v=None, w=None, detrend=True): ...@@ -81,24 +85,30 @@ def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
assert isinstance(spatial_resolution, (int, float)) assert isinstance(spatial_resolution, (int, float))
k = 2 * np.pi * spatial_resolution k = 2 * np.pi * spatial_resolution
if v is not None: assert u.shape==v.shape if v is not None:
if w is not None: assert u.shape==w.shape assert u.shape == v.shape
if w is not None:
if 1 and len(u.shape)==2: assert u.shape == w.shape
assert np.abs(np.mean(u,0)).max()<1
if v is not None: assert np.abs(np.mean(v,0)).max()<1 if 1 and len(u.shape) == 2:
if w is not None: assert np.abs(np.mean(w,0)).max()<1 # assert np.abs(np.mean(u, 0)).max() < 2
# if v is not None:
# assert np.abs(np.mean(v, 0)).max() < 1
# if w is not None:
# assert np.abs(np.mean(w, 0)).max() < 1
if isinstance(k, float): if isinstance(k, float):
k = np.repeat(k,u.shape[1]) k = np.repeat(k, u.shape[1])
else: else:
assert u.shape[1]==k.shape[0] assert u.shape[1] == k.shape[0]
k1_vec = np.array([np.linspace(0, k_ / 2, len(u) / 2)[1:] for k_ in k]).T k1_vec = np.array([np.linspace(0, k_ / 2, len(u) // 2)[1:] for k_ in k]).T
else: else:
assert np.abs(np.mean(u))<1 #assert np.abs(np.mean(u)) < 1
if v is not None: assert np.abs(np.mean(v))<1 if v is not None:
if w is not None: assert np.abs(np.mean(w))<1 assert np.abs(np.mean(v)) < 1
assert isinstance(k,float) if w is not None:
k1_vec = np.linspace(0, k / 2, len(u) / 2)[1:] assert np.abs(np.mean(w)) < 1
assert isinstance(k, float)
k1_vec = np.linspace(0, k / 2, int(len(u) / 2))[1:]
if detrend: if detrend:
u, v, w = detrend_wsp(u, v, w) u, v, w = detrend_wsp(u, v, w)
...@@ -114,7 +124,7 @@ def spectra_from_time_series(sample_frq, Uvw_lst): ...@@ -114,7 +124,7 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
Sample frequency Sample frequency
Uvw_lst : array_like Uvw_lst : array_like
list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...], v and w are optional list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...], v and w are optional
Returns Returns
------- -------
k1, uu[, vv[, ww, uw]] : 2-5 x array_like k1, uu[, vv[, ww, uw]] : 2-5 x array_like
...@@ -129,34 +139,28 @@ def spectra_from_time_series(sample_frq, Uvw_lst): ...@@ -129,34 +139,28 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
assert isinstance(sample_frq, (int, float)) assert isinstance(sample_frq, (int, float))
Uvw_arr = np.array(Uvw_lst) Uvw_arr = np.array(Uvw_lst)
Ncomp = Uvw_arr.shape[1] Ncomp = Uvw_arr.shape[1]
U,v,w = [Uvw_arr[:,i,:].T for i in range(Ncomp)] + [None]*(3-Ncomp) U, v, w = [Uvw_arr[:, i, :].T for i in range(Ncomp)] + [None] * (3 - Ncomp)
k = 2 * np.pi * sample_frq / U.mean(0) k = 2 * np.pi * sample_frq / U.mean(0)
if v is not None: assert np.abs(np.nanmean(v,0)).max()<1, "Max absolute mean of v is %f"%np.abs(np.nanmean(v,0)).max() if v is not None:
if w is not None: assert np.abs(np.nanmean(w,0)).max()<1 assert np.abs(np.nanmean(v, 0)).max() < 1, "Max absolute mean of v is %f" % np.abs(np.nanmean(v, 0)).max()
k1_vec = np.array([np.linspace(0, k_ / 2, U.shape[0] / 2)[1:] for k_ in k]).T if w is not None:
u = U-np.nanmean(U, 0) assert np.abs(np.nanmean(w, 0)).max() < 1
k1_vec = np.array([np.linspace(0, k_ / 2, U.shape[0] // 2)[1:] for k_ in k]).T
u = U - np.nanmean(U, 0)
u, v, w = detrend_wsp(u, v, w) u, v, w = detrend_wsp(u, v, w)
return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]] return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]]
{'ns':( 298, 0.0455033341592, 30.7642421893, 2.32693322102, 0.00538254487239),
'nu':( 853, 0.0355868134195, 71.0190181051, 2.95576644546, 0.00268903502954),
's' :( 367, 0.0300454276184, 24.7375807569, 2.40165270878, 0.00253338624678),
'vs':( 223, 0.0149281609201, 11.7467239752, 3.21842435078, 0.00143453350246),
'n' :( 984, 0.0494970534618, 47.3846412548, 2.83466101838, 0.00513954022345),
'u' :( 696, 0.0258456660845, 78.9511607824, 2.51128584334, 0.00250014140782),
'vu':( 468, 0.0218274041889, 73.0730383576, 2.01636095317, 0.00196337613117)}
def bin_spectrum(x, y, bin_size, min_bin_count=2): def bin_spectrum(x, y, bin_size, min_bin_count=2):
assert min_bin_count > 0 assert min_bin_count > 0
x = x / bin_size x = x / bin_size
low, high = np.floor(np.nanmin(x)), np.ceil(np.nanmax(x)) low, high = np.floor(np.nanmin(x)), np.ceil(np.nanmax(x))
bins = int(high - low) bins = int(high - low)
nbr_in_bins = np.histogram(x, bins, range=(low, high))[0] nbr_in_bins = np.histogram(x, bins, range=(low, high))[0]
if len(x.shape)==2: if len(x.shape) == 2:
min_bin_count*=x.shape[1] min_bin_count *= x.shape[1]
mask = nbr_in_bins >= min_bin_count mask = nbr_in_bins >= min_bin_count
return np.histogram(x, bins, range=(low, high), weights=y)[0][mask] / nbr_in_bins[mask], nbr_in_bins return np.histogram(x, bins, range=(low, high), weights=y)[0][mask] / nbr_in_bins[mask], nbr_in_bins
...@@ -171,10 +175,11 @@ def logbin_spectrum(k1, xx, log10_bin_size=.2, min_bin_count=2): ...@@ -171,10 +175,11 @@ def logbin_spectrum(k1, xx, log10_bin_size=.2, min_bin_count=2):
def logbin_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=0.2, min_bin_count=2): def logbin_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=0.2, min_bin_count=2):
return tuple([logbin_spectrum(k1, xx, log10_bin_size, min_bin_count) for xx in [k1, uu, vv, ww, uw]]) return tuple([logbin_spectrum(k1, xx, log10_bin_size, min_bin_count) for xx in [k1, uu, vv, ww, uw]])
def plot_spectrum(spacial_frq, u, plt=None):
def plot_spectrum(spatial_resolution, u, plt=None):
if plt is None: if plt is None:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
k1, uu = logbin_spectra(*spectra(spacial_frq, u))[:2] k1, uu = logbin_spectra(*spectra(spatial_resolution, u))[:2]
plt.semilogx(k1, k1 * uu, 'b-') plt.semilogx(k1, k1 * uu, 'b-')
...@@ -187,7 +192,8 @@ def detrend_wsp(u, v=None, w=None): ...@@ -187,7 +192,8 @@ def detrend_wsp(u, v=None, w=None):
t = np.arange(dwsp.shape[0]) t = np.arange(dwsp.shape[0])
A = np.vstack([t, np.ones(len(t))]).T A = np.vstack([t, np.ones(len(t))]).T
for i in range(dwsp.shape[1]): for i in range(dwsp.shape[1]):
trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0] m = ~np.isnan(dwsp[:, i])
trend, offset = np.linalg.lstsq(A[:m.sum()], dwsp[:, i][m], rcond=None)[0]
dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
return dwsp.reshape(wsp.shape) return dwsp.reshape(wsp.shape)
return [_detrend(wsp) for wsp in [u, v, w]] return [_detrend(wsp) for wsp in [u, v, w]]
...@@ -197,11 +203,12 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, ...@@ -197,11 +203,12 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2,
if plt is None: if plt is None:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size) bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
def plot(xx, label, color, plt): def plot(xx, label, color, plt):
plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , marker_style + color, label=label) plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2, marker_style + color, label=label)
plot(buu, 'uu', 'r', plt) plot(buu, 'uu', 'r', plt)
plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]') plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
if mean_u ==1: if mean_u == 1:
plt.ylabel(r'Spectral density $k_{1} F(k_{1}) [m^2/s^2]$') plt.ylabel(r'Spectral density $k_{1} F(k_{1}) [m^2/s^2]$')
else: else:
plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [-]$') plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [-]$')
...@@ -209,4 +216,4 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, ...@@ -209,4 +216,4 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2,
plot(bvv, 'vv', 'g', plt) plot(bvv, 'vv', 'g', plt)
if bww is not None: if bww is not None:
plot(bww, 'ww', 'b', plt) plot(bww, 'ww', 'b', plt)
plot(buw, 'uw', 'm', plt) plot(buw, 'uw', 'm', plt)
\ No newline at end of file
wetb/wind/turbulence/tests/test_man_turbulence.py 0 → 100644
View file @ 62f0f2e4
'''
Created on 20. jul. 2017
@author: mmpe
'''
import unittest
from wetb.utils.test_files import get_test_file
from wetb.wind.turbulence import mann_turbulence
from wetb.wind.turbulence.mann_turbulence import parameters2name
from tests import npt
class TestMannTurbulence(unittest.TestCase):
def testLoad(self):
fn = get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat')
self.assertRaises(AssertionError, mann_turbulence.load, fn, (31, 31))
self.assertRaises(AssertionError, mann_turbulence.load, fn, (8192, 32, 32))
u = mann_turbulence.load(fn, (8192, 8, 8))
self.assertEqual(u.shape, (8192, 8 * 8))
def test_loaduvw(self):
fn = 'h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat'
fn_lst = [get_test_file(fn % uvw) for uvw in 'uvw']
u, v, w = mann_turbulence.load_uvw(fn_lst, (8192, 8, 8))
self.assertEqual(u.shape, (8192, 8 * 8))
self.assertTrue(u.shape == v.shape == w.shape == (8192, 8 * 8))
u, v, w = mann_turbulence.load_uvw(fn_lst[0].replace("u.dat", "%s.dat"), (8192, 8, 8))
self.assertEqual(u.shape, (8192, 8 * 8))
self.assertTrue(u.shape == v.shape == w.shape == (8192, 8 * 8))
def test_parameters2name(self):
self.assertEqual(parameters2name((8192, 8, 8), (16384, 32, 32), 0.15, 10, 3.3, 1, 1)[
0], './turb/mann_l10.0_ae0.1500_g3.3_h1_8192x8x8_2.000x4.00x4.00_s0001u.turb')
def test_save(self):
fn = get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat')
u_ref = mann_turbulence.load(fn, (8192, 8, 8))
mann_turbulence.save(u_ref, 'tmp_u.turb')
u = mann_turbulence.load('tmp_u.turb', (8192, 8, 8))
npt.assert_array_equal(u, u_ref)
# def test_name2parameters(self):
# self.assertEqual(name2parameters('./turb/mann_l10.0_ae0.15_g3.3_h1_8192x8x8_2.000x4.00x4.00_s0001u.turb'),
# ((8192, 8, 8), (16384., 32., 32.), 0.15, 10, 3.3, 1, 1))
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
wetb/wind/turbulence/tests/test_mann_parameters.py 0 → 100644
View file @ 62f0f2e4
'''
Created on 05/06/2012
@author: Mads
'''
import os
import sys
import unittest
import matplotlib.pyplot as plt
import numpy as np
from wetb.utils.test_files import get_test_file, move2test_files
from wetb.wind.turbulence.mann_turbulence import load_uvw, fit_mann_parameters,\
fit_mann_parameters_from_time_series
from wetb.utils.timing import print_time
from wetb.wind.turbulence.spectra import spectra_from_time_series, spectra,\
plot_spectra, logbin_spectra
import warnings
from wetb.wind.turbulence import mann_turbulence
from wetb.wind.turbulence.mann_parameters import var2ae, fit_ae, ae2ti
from numpy import spacing
tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
class TestMannTurbulence(unittest.TestCase):
def test_fit_mann_parameters_turbulence_box(self):
# for uvw in 'uvw':
# move2test_files(tfp + 'h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat'%uvw)
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
plt = None
ae, L, G = fit_mann_parameters(fx, u, v, w, plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
def test_spectra_from_timeseries(self):
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
k1, uu, vv, ww, uw = logbin_spectra(*spectra(fx, u, v, w))
U = u + 4
sample_frq = 2
k12, uu2, vv2, ww2, uw2 = logbin_spectra(
*spectra_from_time_series(sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
np.testing.assert_allclose(uu, uu2, 0.02)
U = u + 8
sample_frq = 2
k13, uu3, vv3, ww3, uw3 = logbin_spectra(
*spectra_from_time_series(sample_frq, [(U_[::2], v_[::2], w_[::2]) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
np.testing.assert_allclose(uu[:-3], uu3[:-2], 0.1)
# One set of time series with U=4
U = u + 4
Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
# Another set of time series with U=8 i.e. only every second point to have
# same sample_frq. (nan added to have same length)
U = u + 4
Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
sample_frq = 2
with warnings.catch_warnings():
warnings.simplefilter("ignore")
k14, uu4, vv4, ww4, uw4 = logbin_spectra(*spectra_from_time_series(sample_frq, Uvw_lst))
np.testing.assert_allclose(uu[:-3], uu3[:-2], rtol=0.1)
if 0:
import matplotlib.pyplot as plt
plt.semilogx(k1, uu * k1)
plt.semilogx(k12, uu2 * k12)
plt.semilogx(k13, uu3 * k13)
plt.semilogx(k14, uu4 * k14)
plt.show()
def test_fit_mann_parameters_from_timeseries(self):
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
ae, L, G = fit_mann_parameters(fx, u, v, w)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
#import matplotlib.pyplot as plt
plt = None
U = u + 4
sample_frq = 2
ae, L, G = fit_mann_parameters_from_time_series(
sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)], plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
# One set of time series with U=4
U = u + 4
Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
# Another set of time series with U=8 i.e. only every second point to have
# same sample_frq. (nan added to have same length)
U = u + 4
Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
sample_frq = 2
ae, L, G = fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
def test_var2ae_U(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
dx = 2
for U in [1, 10, 100]:
# should be independent of U
dt = dx / U
T = 16384 / U
self.assertAlmostEqual(var2ae(variance=u.var(), L=10, G=3.3, U=U, T=T, sample_frq=1 / dt), .15, delta=.024)
def test_var2ae_T(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
dx = 2
U = 10
dt = dx / U
for i in np.arange(6):
# reshape to more and shorter series. Variance should decrease while ae should be roughly constant
n = 2**i
u_ = u.T.reshape((u.T.shape * np.array([n, 1 / n])).astype(int)).T
var = u_.var(0).mean()
ae = var2ae(variance=var, L=10, G=3.3, U=U, T=dx * u_.shape[0] / U, sample_frq=1 / dt)
self.assertAlmostEqual(ae, .15, delta=.027)
def test_var2ae_dt(self):
u = mann_turbulence.load(get_test_file("h2a16384_8_8_65536_32_32_0.15_40_4.0u.dat"), (16384, 8, 8))
dx = 4
U = 10
T = u.shape[0] * dx / U
for i in np.arange(9):
# average every neighbouring samples to decrease dt.
# Variance should decrease while ae should be roughly constant
n = 2**i
u_ = u.reshape(u.shape[0] // n, n, u.shape[1]).mean(1)
var = u_.var(0).mean()
ae = var2ae(variance=var, L=40, G=4, U=U, T=T, sample_frq=1 / (n * dx / U))
self.assertAlmostEqual(ae, .15, delta=.041 + i / 100)
def test_fit_ae2var(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
self.assertAlmostEqual(fit_ae(spatial_resolution=2, u=u, L=10, G=3.3), .15, delta=.02)
def test_ae2ti(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
dx = 2
U = 10
dt = dx / U
T = 16384 / U
ae23 = var2ae(variance=u.var(), L=10, G=3.3, U=U, T=T, sample_frq=1 / dt)
ti = u.std() / U
self.assertAlmostEqual(ae2ti(ae23, L=10, G=3.3, U=U, T=T, sample_frq=1 / dt), ti)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
wetb/wind/utils.py
View file @ 62f0f2e4
...@@ -38,14 +38,15 @@ def wsp_dir2uv(wsp, dir, dir_ref=None): ...@@ -38,14 +38,15 @@ def wsp_dir2uv(wsp, dir, dir_ref=None):
v = -np.sin(rad(dir)) * wsp[:] v = -np.sin(rad(dir)) * wsp[:]
return np.array([u, v]) return np.array([u, v])
def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None): def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
"""Convert horizontal wind speed and direction to u,v,w r"""Convert horizontal wind speed and direction to u,v,w
Parameters Parameters
---------- ----------
wsp : array_like wsp : array_like
- if wsp_horizontal is True: Horizontal wind speed, $\sqrt{u^2+v^2}\n - if wsp_horizontal is True: Horizontal wind speed, $\sqrt{u^2+v^2}$\n
- if wsp_horizontal is False: Wind speed, $\sqrt{u^2+v^2+w^2} - if wsp_horizontal is False: Wind speed, $\sqrt{u^2+v^2+w^2}$
dir : array_like dir : array_like
Wind direction Wind direction
tilt : array_like tilt : array_like
...@@ -66,6 +67,7 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None): ...@@ -66,6 +67,7 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
w : array_like w : array_like
v wind component v wind component
""" """
wsp, dir, tilt = wsp[:], dir[:], tilt[:] wsp, dir, tilt = wsp[:], dir[:], tilt[:]
if wsp_horizontal: if wsp_horizontal:
w = tand(tilt) * wsp w = tand(tilt) * wsp
...@@ -76,7 +78,6 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None): ...@@ -76,7 +78,6 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
return np.array([u, v, w]) return np.array([u, v, w])
def xyz2uvw(x, y, z, left_handed=True): def xyz2uvw(x, y, z, left_handed=True):
"""Convert sonic x,y,z measurements to u,v,w wind components """Convert sonic x,y,z measurements to u,v,w wind components
...@@ -108,7 +109,7 @@ def xyz2uvw(x, y, z, left_handed=True): ...@@ -108,7 +109,7 @@ def xyz2uvw(x, y, z, left_handed=True):
SV = cosd(theta) * y - sind(theta) * x SV = cosd(theta) * y - sind(theta) * x
# SUW = cosd(theta) * x + sind(theta) * y # SUW = cosd(theta) * x + sind(theta) * y
# #
# #% rotation around y of tilt # #% rotation around y of tilt
# tilt = deg(np.arctan2(np.mean(z), np.mean(SUW))) # tilt = deg(np.arctan2(np.mean(z), np.mean(SUW)))
# SU = SUW * cosd(tilt) + z * sind(tilt); # SU = SUW * cosd(tilt) + z * sind(tilt);
...@@ -120,8 +121,6 @@ def xyz2uvw(x, y, z, left_handed=True): ...@@ -120,8 +121,6 @@ def xyz2uvw(x, y, z, left_handed=True):
return np.array([SU, SV, SW]) return np.array([SU, SV, SW])
def abvrel2xyz_old(alpha, beta, vrel): def abvrel2xyz_old(alpha, beta, vrel):
"""Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
...@@ -143,27 +142,29 @@ def abvrel2xyz_old(alpha, beta, vrel): ...@@ -143,27 +142,29 @@ def abvrel2xyz_old(alpha, beta, vrel):
z : array_like z : array_like
Wind component in beta plane (positive for negative beta) Wind component in beta plane (positive for negative beta)
""" """
alpha = np.array(alpha, dtype=np.float) alpha = np.array(alpha, dtype=float)
beta = np.array(beta, dtype=np.float) beta = np.array(beta, dtype=float)
vrel = np.array(vrel, dtype=np.float) vrel = np.array(vrel, dtype=float)
sign_vsx = -((np.abs(beta) > np.pi / 2) * 2 - 1) # +1 for |beta| < 90, -1 for |beta|>90 sign_vsx = -((np.abs(beta) > np.pi / 2) * 2 - 1) # +1 for |beta| < 90, -1 for |beta|>90
sign_vsy = np.sign(alpha) #+ for alpha > 0 sign_vsy = np.sign(alpha) # + for alpha > 0
sign_vsz = -np.sign(beta) #- for beta>0 sign_vsz = -np.sign(beta) # - for beta>0
x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2)) x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
m = alpha != 0 m = alpha != 0
y = np.zeros_like(alpha) y = np.zeros_like(alpha)
y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 + 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2)) y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 +
1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
m = beta != 0 m = beta != 0
z = np.zeros_like(alpha) z = np.zeros_like(alpha)
z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2)) z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 +
1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
return x, y, z return x, y, z
def abvrel2xyz(alpha, beta, vrel): def abvrel2xyz(alpha, beta, vrel):
"""Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
...@@ -197,41 +198,42 @@ def abvrel2xyz(alpha, beta, vrel): ...@@ -197,41 +198,42 @@ def abvrel2xyz(alpha, beta, vrel):
z : array_like z : array_like
Wind component in beta plane (positive for negative beta) Wind component in beta plane (positive for negative beta)
""" """
alpha = np.array(alpha, dtype=np.float) alpha = np.array(alpha, dtype=float)
beta = np.array(beta, dtype=np.float) beta = np.array(beta, dtype=float)
vrel = np.array(vrel, dtype=np.float) vrel = np.array(vrel, dtype=float)
sign_vsx = ((np.abs(beta) > np.pi / 2) * 2 - 1) # -1 for |beta| < 90, +1 for |beta|>90 sign_vsx = ((np.abs(beta) > np.pi / 2) * 2 - 1) # -1 for |beta| < 90, +1 for |beta|>90
sign_vsy = -np.sign(alpha) #- for alpha > 0 sign_vsy = -np.sign(alpha) # - for alpha > 0
sign_vsz = np.sign(beta) # for beta>0 sign_vsz = np.sign(beta) # for beta>0
x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2)) x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
m = alpha != 0 m = alpha != 0
y = np.zeros_like(alpha) y = np.zeros_like(alpha)
y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 + 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2)) y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 +
1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
m = beta != 0 m = beta != 0
z = np.zeros_like(alpha) z = np.zeros_like(alpha)
z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2)) z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 +
1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
return np.array([x, y, z]).T return np.array([x, y, z]).T
def detrend_uvw(u, v=None, w=None): def detrend_uvw(u, v=None, w=None):
# def _detrend(wsp): # def _detrend(wsp):
# if wsp is None: # if wsp is None:
# return None # return None
# dwsp = np.atleast_2d(wsp.copy().T).T # dwsp = np.atleast_2d(wsp.copy().T).T
# t = np.arange(dwsp.shape[0]) # t = np.arange(dwsp.shape[0])
# A = np.vstack([t, np.ones(len(t))]).T # A = np.vstack([t, np.ones(len(t))]).T
# for i in range(dwsp.shape[1]): # for i in range(dwsp.shape[1]):
# trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0] # trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0]
# dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend # dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
# return dwsp.reshape(wsp.shape) # return dwsp.reshape(wsp.shape)
def _detrend(y): def _detrend(y):
if y is None: if y is None:
return None return None
return detrend(y) return detrend(y)
return [_detrend(uvw) for uvw in [u, v, w]] return [_detrend(uvw) for uvw in [u, v, w]]
\ No newline at end of file
wetb/wind/weibull.py
View file @ 62f0f2e4
...@@ -103,7 +103,7 @@ def fit(wsp): ...@@ -103,7 +103,7 @@ def fit(wsp):
if __name__ == "__main__": if __name__ == "__main__":
from wetb.wind import weibull from wetb.wind import weibull
from pylab import hist, show, plot from pylab import hist, show, plot
hist(weibull.random(10, 2, 10000), 20, normed=True) hist(weibull.random(10, 2, 10000), 20, density=True)
wsp = np.arange(0, 20, .5) wsp = np.arange(0, 20, .5)
plot(wsp, weibull.pdf(10, 2)(wsp)) plot(wsp, weibull.pdf(10, 2)(wsp))
show() show()