geometry + tests

README.md

@@ -17,4 +17,5 @@ General Time Series Data Format, a binary hdf5 data format for storing time seri
 
 ### [functions](wetb/functions)
 Other functions
+- [geometry](wetb/functions/geometry.py): Different kind of geometry conversion functions
 - [process_exec](wetb/functions/process_exec.py): Run system command in subprocess

wetb/functions/geometry.py

+import numpy as np
+
+def rad(deg):
+    return deg * np.pi / 180
+
+def deg(rad):
+    return rad / np.pi * 180
+
+def sind(dir_deg):
+    return np.sin(rad(dir_deg))
+
+def cosd(dir_deg):
+    return np.cos(rad(dir_deg))
+
+def tand(dir_deg):
+    return np.tan(rad(dir_deg))
+
+def mean_deg(dir):
+    """Mean of angles in degrees
+
+    Parameters
+    ----------
+    dir : array_like
+        Angles in degrees
+
+    Returns
+    -------
+    mean_deg : float
+        Mean angle
+    """
+    return deg(np.arctan2(np.mean(sind(dir[:])), np.mean(cosd(dir[:]))))
+
+def std_deg(dir):
+    """Standard deviation of angles in degrees
+
+    Parameters
+    ----------
+    dir : array_like
+        Angles in degrees
+
+    Returns
+    -------
+    std_deg : float
+        standard deviation
+    """
+    return deg(np.sqrt(1 - (np.mean(sind(dir)) ** 2 + np.mean(cosd(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(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
+    beta : array_like
+        Pitot tube side slip angle
+    vrel : array_like
+        Pitot tube relative velocity
+
+    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 positive 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 * 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 * 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

wetb/functions/tests/test_geometry.py

+'''
+Created on 15/01/2014
+
+@author: MMPE
+'''
+import unittest
+
+import wetb.gtsdf
+import numpy as np
+from wetb.functions.geometry import rad, deg, mean_deg, sind, cosd, std_deg, xyz2uvw, \
+    wsp_dir2uv, wsp_dir_tilt2uvw, tand
+
+
+class Test(unittest.TestCase):
+
+
+    def test_rad(self):
+        self.assertEqual(rad(45), np.pi / 4)
+        self.assertEqual(rad(135), np.pi * 3 / 4)
+
+
+    def test_deg(self):
+        self.assertEqual(45, deg(np.pi / 4))
+        self.assertEqual(135, deg(np.pi * 3 / 4))
+
+    def test_rad_deg(self):
+        for i in [15, 0.5, 355, 400]:
+            self.assertEqual(i, deg(rad(i)), i)
+
+    def test_sind(self):
+        self.assertAlmostEqual(sind(30), .5)
+
+    def test_cosd(self):
+        self.assertAlmostEqual(cosd(60), .5)
+
+    def test_tand(self):
+        self.assertAlmostEqual(tand(30), 0.5773, 3)
+
+
+    def test_mean_deg(self):
+        self.assertEqual(mean_deg(np.array([0, 90])), 45)
+        self.assertAlmostEqual(mean_deg(np.array([350, 10])), 0)
+
+    def test_std_deg(self):
+        self.assertEqual(std_deg(np.array([0, 0, 0])), 0)
+        self.assertAlmostEqual(std_deg(np.array([0, 90, 180, 270])), 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("test_files/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("test_files/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("test_files/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))
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.test_rad']
+    unittest.main()