From e5a653a37bd83d769f859d6ece8bb24fd921dff9 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 21 Dec 2016 11:11:02 +0100
Subject: [PATCH] New stuff
---
wetb/utils/euler.py | 115 ++++++++++++++
wetb/utils/rotation.py | 242 ++++++++++++++++++++++++++++++
wetb/utils/tests/test_rotation.py | 126 ++++++++++++++++
wetb/wind/air_density.py | 221 +++++++++++++++++++++++++++
wetb/wind/masks.py | 13 ++
wetb/wind/weibull.py | 109 ++++++++++++++
6 files changed, 826 insertions(+)
create mode 100644 wetb/utils/euler.py
create mode 100644 wetb/utils/rotation.py
create mode 100644 wetb/utils/tests/test_rotation.py
create mode 100644 wetb/wind/air_density.py
create mode 100644 wetb/wind/masks.py
create mode 100644 wetb/wind/weibull.py
diff --git a/wetb/utils/euler.py b/wetb/utils/euler.py
new file mode 100644
index 0000000..6c2f924
--- /dev/null
+++ b/wetb/utils/euler.py
@@ -0,0 +1,115 @@
+'''
+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:]]
diff --git a/wetb/utils/rotation.py b/wetb/utils/rotation.py
new file mode 100644
index 0000000..2254630
--- /dev/null
+++ b/wetb/utils/rotation.py
@@ -0,0 +1,242 @@
+'''
+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
diff --git a/wetb/utils/tests/test_rotation.py b/wetb/utils/tests/test_rotation.py
new file mode 100644
index 0000000..07f21bc
--- /dev/null
+++ b/wetb/utils/tests/test_rotation.py
@@ -0,0 +1,126 @@
+'''
+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()
diff --git a/wetb/wind/air_density.py b/wetb/wind/air_density.py
new file mode 100644
index 0000000..0af9724
--- /dev/null
+++ b/wetb/wind/air_density.py
@@ -0,0 +1,221 @@
+'''
+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()
diff --git a/wetb/wind/masks.py b/wetb/wind/masks.py
new file mode 100644
index 0000000..a0c8bec
--- /dev/null
+++ b/wetb/wind/masks.py
@@ -0,0 +1,13 @@
+'''
+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)
+
diff --git a/wetb/wind/weibull.py b/wetb/wind/weibull.py
new file mode 100644
index 0000000..43ef018
--- /dev/null
+++ b/wetb/wind/weibull.py
@@ -0,0 +1,109 @@
+"""
+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()
--
GitLab