Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
import math
import numpy as np
NaN = float('nan')
is_nan = lambda a: a != a
def dewpoint(tempC, relhum):
"""
Algorithm from Tom Whittaker tempC is the temperature in degrees Celsius,
relhum is the relative humidity as a percentage.
:param tempC: temperature in celsius
:param relhum: relative humidity as a percentage
"""
if tempC is None or relhum is None:
return NaN
gasconst = 461.5
latheat = 2500800.0
dp = 1.0 / (1.0 / (273.15 + tempC) - gasconst * np.log((0.0 + relhum) / 100) /
(latheat - tempC * 2397.5))
return min(dp - 273.15, tempC)
def relhum(airTempK, dewpointTempK):
"""
Algorithm derived by David Hoese from the above
dewpoint(tempC, relhum) function, both parameters are in Kelvin units.
:param airTempK: air temperature in Kelvin
:param dewpointTempK: dewpoint temp in Kelvin
"""
if airTempK == None or dewpointTempK == None:
return NaN
gas_constant = 461.5
latheat = 2500800.0
# Only one section of the equation
latpart = (latheat - (airTempK - 273.15) * 2397.5)
relativehum = 100 * math.e ** ((latpart / airTempK - latpart / dewpointTempK) / gas_constant)
return relativehum
def potentialtemp(airTempK, pressureMB):
"""
Algorithm from David Hoese to calculate potential temperature.
:param airTempK: air temperature in Kelvin
:param pressureMB: air pressure in millibars
"""
if airTempK == None or pressureMB == None:
return NaN
pT = airTempK * (pressureMB.max() / pressureMB) ** .286
return pT
def altimeter(p, alt):
"""Compute altimeter from pressure and altitude.
Converted from code provided by TomW.
:param p: pressure in hPa.
:param alt: altitude of the measurement in meters.
:returns: altimeter in inHg
"""
n = .190284
c1 = .0065 * pow(1013.25, n) / 288.
c2 = alt / pow((p - .3), n)
ff = pow(1. + c1 * c2, 1. / n)
return ((p - .3) * ff * 29.92 / 1013.25)
def dir2txt(val):
"""Convert degrees [0, 360) to a textual representation.
:param val: decimal degrees
>>> dir2txt(0)
'N'
>>> dir2txt(90)
'E'
>>> dir2txt(180)
'S'
>>> dir2txt(270)
'W'
>>> dir2txt(359)
'N'
"""
assert val >= 0 and val < 360, "'%s' out of range" % val
dirs = ("NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW")
if ((val >= 348.75 and val <= 360) or val >= 0 and val < 11.25): return "N"
# 1/2 degree increment between the directions
i = 11.25;
for dir in dirs:
if val >= i and val < (i + 22.5):
return dir
i += 22.5
def wind_vector_components(windspd, winddir):
"""Decompose scalar or list/array polar wind direction and speed data
into the horizontal and vertical vector components and speed vector.
Inputs can be scalar or arrays.
"""
dir_rad = np.deg2rad(winddir)
spd_arr = np.array(windspd)
V_e = spd_arr * np.sin(dir_rad)
V_n = spd_arr * np.cos(dir_rad)
U_spd = np.sqrt(pow(V_e, 2) + pow(V_n, 2))
return V_e, V_n, U_spd
def wind_vector_degrees(vector_east, vector_north):
"""Re-compose horizontal (east/west) and vertical (north/south) vector
components into wind direction in degrees.
Inputs can be scalar or arrays.
"""
rads = np.arctan2(vector_east, vector_north)
winddir = np.rad2deg(rads)
if isinstance(winddir, np.ndarray):
winddir[np.less(winddir, 0)] += 360
elif winddir < 0:
winddir += 360
return winddir % 360
def mean_wind_vector(windspd, winddir):
V_e, V_n, V_spd = wind_vector_components(windspd, winddir)
avg_dir = wind_vector_degrees(np.mean(V_e), np.mean(V_n))
return avg_dir, np.mean(V_spd)