Newer
Older
from metpy.units import units
from metpy.calc import thickness_hydrostatic
class MyGenericException(Exception):
def __init__(self, message):
self.message = message
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
def haversine_np(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
(lon1, lat1) must be broadcastable with (lon2, lat2).
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
def bin_data_by(a, b, bin_ranges):
nbins = len(bin_ranges)
binned_data = []
for i in range(nbins):
rng = bin_ranges[i]
idxs = (b >= rng[0]) & (b < rng[1])
binned_data.append(a[idxs])
return binned_data
def get_bin_ranges(lop, hip, bin_size=100):
bin_ranges = []
delp = hip - lop
nbins = int(delp/bin_size)
for i in range(nbins):
rng = [lop + i*bin_size, lop + i*bin_size + bin_size]
bin_ranges.append(rng)
return bin_ranges
# t must be monotonic increasing
def get_breaks(t, threshold):
t_0 = t[0:t.shape[0]-1]
t_1 = t[1:t.shape[0]]
d = t_1 - t_0
idxs = np.nonzero(d > threshold)
return idxs
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
def pressure_to_altitude(pres, temp, prof_pres, prof_temp, sfc_pres, sfc_temp, sfc_elev):
if not np.all(np.diff(prof_pres) > 0):
raise MyGenericException("target pressure profile must be monotonic increasing")
if pres < prof_pres[0]:
raise MyGenericException("target pressure less than top of pressure profile")
if temp is None:
temp = np.interp(pres, prof_pres, prof_temp)
i_top = np.argmax(np.extract(prof_pres <= pres, prof_pres)) + 1
pres_s = prof_pres.tolist()
temp_s = prof_temp.tolist()
pres_s = [pres] + pres_s[i_top:]
temp_s = [temp] + temp_s[i_top:]
prof_pres = np.array(pres_s)
prof_temp = np.array(temp_s)
i_bot = prof_pres.shape[0] - 1
if sfc_pres > prof_pres[i_bot]:
pres_s = pres_s + [sfc_pres]
temp_s = temp_s + [sfc_temp]
else:
idx = np.argmax(np.extract(prof_pres <= sfc_pres, prof_pres))
if sfc_temp is None:
sfc_temp = np.interp(sfc_pres, prof_pres, prof_temp)
pres_s = prof_pres.tolist()
temp_s = prof_temp.tolist()
pres_s = pres_s[0:idx] + [sfc_pres]
temp_s = temp_s[0:idx] + [sfc_temp]
prof_pres = np.array(pres_s)
prof_temp = np.array(temp_s)
prof_pres = prof_pres * units.hectopascal
prof_temp = prof_temp * units.kelvin
sfc_elev = sfc_elev * units.meter
z = thickness_hydrostatic(prof_pres, prof_temp) + sfc_elev
return z