New stuff

wetb/utils/euler.py

+'''
+Created on 15/01/2014
+
+@author: MMPE
+'''
+
+import numpy as np
+from wetb.utils.geometry import deg
+
+
+def Ax(angle):
+    cos = np.cos(angle)
+    sin = np.sin(angle)
+    return np.array([[1, 0, 0],
+                   [0, cos, -sin],
+                   [0, sin, cos]])
+
+def Ay(angle):
+    cos = np.cos(angle)
+    sin = np.sin(angle)
+    return np.array([[cos, 0, sin],
+                   [0, 1, 0       ],
+                   [-sin, 0, cos  ]])
+
+def Az(angle):
+    cos = np.cos(angle)
+    sin = np.sin(angle)
+    return np.array([[cos, -sin, 0],
+                   [sin, cos, 0],
+                   [0, 0, 1]])
+
+
+
+
+def euler2A(euler_param):
+    assert len(euler_param) == 4
+    e = euler_param
+    return np.array([[e[0] ** 2 + e[1] ** 2 - e[2] ** 2 - e[3] ** 2, 2 * (e[1] * e[2] + e[0] * e[3]) , 2 * (e[1] * e[3] - e[0] * e[2]) ],
+                         [2 * (e[1] * e[2] - e[0] * e[3]), e[0] ** 2 - e[1] ** 2 + e[2] ** 2 - e[3] ** 2, 2 * (e[2] * e[3] + e[0] * e[1]) ],
+                         [2 * (e[1] * e[3] + e[0] * e[2]), 2 * (e[2] * e[3] - e[0] * e[1]), e[0] ** 2 - e[1] ** 2 - e[2] ** 2 + e[3] ** 2]]).T
+
+def A2euler(A):
+    # method from http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
+    sqrt = np.sqrt
+    (m00, m01, m02), (m10, m11, m12), (m20, m21, m22) = A
+    tr = m00 + m11 + m22
+
+    if (tr > 0):
+        S = sqrt(tr + 1.0) * 2  # // S=4*qw
+        qw = 0.25 * S
+        qx = (m21 - m12) / S
+        qy = (m02 - m20) / S
+        qz = (m10 - m01) / S
+    elif ((m00 > m11) and (m00 > m22)):
+        S = sqrt(1.0 + m00 - m11 - m22) * 2  # // S=4*qx
+        qw = (m21 - m12) / S
+        qx = 0.25 * S
+        qy = (m01 + m10) / S
+        qz = (m02 + m20) / S
+    elif (m11 > m22):
+        S = sqrt(1.0 + m11 - m00 - m22) * 2  # // S=4*qy
+        qw = (m02 - m20) / S
+        qx = (m01 + m10) / S
+        qy = 0.25 * S
+        qz = (m12 + m21) / S
+    else:
+        S = sqrt(1.0 + m22 - m00 - m11) * 2  # // S=4*qz
+        qw = (m10 - m01) / S
+        qx = (m02 + m20) / S
+        qy = (m12 + m21) / S
+        qz = 0.25 * S
+    return np.array([qw, qx, qy, qz])
+
+
+
+#def A2xyz(A):
+#    if abs(A[2, 0]) != 1:
+#        y = -np.arcsin(A[2, 0])
+#        x = np.arctan2(A[2, 1] / np.cos(y), A[2, 2] / np.cos(y))
+#        z = np.arctan2(A[1, 0] / np.cos(y), A[0, 0] / np.cos(y))
+#    else:
+#        z = 0
+#        if A[2, 0] == -1:
+#            y = np.pi / 2
+#            x = z + np.arctan(A[0, 1], A[0, 2])
+#        else:
+#            y = -np.pi / 2
+#            x = -z + np.arctan(-A[0, 1], -A[0, 2])
+#    return np.array([x, y, z])
+#
+#def zxz2euler(z1, x, z2):
+#    return np.array([np.cos(.5 * (z1 + z2)) * np.cos(.5 * x),
+#                     np.cos(.5 * (z1 - z2)) * np.sin(.5 * x),
+#                     np.sin(.5 * (z1 - z2)) * np.sin(.5 * x),
+#                     np.sin(.5 * (z1 + z2)) * np.cos(.5 * x)])
+#
+#def xyz2A(x, y, z):
+#    return np.dot(Ax(x), np.dot(Ay(y), Az(z)))
+
+#def euler2xyz(euler):
+#    return A2xyz(euler2A(euler))
+
+#def A2euler(A):
+#    return xyz2euler(*A2xyz(A))
+
+def euler2angle(euler):
+    if euler[0] > 1:
+        euler[0] = 1
+    if euler[0] < -1:
+        euler[0] = -1
+
+    return np.arccos(euler[0]) * 2
+
+def euler2gl(euler):
+    return np.r_[deg(euler2angle(euler)), euler[1:]]

wetb/utils/rotation.py

+'''
+Created on 21/07/2014
+
+@author: mmpe
+'''
+import numpy as np
+
+
+def transformation_matrix(angles, xyz):
+    """Create Transformation matrix(es)
+    !!!Note that the columns of the returned matrix(es) is original(unrotate) xyz-axes in rotated coordinates\n
+    !!!Multiplying a this matrix by a vector rotates the vector -angle radians (in right handed terminology).
+
+    Parameters
+    ----------
+    angles : int, float or array_like
+        Angle(s) for rotation matrix(es). \n
+        One rotation matrix will be returned for each angle
+    xyz : {0,1,2}
+        - 0: rotation around x(first) axis\n
+        - 1: rotation around y(second) axis\n
+        - 2: rotation around z(third) axis\n
+
+    Returns
+    -------
+    transformation_matrix : array_like, shape = (no_angles, 3,3)
+        Rotation matrix(es)
+    """
+    indexes = [(0, 4, 8, 5, 7), (4, 0, 8, 6, 2), (8, 0, 4, 1, 3)]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
+    if isinstance(angles, (int, float)):
+        n = 1
+    else:
+        n = len(angles)
+    m = np.zeros(n * 9, np.float)
+    cosx = np.cos(angles)
+    sinx = np.sin(angles)
+    m[indexes[xyz][0]::9] = 1
+    m[indexes[xyz][1]::9] = cosx
+    m[indexes[xyz][2]::9] = cosx
+    m[indexes[xyz][3]::9] = sinx
+    m[indexes[xyz][4]::9] = -sinx
+    return m.reshape(n, 3, 3)
+
+def rotmat(angles, xyz):
+    """Create rotation matrix(es)
+    !!!Note that the columns of the returned matrix(es) is rotated xyz-axes in original(unrotated) coordinates\n
+    Multiplying a this matrix by a vector rotates the vector +angle radians (in right handed terminology).
+    Parameters
+    ----------
+    angles : int, float or array_like
+        Angle(s) for rotation matrix(es). \n
+        Note that
+        One rotation matrix will be returned for each angle
+    xyz : {0,1,2}
+        - 0: rotation around x(first) axis\n
+        - 1: rotation around y(second) axis\n
+        - 2: rotation around z(third) axis\n
+
+    Returns
+    -------
+    transformation_matrix : array_like, shape = (no_angles, 3,3)
+        Rotation matrix(es)
+    """
+    indexes = [(0, 4, 8, 5, 7), (4, 0, 8, 6, 2), (8, 0, 4, 1, 3)]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
+    if isinstance(angles, (int, float)):
+        n = 1
+    else:
+        n = len(angles)
+    m = np.zeros(n * 9, np.float)
+    cosx = np.cos(angles)
+    sinx = np.sin(angles)
+    m[indexes[xyz][0]::9] = 1
+    m[indexes[xyz][1]::9] = cosx
+    m[indexes[xyz][2]::9] = cosx
+    m[indexes[xyz][3]::9] = -sinx
+    m[indexes[xyz][4]::9] = sinx
+    return m.reshape(n, 3, 3)
+
+def mdot(m1, m2):
+    """Multiplication of matrix pairs
+
+    Parameters
+    ----------
+    m1 : array_like shape = (i, n, m)
+        First matrix set, i.e. i matrixes, size n x m\n
+        i must be equal to j or 1
+    m2 : array_like, shape = (j,m,p)
+        Second matrix set, i.e. j matrixes, size m x p\n
+        j must be equal to i or 1
+
+    Returns
+    -------
+    mdot : array_like, shape(max(i,j),n,p)
+        Matrix products
+    """
+    if m1.shape[0] == 1 and m2.shape[0] == 1:
+        return np.array([np.dot(m1[0], m2[0])])
+    elif m1.shape[0] > 1 and m2.shape[0] == 1:
+        mprod = np.empty_like(m1)
+        for i in range(m1.shape[0]):
+            mprod[i, :] = np.dot(m1[i, :], m2[0])
+    elif m1.shape[0] == 1 and m2.shape[0] > 1:
+        mprod = np.empty_like(m2)
+        for i in range(m2.shape[0]):
+            mprod[i, :] = np.dot(m1[0], m2[i, :])
+    elif m1.shape[0] > 1 and m2.shape[0] > 1 and m1.shape[0] == m2.shape[0]:
+        mprod = np.empty_like(m1)
+        for i in range(m1.shape[0]):
+            mprod[i, :] = np.dot(m1[i, :], m2[i, :])
+    else:
+        raise Exception("m1 and m2 must have same first dimension or one of the matrices must have first dimension 1")
+    return mprod
+
+
+def dots(rotmats, v):
+    """Rotate vector(s) v by rotation matrix(es)
+
+    Parameters
+    ----------
+    rotmats : array_like, shape=(i,3,3) or list of array_like, shape=(i,3,3)
+        if list of array_like, rotation matrixes are recursively reduced by multiplication\n
+        i must be 1 or equal to j
+    vector : array_like, shape=(j,3)
+        vectors to rotate
+
+    Returns
+    -------
+    dots : array_like, shape=(j,3)
+    """
+    if isinstance(rotmats, list):
+        if len(rotmats) > 1:
+            rotmats_red = [mdot(rotmats[0], rotmats[1])]
+            rotmats_red.extend(rotmats[2:])
+            return dots(rotmats_red, v)
+        else:
+            m = rotmats[0]
+    else:
+        m = rotmats
+    if len(v.shape) == 1:
+        v = np.array([v]).T
+
+    if m.shape[0] == 1:
+        return np.dot(m[0], v)
+    else:
+        if m.shape[0] != v.shape[1]:
+            raise Exception("m must have same dimension as v has number of columns")
+        res = np.zeros_like(v, dtype=np.float)
+        for i in range(v.shape[1]):
+            res[:, i] = np.dot(m[i], v[:, i])
+        return res
+
+
+def rotate(rotmats, v):
+    global x, y, z
+    v = np.array(v)
+    if isinstance(rotmats, (list, tuple)):
+        if len(rotmats) > 1:
+            rotmats_red = [mdot(rotmats[0], rotmats[1])]
+            rotmats_red.extend(rotmats[2:])
+            return rotate(rotmats_red, v)
+        else:
+            rotmats = rotmats[0]
+    assert rotmats.shape[-1] == rotmats.shape[-2] == v.shape[-1] == 3, "rotmats is %d x %x and v is i x %d, but must be 3" % (rotmats.shape[-1], rotmats.shape[-2], v.shape[-1])
+
+#    if isinstance(v, tuple):
+#        v = np.array([v])
+    n = (1, v.shape[0])[len(v.shape) == 2]
+
+    if rotmats.shape[0] == 1:
+        return np.dot(rotmats[0], v.T).T
+    else:
+        if v.shape[0] != rotmats.shape[0]:
+            raise Exception("V and rotmats must have same first dimension")
+        v = v.T
+        v_rot = np.zeros_like(v, dtype=np.float)
+        for i in range(n):
+            v_rot[:, i] = np.dot(rotmats[i], v[:, i])
+        return v_rot.T
+
+def rotate_x(v, angle):
+    """Rotate vector(s) around x axis
+
+    Parameters
+    ---------
+    v : array_like, shape=(3,) or shape=(1-N,3)
+        Vector(s) to rotate
+    angle : int, float or array_like shape=(1,) or shape=(N,)
+        Angle(s) [rad] to rotate
+
+    Returns
+    -------
+    y : array_like
+        Rotated vector(s)
+    """
+    return _rotate(v, angle, lambda x, y, z, cos, sin : [x, cos * y - sin * z, sin * y + cos * z])
+
+def rotate_y(v, angle):
+    """Rotate vector(s) around y axis
+
+    Parameters
+    ---------
+    v : array_like, shape=(3,) or shape=(1-N,3)
+        Vector(s) to rotate
+    angle : int, float or array_like shape=(1,) or shape=(N,)
+        Angle(s) [rad] to rotate
+
+    Returns
+    -------
+    y : array_like
+        Rotated vector(s)
+    """
+    return _rotate(v, angle, lambda x, y, z, cos, sin : [cos * x + sin * z, y, -sin * x + cos * z])
+
+def rotate_z(v, angle):
+    """Rotate vector(s) around z axis
+
+    Parameters
+    ---------
+    v : array_like, shape=(3,) or shape=(1-N,3)
+        Vector(s) to rotate
+    angle : int, float or array_like shape=(1,) or shape=(N,)
+        Angle(s) [rad] to rotate
+
+    Returns
+    -------
+    y : array_like
+        Rotated vector(s)
+    """
+    return _rotate(v, angle, lambda x, y, z, cos, sin : [cos * x - sin * y, sin * x + cos * y, z])
+
+def _rotate(v, angle, rotfunc):
+    angle = np.atleast_1d(angle)
+    cos, sin = np.cos(angle), np.sin(angle)
+    v = np.array(v)
+    if len(v.shape) == 1:
+        assert angle.shape[0] == 1
+        assert v.shape[0] == 3
+        return np.array(rotfunc(v[0], v[1], v[2], cos, sin), dtype=np.float).T
+    else:
+        assert angle.shape[0] == 1 or angle.shape[0] == v.shape[0]
+        assert v.shape[1] == 3
+        return np.array(rotfunc(v[:, 0], v[:, 1], v[:, 2], cos, sin), dtype=np.float).T

wetb/utils/tests/test_rotation.py

+'''
+Created on 15/01/2014
+
+@author: MMPE
+'''
+import unittest
+
+import numpy as np
+from wetb.utils.geometry import rad
+from wetb.utils.rotation import transformation_matrix, mdot, dots, rotate, rotate_x, \
+    rotmat, rotate_y, rotate_z
+
+x, y, z = 0, 1, 2
+class Test(unittest.TestCase):
+
+
+    def test_transformation_matrix(self):
+        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), x), [[[ 1., 0. , 0.],
+                                                                                  [ 0., 0.8660254, 0.5],
+                                                                                  [ 0., -0.5 , 0.8660254]]])
+        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), y), [[[ 0.8660254, 0. , -0.5],
+                                                                     [ 0., 1., 0 ],
+                                                                     [ 0.5, 0        , 0.8660254]]])
+        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), z), [[[ 0.8660254, 0.5 , 0.],
+                                                                                  [ -0.5, 0.8660254, 0],
+                                                                                  [ 0., 0 , 1]]])
+
+    def test_rotation_matrix(self):
+        np.testing.assert_array_almost_equal(rotmat(rad(30), x), [[[ 1., 0. , 0.],
+                                                                   [ 0., 0.8660254, -0.5],
+                                                                   [ 0., 0.5 , 0.8660254]]])
+        np.testing.assert_array_almost_equal(rotmat(rad(30), y), [[[ 0.8660254, 0. , 0.5],
+                                                                   [ 0., 1., 0 ],
+                                                                   [ -0.5, 0        , 0.8660254]]])
+        np.testing.assert_array_almost_equal(rotmat(rad(30), z), [[[ 0.8660254, -0.5 , 0.],
+                                                                   [ 0.5, 0.8660254, 0],
+                                                                   [ 0., 0 , 1]]])
+    def test_rotation_matrixes(self):
+        np.testing.assert_array_almost_equal(rotmat([rad(30), rad(60)], 0), [[[ 1., 0. , 0.],
+                                                                     [ 0., 0.8660254, -0.5],
+                                                                     [ 0., 0.5 , 0.8660254]],
+                                                                             [[ 1., 0. , 0.],
+                                                                     [ 0., 0.5, -0.8660254, ],
+                                                                     [ 0., +0.8660254, 0.5 ]]])
+
+    def test_mdot(self):
+        x, y, z = 0, 1, 2
+        m1 = transformation_matrix(np.pi, x)
+        m2 = transformation_matrix([np.pi, np.pi / 2], y)
+
+        np.testing.assert_array_almost_equal(m1, [[[1, 0, 0], [0, -1, 0], [0, 0, -1]]])
+        np.testing.assert_array_almost_equal(mdot(m1, m1), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]])
+        np.testing.assert_array_almost_equal(mdot(m1, m2), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, -1], [0, -1, 0], [-1, 0, 0]]])
+        np.testing.assert_array_almost_equal(mdot(m2, m2), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[-1, 0, 0], [0, 1, 0], [0, 0, -1]]])
+        np.testing.assert_array_almost_equal(mdot(m2, m1), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, 1], [0, -1, 0], [1, 0, 0]]])
+
+
+    def test_dots(self):
+        x, y, z = 0, 1, 2
+        v1 = np.array([1, 2, 3])
+        v2 = np.array([1, 2, 4])
+
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, x), v1).T, [[1, -3, 2]])
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , x), v1).T, [[1, -2, -3]])
+
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, y), v1).T, [[3, 2, -1]])
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , y), v1).T, [[-1, 2, -3]])
+
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, z), v1).T, [[-2, 1, 3]])
+        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , z), v1).T, [[-1, -2, 3]])
+
+        v = np.array([v1, v2]).T
+
+        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], x), v).T, [[1, -3, 2], [1, -2, -4]])
+        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], y), v).T, [[3, 2, -1], [-1, 2, -4]])
+        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], z), v).T, [[-2, 1, 3], [-1, -2, 4]])
+        np.testing.assert_array_almost_equal(dots([transformation_matrix(np.pi / 2, z), transformation_matrix(np.pi / 2, y), transformation_matrix(np.pi / 2, x)], v1).T, [[3, -2, 1]])
+
+
+    def test_rotate(self):
+        x, y, z = 0, 1, 2
+        v = np.array([1, 2, 3])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, x), v), [1, -3, 2])
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, x), v), [1, -2, -3])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, y), v), [3, 2, -1])
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, y), v), [-1, 2, -3])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, z), v), [-2, 1, 3])
+        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, z), v), [-1, -2, 3])
+
+        v = np.array([[1, 2, 3], [1, 2, 4]])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[0], [1, -3, 2])
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[1], [1, -2, -4])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[0], [3, 2, -1])
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[1], [-1, 2, -4])
+
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[0], [-2, 1, 3])
+        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[1], [-1, -2, 4])
+
+
+    def test_rotate_xyz(self):
+        x, y, z = 0, 1, 2
+        v = np.array([1, 2, 3])
+
+        np.testing.assert_array_almost_equal(rotate_x(v, np.pi / 2), [1, -3, 2])
+        np.testing.assert_array_almost_equal(rotate_x(v, -np.pi), [1, -2, -3])
+
+        np.testing.assert_array_almost_equal(rotate_y(v, np.pi / 2), [3, 2, -1])
+        np.testing.assert_array_almost_equal(rotate_y(v, np.pi), [-1, 2, -3])
+
+        np.testing.assert_array_almost_equal(rotate_z(v, np.pi / 2), [-2, 1, 3])
+        np.testing.assert_array_almost_equal(rotate_z(v, np.pi), [-1, -2, 3])
+
+        v = np.array([[1, 2, 3], [4, 5, 6]])
+        np.testing.assert_array_almost_equal(rotate_z(v, np.pi / 2), [[-2, 1, 3], [-5, 4, 6]])
+        np.testing.assert_array_almost_equal(rotate_z(v, [np.pi / 2]), [[-2, 1, 3], [-5, 4, 6]])
+        np.testing.assert_array_almost_equal(rotate_z(v, [-np.pi / 2, np.pi]), [[2, -1, 3], [-4, -5, 6]])
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.test_rad']
+    unittest.main()

wetb/wind/air_density.py

+'''
+Created on 08/02/2016
+
+@author: mmpe
+'''
+
+import numpy as np
+
+
+
+
+def saturated_vapor_pressure(T):
+    """Calculate pressure of saturated water vapor at specified temperature as described at
+    http://wahiduddin.net/calc/density_altitude.htm
+    Parameters
+    ---------
+    t : float
+        Temperature [C]
+
+    Returns
+    -------
+    float
+        Pressure [mb]
+    """
+
+    eso = 6.1078
+    c0 = 0.99999683
+    c1 = -0.90826951 * 10 ** -2
+    c2 = 0.78736169 * 10 ** -4
+    c3 = -0.61117958 * 10 ** -6
+    c4 = 0.43884187 * 10 ** -8
+    c5 = -0.29883885 * 10 ** -10
+    c6 = 0.21874425 * 10 ** -12
+    c7 = -0.17892321 * 10 ** -14
+    c8 = 0.11112018 * 10 ** -16
+    c9 = -0.30994571 * 10 ** -19
+    p = (c0 + T * (c1 + T * (c2 + T * (c3 + T * (c4 + T * (c5 + T * (c6 + T * (c7 + T * (c8 + T * (c9))))))))))
+    return eso / p ** 8
+
+def saturated_vapor_pressure2(t):
+    """Calculate pressure of saturated water vapor at specified temperature as described at
+    http://wahiduddin.net/calc/density_altitude.htm
+    Parameters
+    ---------
+    t : float
+        Temperature [C]
+
+    Returns
+    -------
+    float
+        Pressure [mb]
+    """
+    c0 = 6.1078
+    c1 = 7.5
+    c2 = 237.3
+    return c0 * 10 ** ((c1 * t) / (c2 + t))
+
+def saturated_vapor_pressure3(t):
+    """Calculate pressure of saturated water vapor at specified temperature according to
+
+    The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use
+    W. Wagner and A. Pruss
+    J. Phys. Chem. Ref. Data 31, 387 (2002); http://dx.doi.org/10.1063/1.1461829
+
+    Parameters
+    ---------
+    t : float
+        Temperature [C]
+
+    Returns
+    -------
+    float
+        Pressure [mb]
+    """
+
+    T = t + 273.15
+    Tc = 647.096
+    Pc = 220640
+    v = 1 - T / Tc
+    C1 = -7.85951783
+    C2 = 1.84408259
+    C3 = -11.7866497
+    C4 = 22.6807411
+    C5 = -15.9618719
+    C6 = 1.80122502
+    return Pc * np.exp(Tc / T * (C1 * v + C2 * v ** 1.5 + C3 * v ** 3 + C4 * v ** 3.5 + C5 * v ** 4 + C6 * v ** 7.5))
+
+def saturated_vapor_pressure4(t):
+
+    """Calculate the saturated vapor pressure as described in
+    http://www.vaisala.com/Vaisala%20Documents/Application%20notes/Humidity_Conversion_Formulas_B210973EN-F.pdf
+
+    Parameters
+    ---------
+    t : float
+        Temperature [C]
+
+    Returns
+    -------
+    Saturated vapor pressure [mBar]
+    """
+    A = 6.116441
+    m = 7.591386
+    Tn = 240.7263
+    return A * 10 ** ((m * t) / (t + Tn))
+
+def saturated_vapor_pressure_IEC(t):
+    """Calculate the saturated vapor pressure according to IEC 61400-12-1
+
+    Parameters
+    ---------
+    t : float
+        Temperature [C]
+
+    Returns
+    -------
+    Saturated vapor pressure [mBar]
+    """
+
+    T = t + 273.15
+    return 0.000000205 * np.exp(0.0631846 * T)
+
+
+
+def drew_point(t, RH):
+    A = 6.116441
+    m = 7.591386
+    Tn = 240.7263
+    Pw = saturated_vapor_pressure4(t) * RH / 100
+    return Tn / (m / np.log10(Pw / A) - 1)
+
+def air_density(P, t, rh=0, saturated_vapor_pressure_function=saturated_vapor_pressure):
+    """Calculate the density of atmospheric air at specified pressure, temperature and humidity
+    source: http://wahiduddin.net/calc/density_altitude.htm
+    Equivalent to formulation in IEC61400-12-1 if used with the saturated_vapor_pressure_IEC function
+
+
+    Parameters
+    ---------
+    P : float
+        Atmospheric pressure [mb]=[hPa]
+    t : float
+        Temperature [C]
+    rh : float
+        Relative humidity [%]
+    saturated_vapor_pressure_function : function
+        Function, f(t)->P, that takes the temperature in celcius as input and
+        returns the saturated vapor pressure in mbar
+
+    Returns
+    -------
+    float
+        Density of atmospheric air
+    """
+    Pv = saturated_vapor_pressure_function(t) * rh / 100
+    Pd = P - Pv
+    Rv = 461.4964
+    Rd = 287.0531
+    Tk = t + 273.15
+    return (Pd * 100 / (Rd * Tk)) + (Pv * 100 / (Rv * Tk))
+
+def R(rh=0, t=15, P=1014):
+    """Specific gas constant ~287.058 J/(kg K) for dry air
+    
+    Parameters
+    ---------
+    rh : float
+        Relative humidity [%]
+    t : float
+        Temperature [C]
+    P : float
+        pressure [hPa]
+    
+    Returns
+    -------
+    Specific gas constant
+    """
+    assert np.all((900<P)&(P<1100)), "Pressure outside range 900 to 1100"
+    assert np.all((-50<t)&(t<100)), "Temperature outside range -50 to 100"
+    Tk = t + 273.15
+    return P * 100 / (air_density(P, t, rh) * Tk)
+
+if __name__ == "__main__":
+    pass
+    
+#     #print (air_density(1013, 65, 50))
+#     #print (drew_point(40, 50))
+#     import matplotlib.pyplot as plt
+#     if 0:
+#         t = np.arange(-20, 50)
+#         plt.plot(t, saturated_vapor_pressure(t))
+#         plt.plot(t, saturated_vapor_pressure2(t), "--")
+#         plt.plot(t, saturated_vapor_pressure3(t), ":")
+#         plt.show()
+# 
+#     if 0:
+#         t = np.arange(5, 25)
+#         plt.xlabel("Temperature [C]")
+#         plt.ylabel("Density [kg/m3]")
+#         plt.plot(t, air_density(990, t, 50), label='P: 990 hPa')
+#         plt.plot(t, air_density(1018, t, 50), label='P: 1018 hPa')
+#         plt.legend()
+#         plt.show()
+# 
+#     p = 1013
+#     t = 20
+#     rh = 70
+#     print (R(p, t, rh) * (t + 273.15))
+#     print ((p * 100) / air_density(p, t, rh))
+#     
+#     print (R(p,10,80))
+#     print (R(p,100,80))
+# 
+#     if 0:
+#         rh = np.arange(0, 100)
+#         plt.xlabel("Rh [%]")
+#         plt.ylabel("R")
+#         plt.plot(rh, R(990, 20, rh), label='P: 990 hPa, 20 t')
+# 
+#         plt.legend()
+#         plt.show()

wetb/wind/masks.py

+'''
+Created on 06/11/2015
+
+@author: MMPE
+'''
+
+
+def wsp_mask(x, wsp, pm_tolerance):
+    return (wsp - pm_tolerance <= x) & (x < wsp + pm_tolerance)
+
+def wdir_mask(x, wdir, pm_tolerance):
+    return ((x - wdir + pm_tolerance) % 360 >= 0) & ((x - wdir + pm_tolerance) % 360 < 2 * pm_tolerance)
+

wetb/wind/weibull.py

+"""
+Contents
+--------
+- `fit <#wetb.wind.weibull.fit>`_: Fit a weibull distribution, in terms of the parameters k and A, to the provided wind speeds
+- `pdf <#wetb.wind.weibull.pdf>`_: Create Weibull pdf function
+- `random <#wetb.wind.weibull.random>`_: Create a list of n random Weibull distributed values
+"""
+import numpy as np
+import math
+gamma = math.gamma
+def pdf(A, k):
+    """Create Weibull pdf function
+
+    Parameters
+    ----------
+    A : float
+        Scale parameter
+    k : float
+        Shape parameter
+
+    Returns
+    -------
+    pdf-function
+
+    Examples
+    --------
+    >>> from wetb.wind import weibull
+    >>> pdf_func = weibull.pdf(3,4)
+    >>> pdf_func(5)
+    0.131007116969
+    >>> from pylab import arange, plot, show
+    >>> wsp = arange(20)
+    >>> plot(wsp, pdf_func(wsp))
+    >>> show()
+
+    """
+
+    return lambda x: k * x ** (k - 1) / A ** k * np.exp(-(x / A) ** k)
+
+def cdf(A,k):
+    return lambda x: 1-np.exp(-(x/A)**k)
+
+def random(A, k, n):
+    """Create a list of n random Weibull distributed values
+
+    Parameters
+    ----------
+    A : float
+        Scale parameter
+    k : float
+        Shape parameter
+    n : int
+        Number of values
+
+    Returns
+    -------
+    x : array_like, shape (n,)
+        n random Weibull distributed values
+
+    Examples
+    --------
+    >>> from wetb.wind import weibull
+    >>> from pylab import hist, show
+    >>> hist(weibull.random(4,2,1000), 20)
+    >>> show()
+    """
+    return A * np.random.weibull(k, n)
+
+def fit(wsp):
+    """Fit a weibull distribution, in terms of the parameters k and A, to the provided wind speeds
+
+    Parameters
+    ----------
+    wsp : array_like
+        Wind speeds
+
+    Returns
+    -------
+    A : float
+        Scale parameter
+    k : float
+        Shape parameter
+
+    Examples
+    --------
+    >>> from wetb.wind import weibull
+    >>> A,k = weibull.fit(wsp_lst)
+    """
+    res_pr_ms = 2  # number of wind speed bins pr m/s
+
+    pdf, x = np.histogram(wsp, bins=np.arange(0, np.ceil(np.nanmax(wsp)), 1 / res_pr_ms))
+    x = (x[1:] + x[:-1]) / 2
+    N = np.sum(~np.isnan(wsp))
+    pdf = pdf / N * res_pr_ms
+
+    m = lambda n : np.sum(pdf * x ** n / res_pr_ms)
+    from scipy.optimize import newton
+    func = lambda k : gamma(1 / k + 1) ** 2 / gamma(2 / k + 1) - m(1) ** 2 / m(2)
+    k = newton(func, 1)
+    A = m(1) / gamma(1 / k + 1)
+    return A, k
+
+if __name__ == "__main__":
+    from wetb.wind import weibull
+    from pylab import hist, show, plot
+    hist(weibull.random(10, 2, 10000), 20, normed=True)
+    wsp = np.arange(0, 20, .5)
+    plot(wsp, weibull.pdf(10, 2)(wsp))
+    show()