From 4cb24d2b4667888dfc4d4747d42876e20d8ae6bc Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 21 Dec 2016 11:08:50 +0100
Subject: [PATCH] wsp and pitot mapping moved to wind.dir_mapping.py
---
.gitignore | 2 +
wetb/utils/geometry.py | 204 +---------------------------
wetb/utils/tests/test_geometry.py | 68 +---------
wetb/wind/dir_mapping.py | 214 ++++++++++++++++++++++++++++++
4 files changed, 219 insertions(+), 269 deletions(-)
create mode 100644 wetb/wind/dir_mapping.py
diff --git a/.gitignore b/.gitignore
index 5c2598f..1e37b33 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,3 +20,5 @@ wetb/hawc2/ascii2bin/tests/test_files/Hawc2ascii_bin.dat
wetb/prepost/tests/data/demo_dlc/remote*
/wetb/fatigue_tools/rainflowcounting/compile.py
+/docs/api
+/htmlcov
diff --git a/wetb/utils/geometry.py b/wetb/utils/geometry.py
index 3a7761c..211a811 100644
--- a/wetb/utils/geometry.py
+++ b/wetb/utils/geometry.py
@@ -87,207 +87,5 @@ def std_rad(dir):
"""
return np.sqrt(1 - (np.nanmean(np.sin(dir)) ** 2 + np.nanmean(np.cos(dir)) ** 2))
-def wsp_dir2uv(wsp, dir, dir_ref=None):
- """Convert horizontal wind speed and direction to u,v
-
- Parameters
- ----------
- wsp : array_like
- Horizontal wind speed
- dir : array_like
- Wind direction
- dir_ref : int or float, optional
- Reference direction\n
- If None, default, the mean direction is used as reference
-
- Returns
- -------
- u : array_like
- u wind component
- v : array_like
- v wind component
- """
- if dir_ref is None:
- dir = dir[:] - mean_deg(dir[:])
- else:
- dir = dir[:] - dir_ref
- u = np.cos(rad(dir)) * wsp[:]
- v = -np.sin(rad(dir)) * wsp[:]
- return np.array([u, v])
-
-def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
- """Convert horizontal wind speed and direction to u,v,w
-
- Parameters
- ----------
- wsp : array_like
- - 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}
- dir : array_like
- Wind direction
- tilt : array_like
- Wind tilt
- wsp_horizontal : bool
- See wsp
- dir_ref : int or float, optional
- Reference direction\n
- If None, default, the mean direction is used as reference
-
-
- Returns
- -------
- u : array_like
- u wind component
- v : array_like
- v wind component
- w : array_like
- v wind component
- """
- wsp, dir, tilt = wsp[:], dir[:], tilt[:]
- if wsp_horizontal:
- w = tand(tilt) * wsp
- u, v = wsp_dir2uv(wsp, dir, dir_ref)
- else:
- w = sind(tilt) * wsp
- u, v = wsp_dir2uv(np.sqrt(wsp ** 2 - w ** 2), dir, dir_ref)
- return np.array([u, v, w])
-
-
-
-def xyz2uvw(x, y, z, left_handed=True):
- """Convert sonic x,y,z measurements to u,v,w wind components
-
- Parameters
- ----------
- x : array_like
- Sonic x component
- y : array_like
- Sonic x component
- z : array_like
- Sonic x component
- left_handed : boolean
- if true (default), xyz are defined in left handed coodinate system (default for some sonics)
- if false, xyz are defined in normal right handed coordinate system
-
- Returns
- -------
- u : array_like
- u wind component
- v : array_like
- v wind component
- w : array_like
- w wind component
- """
- x, y, z = map(np.array, [x, y, z])
- if left_handed:
- y *= -1
- theta = deg(np.arctan2(np.mean(y), np.mean(x)))
- SV = cosd(theta) * y - sind(theta) * x
- SUW = cosd(theta) * x + sind(theta) * y
-
- #% rotation around y of tilt
- tilt = deg(np.arctan2(np.mean(z), np.mean(SUW)))
- SU = SUW * cosd(tilt) + z * sind(tilt);
- SW = z * cosd(tilt) - SUW * sind(tilt);
-
- return np.array([SU, SV, SW])
-
def rpm2rads(rpm):
- return rpm * 2 * np.pi / 60
-
-
-def abvrel2xyz_old(alpha, beta, vrel):
- """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
-
- Parameters
- ----------
- alpha : array_like
- Pitot tube angle of attack [rad]. Zero: Parallel to pitot tube. Positive: Flow from wind side (pressure side)
- beta : array_like
- Pitot tube side slip angle [rad]. Zero: Parallel to pitot tube. Positive: Flow from root side
- vrel : array_like
- Pitot tube relative velocity. Positive: flow towards pitot tube
-
- Returns
- -------
- x : array_like
- Wind component towards pitot tube (positive for postive vrel and -90<beta<90)
- y : array_like
- Wind component in alpha plane (positive for positive alpha)
- z : array_like
- Wind component in beta plane (positive for negative beta)
- """
- alpha = np.array(alpha, dtype=np.float)
- beta = np.array(beta, dtype=np.float)
- vrel = np.array(vrel, dtype=np.float)
-
- 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_vsz = -np.sign(beta) #- for beta>0
-
-
- x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
-
- m = alpha != 0
- 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))
-
- m = beta != 0
- 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))
-
- return x, y, z
-
-def abvrel2xyz(alpha, beta, vrel):
- """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
-
- x : parallel to pitot tube, direction pitot tube root to tip, i.e. normal flow gives negative x\n
- y : component in alpha plane
- z : component in beta plane
-
- For typical usage where pitot tube is mounted on leading edge:\n
- x: Opposite rotational direction\n
- y: Direction of mean wind\n
- z: From blade root to tip\n
-
-
- Parameters
- ----------
- alpha : array_like
- Pitot tube angle of attack [rad]. Zero for flow towards pitot tube. Positive around z-axis. I.e.
- negative alpha (normal flow) gives positive y component
- beta : array_like
- Pitot tube side slip angle [rad]. Zero for flow towards pitot tube. Positive around y-axis. I.e.
- Positive beta (normal flow due to expansion and position in front of blade) gives positive z
- vrel : array_like
- Pitot tube relative velocity. Positive: flow towards pitot tube
-
- Returns
- -------
- x : array_like
- Wind component away from pitot tube (positive for postive vrel and -90<beta<90)
- y : array_like
- Wind component in alpha plane (positive for positive alpha)
- z : array_like
- Wind component in beta plane (positive for negative beta)
- """
- alpha = np.array(alpha, dtype=np.float)
- beta = np.array(beta, dtype=np.float)
- vrel = np.array(vrel, dtype=np.float)
-
- 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_vsz = np.sign(beta) # for beta>0
-
-
- x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
-
- m = alpha != 0
- 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))
-
- m = beta != 0
- 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))
-
- return np.array([x, y, z]).T
+ return rpm * 2 * np.pi / 60
\ No newline at end of file
diff --git a/wetb/utils/tests/test_geometry.py b/wetb/utils/tests/test_geometry.py
index 1ae8283..a624ce7 100644
--- a/wetb/utils/tests/test_geometry.py
+++ b/wetb/utils/tests/test_geometry.py
@@ -13,8 +13,7 @@ import unittest
import wetb.gtsdf
import numpy as np
-from wetb.utils.geometry import rad, deg, mean_deg, sind, cosd, std_deg, xyz2uvw, \
- wsp_dir2uv, wsp_dir_tilt2uvw, tand, abvrel2xyz
+from wetb.utils.geometry import rad, deg, mean_deg, sind, cosd, std_deg, tand
import os
@@ -68,70 +67,7 @@ class TestGeometry(unittest.TestCase):
def test_std_deg_nan(self):
self.assertAlmostEqual(std_deg(np.array([0, 90, 180, 270, np.nan])), 57.296, 2)
- def test_wspdir2uv(self):
- u, v = wsp_dir2uv(np.array([1, 1, 1]), np.array([30, 0, 330]))
- np.testing.assert_array_almost_equal(u, [0.8660, 1, 0.8660], 3)
- np.testing.assert_array_almost_equal(v, [-0.5, 0, 0.5], 3)
-
- def test_wspdir2uv_dir_ref(self):
- u, v = wsp_dir2uv(np.array([1, 1, 1]), np.array([30, 0, 330]), 30)
- np.testing.assert_array_almost_equal(u, [1, 0.8660, .5], 3)
- np.testing.assert_array_almost_equal(v, [0, 0.5, .8660], 3)
-
- def test_xyz2uvw(self):
- u, v, w = xyz2uvw([1, 1, 0], [0, 1, 1], 0, left_handed=False)
- np.testing.assert_almost_equal(u, [np.sqrt(1 / 2), np.sqrt(2), np.sqrt(1 / 2)])
- np.testing.assert_almost_equal(v, [-np.sqrt(1 / 2), 0, np.sqrt(1 / 2)])
-
-
- u, v, w = xyz2uvw([1, 1, 0], [0, 1, 1], 0, left_handed=True)
- np.testing.assert_almost_equal(u, [np.sqrt(1 / 2), np.sqrt(2), np.sqrt(1 / 2)])
- np.testing.assert_almost_equal(v, [np.sqrt(1 / 2), 0, -np.sqrt(1 / 2)])
-
- u, v, w = xyz2uvw(np.array([-1, -1, -1]), np.array([-0.5, 0, .5]), np.array([0, 0, 0]), left_handed=False)
- np.testing.assert_array_almost_equal(u, np.array([1, 1, 1]))
- np.testing.assert_array_almost_equal(v, np.array([.5, 0, -.5]))
- np.testing.assert_array_almost_equal(w, np.array([0, 0, 0]))
-
- u, v, w = xyz2uvw(np.array([.5, cosd(30), 1]), np.array([sind(60), sind(30), 0]), np.array([0, 0, 0]), left_handed=False)
- np.testing.assert_array_almost_equal(u, np.array([sind(60), 1, sind(60)]))
- np.testing.assert_array_almost_equal(v, np.array([.5, 0, -.5]))
- np.testing.assert_array_almost_equal(w, np.array([0, 0, 0]))
-
- u, v, w = xyz2uvw(np.array([.5, cosd(30), 1]), np.array([0, 0, 0]), np.array([sind(60), sind(30), 0]), left_handed=False)
- np.testing.assert_array_almost_equal(u, np.array([sind(60), 1, sind(60)]))
- np.testing.assert_array_almost_equal(v, np.array([0, 0, 0]))
- np.testing.assert_array_almost_equal(w, np.array([.5, 0, -.5]))
-
-
- def test_wspdir2uv2(self):
- time, data, info = wetb.gtsdf.load(self.tfp + "SonicDataset.hdf5")
- stat, x, y, z, temp, wsp, dir, tilt = data[2:3].T #xyz is left handed
- np.testing.assert_array_almost_equal(xyz2uvw(*wsp_dir2uv(wsp, dir), z=0), xyz2uvw(x, y, 0))
-
- def test_wspdirtil2uvw(self):
- time, data, info = wetb.gtsdf.load(self.tfp + "SonicDataset.hdf5")
- stat, x, y, z, temp, wsp, dir, tilt = data[3:6].T #xyz is left handed
- wsp = np.sqrt(wsp ** 2 + z ** 2)
- np.testing.assert_array_almost_equal(xyz2uvw(*wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=False), left_handed=False), xyz2uvw(x, y, z))
-
- def test_wspdirtil2uvw_horizontal_wsp(self):
- time, data, info = wetb.gtsdf.load(self.tfp + "SonicDataset.hdf5")
- stat, x, y, z, temp, wsp, dir, tilt = data[:].T #xyz is left handed
- np.testing.assert_array_almost_equal(xyz2uvw(*wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=True), left_handed=False), xyz2uvw(x, y, z))
-
- np.testing.assert_array_almost_equal(wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=True, dir_ref=180), np.array([x, -y, z]), 5)
- np.testing.assert_array_almost_equal(xyz2uvw(*wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal=True), left_handed=False), xyz2uvw(x, y, z))
-
- def test_abvrel2xyz(self):
- abvrel = np.array([(0., 0, 1), (30, 0, 1), (-30, 0, 1), (0, 30, 1), (0, -30, 1), (30, 30, 1), (30, 30, 2)])
- abvrel[:, :2] = rad(abvrel[:, :2])
- for (x, y, z), (a, b, vrel) in zip(abvrel2xyz(*abvrel.T), abvrel):
- #print (deg(a), deg(b), vrel, x, y, z)
- self.assertAlmostEqual(a, np.arctan(y / x))
- self.assertAlmostEqual(b, np.arctan(-z / x))
- self.assertAlmostEqual(vrel, np.sqrt(x ** 2 + y ** 2 + z ** 2))
-
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.test_rad']
diff --git a/wetb/wind/dir_mapping.py b/wetb/wind/dir_mapping.py
new file mode 100644
index 0000000..4da160f
--- /dev/null
+++ b/wetb/wind/dir_mapping.py
@@ -0,0 +1,214 @@
+'''
+Created on 19. dec. 2016
+
+@author: mmpe
+'''
+
+from wetb.utils.geometry import mean_deg, rad, tand, sind, deg, cosd
+
+import numpy as np
+
+
+def wsp_dir2uv(wsp, dir, dir_ref=None):
+ """Convert horizontal wind speed and direction to u,v
+
+ Parameters
+ ----------
+ wsp : array_like
+ Horizontal wind speed
+ dir : array_like
+ Wind direction
+ dir_ref : int or float, optional
+ Reference direction\n
+ If None, default, the mean direction is used as reference
+
+ Returns
+ -------
+ u : array_like
+ u wind component
+ v : array_like
+ v wind component
+ """
+ if dir_ref is None:
+ dir = dir[:] - mean_deg(dir[:])
+ else:
+ dir = dir[:] - dir_ref
+ u = np.cos(rad(dir)) * wsp[:]
+ v = -np.sin(rad(dir)) * wsp[:]
+ return np.array([u, v])
+
+def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
+ """Convert horizontal wind speed and direction to u,v,w
+
+ Parameters
+ ----------
+ wsp : array_like
+ - 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}
+ dir : array_like
+ Wind direction
+ tilt : array_like
+ Wind tilt
+ wsp_horizontal : bool
+ See wsp
+ dir_ref : int or float, optional
+ Reference direction\n
+ If None, default, the mean direction is used as reference
+
+
+ Returns
+ -------
+ u : array_like
+ u wind component
+ v : array_like
+ v wind component
+ w : array_like
+ v wind component
+ """
+ wsp, dir, tilt = wsp[:], dir[:], tilt[:]
+ if wsp_horizontal:
+ w = tand(tilt) * wsp
+ u, v = wsp_dir2uv(wsp, dir, dir_ref)
+ else:
+ w = sind(tilt) * wsp
+ u, v = wsp_dir2uv(np.sqrt(wsp ** 2 - w ** 2), dir, dir_ref)
+ return np.array([u, v, w])
+
+
+
+def xyz2uvw(x, y, z, left_handed=True):
+ """Convert sonic x,y,z measurements to u,v,w wind components
+
+ Parameters
+ ----------
+ x : array_like
+ Sonic x component
+ y : array_like
+ Sonic x component
+ z : array_like
+ Sonic x component
+ left_handed : boolean
+ if true (default), xyz are defined in left handed coodinate system (default for some sonics)
+ if false, xyz are defined in normal right handed coordinate system
+
+ Returns
+ -------
+ u : array_like
+ u wind component
+ v : array_like
+ v wind component
+ w : array_like
+ w wind component
+ """
+ x, y, z = map(np.array, [x, y, z])
+ if left_handed:
+ y *= -1
+ theta = deg(np.arctan2(np.mean(y), np.mean(x)))
+ SV = cosd(theta) * y - sind(theta) * x
+ SUW = cosd(theta) * x + sind(theta) * y
+
+ #% rotation around y of tilt
+ tilt = deg(np.arctan2(np.mean(z), np.mean(SUW)))
+ SU = SUW * cosd(tilt) + z * sind(tilt);
+ SW = z * cosd(tilt) - SUW * sind(tilt);
+
+ return np.array([SU, SV, SW])
+
+
+
+
+def abvrel2xyz_old(alpha, beta, vrel):
+ """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
+
+ Parameters
+ ----------
+ alpha : array_like
+ Pitot tube angle of attack [rad]. Zero: Parallel to pitot tube. Positive: Flow from wind side (pressure side)
+ beta : array_like
+ Pitot tube side slip angle [rad]. Zero: Parallel to pitot tube. Positive: Flow from root side
+ vrel : array_like
+ Pitot tube relative velocity. Positive: flow towards pitot tube
+
+ Returns
+ -------
+ x : array_like
+ Wind component towards pitot tube (positive for postive vrel and -90<beta<90)
+ y : array_like
+ Wind component in alpha plane (positive for positive alpha)
+ z : array_like
+ Wind component in beta plane (positive for negative beta)
+ """
+ alpha = np.array(alpha, dtype=np.float)
+ beta = np.array(beta, dtype=np.float)
+ vrel = np.array(vrel, dtype=np.float)
+
+ 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_vsz = -np.sign(beta) #- for beta>0
+
+
+ x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
+
+ m = alpha != 0
+ 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))
+
+ m = beta != 0
+ 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))
+
+ return x, y, z
+
+def abvrel2xyz(alpha, beta, vrel):
+ """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
+
+ x : parallel to pitot tube, direction pitot tube root to tip, i.e. normal flow gives negative x\n
+ y : component in alpha plane
+ z : component in beta plane
+
+ For typical usage where pitot tube is mounted on leading edge:\n
+ x: Opposite rotational direction\n
+ y: Direction of mean wind\n
+ z: From blade root to tip\n
+
+
+ Parameters
+ ----------
+ alpha : array_like
+ Pitot tube angle of attack [rad]. Zero for flow towards pitot tube. Positive around z-axis. I.e.
+ negative alpha (normal flow) gives positive y component
+ beta : array_like
+ Pitot tube side slip angle [rad]. Zero for flow towards pitot tube. Positive around y-axis. I.e.
+ Positive beta (normal flow due to expansion and position in front of blade) gives positive z
+ vrel : array_like
+ Pitot tube relative velocity. Positive: flow towards pitot tube
+
+ Returns
+ -------
+ x : array_like
+ Wind component away from pitot tube (positive for postive vrel and -90<beta<90)
+ y : array_like
+ Wind component in alpha plane (positive for positive alpha)
+ z : array_like
+ Wind component in beta plane (positive for negative beta)
+ """
+ alpha = np.array(alpha, dtype=np.float)
+ beta = np.array(beta, dtype=np.float)
+ vrel = np.array(vrel, dtype=np.float)
+
+ 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_vsz = np.sign(beta) # for beta>0
+
+
+ x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
+
+ m = alpha != 0
+ 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))
+
+ m = beta != 0
+ 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))
+
+ return np.array([x, y, z]).T
--
GitLab