-
Notifications
You must be signed in to change notification settings - Fork 62
/
coldest91.py
88 lines (75 loc) · 2.13 KB
/
coldest91.py
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
import mx.DateTime
import numpy
from scipy import stats
import iemdb, iemplot
COOP = iemdb.connect("coop", bypass=True)
ccursor = COOP.cursor()
def getc(v):
if v > 23.35:
return 'r'
if v > 21.3956:
return 'y'
if v > 18.82:
return 'g'
return 'b'
vals = []
colors = []
mins = []
for yr in range(1900,2010):
data = numpy.zeros( (366,) )
ccursor.execute("""
SELECT day, (high + low)/2.0 from alldata where stationid = 'ia0200' and
day >= '%s-06-01' and day < '%s-06-01' ORDER by day ASC
""" % (yr, yr+1))
i= 0
for row in ccursor:
data[i] = row[1]
i += 1
minv = 100.
pos = 0
for i in range(0,275):
val = numpy.average(data[i:i+91])
if val < minv:
minv = val
pos = i+90
mins.append( minv )
vals.append( pos + 1)
colors.append( getc(minv) )
h_slope, intercept, h_r_value, p_value, std_err = stats.linregress(numpy.arange(0,110), vals)
vals = numpy.array(vals)
labels = []
xticks = []
sts = mx.DateTime.DateTime(2000,6,1)
for i in range(150,330,1):
ts = sts + mx.DateTime.RelativeDateTime(days=(i-1))
if ts.day in [1,]:
labels.append( ts.strftime("%b %d") )
xticks.append( i )
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
v91 = vals - 91
vals[:] = 90
bars = ax.barh(numpy.arange(1900,2010) - 0.3, vals, left=(v91), facecolor='b', edgecolor='b')
i = 0
for bar in bars:
bar.set_facecolor( getc(mins[i]) )
bar.set_edgecolor( getc(mins[i]) )
i += 1
ax.plot([intercept,intercept+(110.0*h_slope)],[1900,2010], c='#000000')
ax.set_xticklabels( labels )
ax.set_xticks( xticks )
ax.set_ylim(1899.5,2025.5)
p1 = plt.Rectangle((0, 0), 1, 1, fc="r")
p2 = plt.Rectangle((0, 0), 1, 1, fc="y")
p3 = plt.Rectangle((0, 0), 1, 1, fc="g")
p4 = plt.Rectangle((0, 0), 1, 1, fc="b")
ax.legend([p1,p2,p3,p4], ["Warmest 25%","","","Coldest 25%"], ncol=4)
leg = plt.gca().get_legend()
ltext = leg.get_texts()
plt.setp(ltext, fontsize='small')
ax.set_title("Ames [1900-2010] Coldest 91 day period")
ax.grid(True)
fig.savefig('test.ps')
import iemplot
iemplot.makefeature('test')