Permalink
Browse files

sync feature scripts to upstream, hope others like it

  • Loading branch information...
1 parent 3165012 commit fd51a614734153d493b9a8196dee2ceb2127bab1 @akrherz committed Sep 27, 2012
Showing with 994 additions and 540 deletions.
  1. +37 −0 scripts/feature/acc_sigma.py
  2. +14 −11 scripts/feature/asos_1minute_precip.py
  3. +3 −3 scripts/feature/consec90.py
  4. +75 −0 scripts/feature/corn_cond.py
  5. +40 −0 scripts/feature/daily_temp_departure.py
  6. +27 −30 scripts/feature/days_since_T2.py
  7. +31 −0 scripts/feature/days_this_year.py
  8. +0 −65 scripts/feature/days_warmer.py
  9. +49 −0 scripts/feature/doy_accum_warnings.py
  10. +0 −47 scripts/feature/dwpt_bywarn.py
  11. +53 −72 scripts/feature/fair_wx.py
  12. +16 −4 scripts/feature/gdd_departure_since.py
  13. +31 −41 scripts/feature/heat_index_cnt.py
  14. +4 −4 scripts/feature/isuag_soilt.py
  15. +42 −0 scripts/feature/lamkey.py
  16. +39 −0 scripts/feature/last_twoday.py
  17. +38 −0 scripts/feature/longterm.py
  18. +35 −0 scripts/feature/max_daily_precip.py
  19. +51 −0 scripts/feature/max_days_over.py
  20. +43 −0 scripts/feature/max_days_over2.py
  21. +4 −4 scripts/feature/month_depatures.py
  22. +39 −0 scripts/feature/month_total.py
  23. +99 −0 scripts/feature/nhc_forecast_dprog_dt.py
  24. +0 −49 scripts/feature/percentile_monrain.py
  25. +7 −7 scripts/feature/percentile_springhigh.py
  26. +0 −61 scripts/feature/percentile_yrhigh.py
  27. +0 −52 scripts/feature/percentile_yrrain.py
  28. +5 −5 scripts/feature/plot_nam.py
  29. +61 −0 scripts/feature/q-q.py
  30. +29 −0 scripts/feature/weather_feels_like.py
  31. +24 −20 scripts/feature/xy2.py
  32. +21 −18 scripts/feature/yearly_precip_range.py
  33. +38 −24 scripts/feature/yearly_range.py
  34. +39 −23 scripts/feature/ytd_lines.py
View
37 scripts/feature/acc_sigma.py
@@ -0,0 +1,37 @@
+import iemdb
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+import mx.DateTime
+import numpy
+
+sigma = []
+sigma88 = []
+for i in range(1,224):
+ end = mx.DateTime.DateTime(2012,8,11)
+ start = end - mx.DateTime.RelativeDateTime(days=i)
+ print i, start
+ ccursor.execute("""
+ SELECT year, sum(precip) from alldata_ia where station = 'IA0000' and
+ sday BETWEEN '%s' and '%s' GROUP by year
+ """ % (start.strftime("%m%d"), end.strftime("%m%d")))
+ data = numpy.zeros( (ccursor.rowcount,), 'f')
+ for row in ccursor:
+ data[row[0]-1893] = row[1]
+
+ avg = numpy.average(data)
+ stddev = numpy.std(data)
+
+ sigma.insert(0, (data[2012-1893] - avg) / stddev )
+ sigma88.insert(0, (data[1988-1893] - avg) / stddev )
+
+
+import matplotlib.pyplot as plt
+
+(fig, ax) = plt.subplots(1,1)
+
+ax.plot(numpy.arange(224), sigma)
+ax.plot(numpy.arange(224), sigma88)
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
25 scripts/feature/asos_1minute_precip.py
@@ -13,21 +13,24 @@
sprec = numpy.zeros( (3000,), 'f')
acursor.execute("""
- SELECT extract(epoch from (valid at time zone 'EDT')), tmpf, dwpf, drct,
+ SELECT extract(epoch from (valid at time zone 'CDT')), tmpf, dwpf, drct,
sknt, pres1, gust_sknt, precip,
- valid at time zone 'CDT' from t2012_1minute WHERE station = 'DVN'
- and valid BETWEEN '2012-08-04 00:00-05' and '2012-08-04 23:59-05'
+ valid at time zone 'CDT' from t2012_1minute WHERE station = 'LWD'
+ and valid BETWEEN '2012-08-25 00:00-05' and '2012-08-26 23:59-05'
ORDER by valid ASC
""")
+tot = 0
for row in acursor:
#sgust.append( row[6] * 1.15)
offset = row[8].hour * 60 + row[8].minute
+ if row[8].day == 26:
+ offset += 1440
if row[7] > 0:
print offset, row[8], row[7]
- #if row[8].day == 20:
- # offset += 1440
+ tot += row[7]
sprec[offset] = float(row[7] or 0)
+print tot
acc = numpy.zeros( (3000,), 'f')
rate15 = numpy.zeros( (3000,), 'f')
@@ -37,7 +40,7 @@
acc[i] = acc[i-1] + sprec[i]
rate15[i] = numpy.sum(sprec[i-14:i+1]) * 4
rate60[i] = numpy.sum(sprec[i-59:i+1])
- svalid[i] = mx.DateTime.DateTime(2012,8,4, 0) + mx.DateTime.RelativeDateTime(minutes=i)
+ svalid[i] = mx.DateTime.DateTime(2012,8,25, 0) + mx.DateTime.RelativeDateTime(minutes=i)
print acc[818:843]
#print rate60[818:912]
@@ -46,9 +49,9 @@
# print mx.DateTime.DateTime(2012,8,4, 0) + mx.DateTime.RelativeDateTime(minutes=i)
# Figure out ticks
-sts = mx.DateTime.DateTime(2012,8,4, 13, 30)
-ets = mx.DateTime.DateTime(2012,8,4, 14, 30)
-interval = mx.DateTime.RelativeDateTime(minutes=10)
+sts = mx.DateTime.DateTime(2012,8,25, 13, 30)
+ets = mx.DateTime.DateTime(2012,8,26, 14, 30)
+interval = mx.DateTime.RelativeDateTime(minutes=240)
now = sts
xticks = []
xlabels = []
@@ -86,8 +89,8 @@
ax.set_xlim(min(xticks), max(xticks))
ax.legend(loc=2, prop=prop, ncol=2)
ax.set_ylim(0,12)
-ax.set_xlabel("4 August 2012 (CDT)")
-ax.set_title("4 August 2012 Davenport, IA (KDVN) One Minute Rainfall\n2.38 inches between 1:37 and 2:03 PM")
+ax.set_xlabel("25-26 August 2012 (CDT)")
+ax.set_title("25-26 August 2012 Lamoni, IA (KLWD) One Minute Rainfall\n8.07 inches total")
#ax.set_ylim(0,361)
#ax.set_yticks((0,90,180,270,360))
#ax.set_yticklabels(('North','East','South','West','North'))
View
6 scripts/feature/consec90.py
@@ -7,21 +7,21 @@
ccursor.execute("""SELECT year, day, high, low, precip from alldata_ia
- WHERe station = 'IA2203' and day > '1880-01-01'
+ WHERe station = 'IA8706' and day > '1880-01-01'
ORDER by day ASC""")
running = 0
biggest = numpy.zeros( (2013-1880), 'f')
for row in ccursor:
yr = int(row[0])
- if row[2] < 97:
+ if row[2] < 90:
if running > biggest[yr-1880]:
biggest[yr-1880] = running
running = 0
else:
running += 1
#biggest[2012-1880] += 1
-
+print biggest[1921-1880]
import matplotlib.pyplot as plt
View
75 scripts/feature/corn_cond.py
@@ -0,0 +1,75 @@
+import numpy
+import numpy.ma
+import mx.DateTime.ISO
+import Image
+
+def get(state):
+ condition = numpy.ma.zeros((2013-1986,52), 'f')
+ juliandate = numpy.ma.zeros((2013-1986,52), 'f')
+
+ for line in open('corn_condition.csv').readlines()[1:]:
+ tokens = line.replace('"','').split(",")
+ year = int(tokens[1])
+ if tokens[5] != state:
+ continue
+ week = int(tokens[2][6:8])
+ date = mx.DateTime.strptime(tokens[3], '%Y-%m-%d')
+ condcode = tokens[15].strip()
+ condval = tokens[19].replace('"', '').strip()
+ if condval == "":
+ continue
+ #print year, week, date, condcode, condval, float(date.strftime("%j"))
+ if condcode in ['POOR', 'VERY POOR']:
+ condition[year-1986,week-1] += float(condval)
+ juliandate[year-1986,week-1] = float(date.strftime("%j"))
+
+ juliandate.mask = numpy.where(juliandate == 0, True, False)
+ if state == 'ILLINOIS':
+ print numpy.max(condition, 1)
+ return condition, juliandate
+
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.font_manager
+prop = matplotlib.font_manager.FontProperties(size=10)
+
+states = ['IOWA', 'ILLINOIS', 'MINNESOTA', 'WISCONSIN', 'MISSOURI', 'INDIANA']
+fig, ax = plt.subplots(3,2, sharex=True, sharey=True)
+
+i = 0
+for row in range(3):
+ for col in range(2):
+ colors = ['black', 'green', 'blue','red']
+ condition, juliandate = get(states[i])
+
+ for yr in range(1986,2013):
+ if yr in [1993,2012,1988,2005]:
+ ax[row,col].plot( juliandate[yr-1986,:], condition[yr-1986,:], c=colors.pop(),
+ label='%s' % (yr,), lw=3, zorder=2)
+ else:
+ ax[row,col].plot( juliandate[yr-1986,:], condition[yr-1986,:], c='tan',
+ zorder=1)
+ if row == 1 and col == 1:
+ ax[row,col].legend(ncol=4, loc=(-0.75,1.02), prop=prop)
+ ax[row, col].set_xticks( (121,152, 182,213,244,274,305,335,365) )
+ ax[row, col].set_xticklabels( ('May', 'Jun', 'Jul','Aug','Sep','Oct','Nov','Dec') )
+ ax[row, col].set_xlim(120,310)
+ ax[row, col].grid(True)
+ ax[row, col].set_ylim(0,100)
+ if col == 0:
+ ax[row, col].set_ylabel("Coverage [%]")
+ ax[row, col].set_title(" ")
+ ax[row, col].text(126, 85, states[i].capitalize(), size=16)
+ i += 1
+
+#ax[0,0].set_title("USDA Weekly Crop Progress Report (1979-2012)\nIowa Corn Planting Progress (6 years highlighted)")
+
+fig.text(0.5, .91, 'USDA Weekly Corn Crop Condition Report (1986-2012)\nPercentage of State in Poor & Very Poor Condition\n(2012 thru 12 August)', ha='center')
+
+logo = Image.open('../../htdocs/images/logo_small_white.png')
+ax3 = plt.axes([0.10,0.91,0.1,0.1], frameon=False, axisbg=(0.4471,0.6235,0.8117), yticks=[], xticks=[])
+ax3.imshow(logo, origin='lower')
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
40 scripts/feature/daily_temp_departure.py
@@ -0,0 +1,40 @@
+import iemdb
+import numpy
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+
+ccursor.execute("""select day, a.high, c.high, a.low, c.low from alldata_ia a, climate81 c where c.station = a.station and c.station = 'IA2203' and year = 2012 and a.sday = to_char(c.valid, 'mmdd') ORDER by day ASC""")
+
+diff = []
+ldiff = []
+for row in ccursor:
+ diff.append( row[1] - row[2] )
+ ldiff.append( row[3] - row[4] )
+
+import matplotlib.pyplot as plt
+
+(fig, ax) = plt.subplots(2,1, sharex=True)
+
+bars = ax[0].bar(numpy.arange(1,len(diff)+1)-0.5, diff, fc='b', ec='b')
+for bar in bars:
+ if bar.get_y() == 0:
+ bar.set_facecolor('r')
+ bar.set_edgecolor('r')
+ax[0].grid(True)
+ax[0].set_title("1 Jan - 12 Aug 2012 Des Moines Daily Temperature Departure")
+ax[0].set_ylabel("High Temperature $^{\circ}\mathrm{F}$")
+ax[1].set_ylabel("Low Temperature $^{\circ}\mathrm{F}$")
+
+bars = ax[1].bar(numpy.arange(1,len(diff)+1)-0.5, ldiff, fc='b', ec='b')
+for bar in bars:
+ if bar.get_y() == 0:
+ bar.set_facecolor('r')
+ bar.set_edgecolor('r')
+ax[1].set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax[1].set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax[1].set_xlim(0,len(diff)+1)
+ax[1].grid(True)
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
57 scripts/feature/days_since_T2.py
@@ -7,50 +7,47 @@
import mx.DateTime
now = mx.DateTime.now()
-from pyIEM import iemdb
-i = iemdb.iemdb()
-iem = i['iem']
-coop = i['coop']
-
+import iemdb
+IEM = iemdb.connect("iem", bypass=True)
+COOP = iemdb.connect("coop", bypass=True)
+icursor = IEM.cursor()
+ccursor = COOP.cursor()
# Compute normal from the climate database
sql = """
-SELECT c.name, c.network, station, min(max_tmpf), x(s.geom) as lon, y(s.geom) as lat, c.climate_site from summary_2011 s, stations c where s.station = c.id
-and s.day in ('2011-06-06','2011-06-07','2011-06-08') and s.network in ('IA_ASOS', 'AWOS') and max_tmpf > 50 and c.network = s.network
-GROUP by c.name, c.network, station, lat, lon, c.climate_site ORDER by station ASC
+ SELECT c.name, c.network, c.id, max_tmpf - min_tmpf, x(c.geom) as lon,
+ y(c.geom) as lat, c.climate_site from summary_2012 s, stations c where s.iemid = c.iemid
+ and s.day = '2012-09-24' and c.network in ('IA_ASOS', 'AWOS')
+ ORDER by id ASC
"""
lats = []
lons = []
vals = []
-rs = iem.query(sql).dictresult()
-for i in range(len(rs)):
- station = rs[i]['station']
- max_tmpf = rs[i]['min']
- cid = rs[i]['climate_site']
- rs2 = coop.query("""
- SELECT year, high from alldata where stationid = '%s' and day < '2011-06-06' ORDER by day DESC
- """ % (cid.lower(),)).dictresult()
+icursor.execute( sql )
+for row in icursor:
+ station = row[2]
+ ts = mx.DateTime.strptime('20120924', "%Y%m%d")
+ max_tmpf = row[3]
+ cid = row[6]
+ ccursor.execute("""
+ SELECT max(day) from alldata_ia where station = '%s'
+ and high - low >= %.0f
+ """ % (cid, max_tmpf))
running = 0
- for j in range(len(rs2)):
- if rs2[j]['high'] >= max_tmpf:
- running += 1
- else:
- running = 0
- if running > 2:
- vals.append( rs2[j]['year'] )
- break
-
- lats.append( rs[i]['lat'] + (0.0001 * i))
- lons.append( rs[i]['lon'] )
- print '%5s %4s %-20s %3i %s' % (station, rs[i]['network'][-4:], rs[i]['name'], max_tmpf, vals[-1])
+ row2 = ccursor.fetchone()
+ ts2 = mx.DateTime.strptime(row2[0].strftime("%Y%m%d"), "%Y%m%d")
+ vals.append( "%.0f~C~%s" % (max_tmpf, ts2.year,) )
+ lats.append( row[5] )
+ lons.append( row[4] )
+ print '%5s %s' % (station, max_tmpf)
cfg = {
'wkColorMap': 'BlAqGrYeOrRe',
'nglSpreadColorStart': 2,
'nglSpreadColorEnd' : -1,
- '_title' : "Year of Previous Date as Warm as 6 June 2011 (High Temperature)",
- #'_valid' : "%s" % ( now.strftime("%d %b %Y"), ),
+ '_title' : "Amount/Year of Previous Date with as Large Hi-Lo Diff",
+ '_valid' : "since temperature change on 24 September 2012, values in F",
#'lbTitleString' : "[days]",
'_showvalues' : True,
'_format' : '%s',
View
31 scripts/feature/days_this_year.py
@@ -0,0 +1,31 @@
+import iemdb
+import numpy
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+ccursor2 = COOP.cursor()
+
+ccursor.execute("""SELECT year, max(high) from alldata_ia WHERE
+ station = 'IA2203' and year < 2012 GROUP by year ORDER by year ASC""")
+
+years = []
+days = []
+for row in ccursor:
+ years.append( row[0] )
+ ccursor2.execute("""SELECt * from alldata_ia where
+ station = 'IA2203' and year = 2012 and high > %s""" % (row[1],))
+ days.append( ccursor2.rowcount )
+
+import matplotlib.pyplot as plt
+
+(fig, ax) = plt.subplots(1,1)
+
+ax.bar(numpy.array(years)-0.5, days, fc='r', ec='r')
+ax.set_xlim(1886,2013)
+ax.grid(True)
+ax.set_ylabel("2012 Days Warmer than Given Year")
+ax.set_xlabel("Chart average: %.1f days" % (numpy.average(numpy.array(days)),))
+ax.set_title("Des Moines: 1 Jan - 1 Aug 2012 Days with warmer high temperature\nthan warmest daily high temperature for that year (1886-2011)")
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
65 scripts/feature/days_warmer.py
@@ -1,65 +0,0 @@
-# Month percentile
-
-import sys, os, random
-sys.path.append("../lib/")
-import iemplot
-
-import mx.DateTime
-now = mx.DateTime.now()
-
-from pyIEM import iemdb, stationTable
-st = stationTable.stationTable("/mesonet/TABLES/coopClimate.stns")
-i = iemdb.iemdb()
-iem = i['iem']
-mesosite = i['mesosite']
-
-xref = {}
-rs = mesosite.query("SELECT id, climate_site, x(geom) as lon, y(geom) as lat from stations WHERE network in ('IA_ASOS', 'AWOS')").dictresult()
-for i in range(len(rs)):
- xref[ rs[i]['id'] ] = rs[i]
-
-
-# Extract normals
-lats = []
-lons = []
-vals = []
-rs = iem.query("""SELECT station, max_tmpf from summary_2010
- WHERE day = '2010-04-01' and network in ('IA_ASOS', 'AWOS') and
- max_tmpf > 40 and station not in ('AXA','MXO','BNW','CWI','CKP','I75','CNC') """).dictresult()
-for i in range(len(rs)):
- stid = rs[i]['station']
- lats.append( xref[ rs[i]['station'].upper() ]['lat'] )
- lons.append( xref[ rs[i]['station'].upper() ]['lon'] )
- ob = rs[i]['max_tmpf']
- # Find obs
- rs2 = iem.query("""SELECT count(*) as days from summary_2009 where
- station = '%s' and max_tmpf > %s
- """ % (stid, ob)
- ).dictresult()
- print "[%s] 2010 high: %s rank: %s" % (stid, ob, rs2[0]['days'])
- vals.append( float(rs2[0]['days']) )
-
-cfg = {
- 'wkColorMap': 'BlAqGrYeOrRe',
- 'nglSpreadColorStart': 2,
- 'nglSpreadColorEnd' : -1,
- '_title' : "Days in 2009 warmer than 1 Apr 2010",
- '_valid' : "",
- 'lbTitleString' : "[days]",
- '_showvalues' : True,
- '_format' : '%.0f',
- 'pmLabelBarHeightF' : 0.6,
- 'pmLabelBarWidthF' : 0.1,
- 'lbLabelFontHeightF' : 0.025
-}
-# Generates tmp.ps
-fp = iemplot.simple_contour(lons, lats, vals, cfg)
-
-os.system("convert -rotate -90 -trim -border 5 -bordercolor '#fff' -resize 900x700 -density 120 +repage %s.ps %s.png" % (fp,fp))
-if os.environ['USER'] == 'akrherz':
- os.system("xv %s.png" % (fp,))
- sys.exit()
-os.system("convert -rotate -90 -trim -border 5 -bordercolor '#fff' -resize 320x210 -density 120 +repage tmp.ps tmp.png")
-os.system("/home/ldm/bin/pqinsert -p 'plot c 000000000000 lsr_snowfall_thumb.png bogus png' tmp.png")
-os.remove("tmp.png")
-os.remove("tmp.ps")
View
49 scripts/feature/doy_accum_warnings.py
@@ -0,0 +1,49 @@
+import numpy
+import iemdb
+POSTGIS = iemdb.connect('postgis', bypass=True)
+pcursor = POSTGIS.cursor()
+
+counts = numpy.zeros( (5,366), 'f')
+accum = numpy.zeros( (5,366), 'f')
+counts2 = numpy.zeros( (5,366), 'f')
+accum2 = numpy.zeros( (5,366), 'f')
+
+for year in range(2008,2013):
+ pcursor.execute("""SELECT extract(doy from issue) as doy,
+ count(*), sum(case when wfo in ('DMX','DVN','ARX','FSD','OAX') then 1 else 0 end) from
+ sbw_%s WHERE phenomena in ('SV','TO') and status = 'NEW' and significance = 'W' GROUP by doy ORDER by doy ASC""" % (year,))
+ for row in pcursor:
+ counts[year-2008, int(row[0])-1] = row[1]
+ counts2[year-2008, int(row[0])-1] = row[2]
+
+for i in range(366):
+ accum[:,i] = numpy.sum(counts[:,0:i+1], 1)
+ accum2[:,i] = numpy.sum(counts2[:,0:i+1], 1)
+
+import matplotlib.pyplot as plt
+
+(fig, ax) = plt.subplots(2,1, sharex=True)
+
+for i in range(5):
+ e = 366
+ if i == 4:
+ e = 250
+ ax[0].plot(numpy.arange(1,e+1), accum[i,:e], linewidth=2, label='%s' % (2008+i,))
+ ax[1].plot(numpy.arange(1,e+1), accum2[i,:e], linewidth=2, label='%s' % (2008+i,))
+
+ax[1].set_xlim(1,367)
+ax[1].set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax[1].set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax[1].grid(True)
+ax[0].grid(True)
+import matplotlib.font_manager
+prop = matplotlib.font_manager.FontProperties(size=10)
+ax[0].legend(loc=4, ncol=3, prop=prop)
+ax[0].set_title("NWS Severe T'Storm + Tornado Warnings\nUnited States 2008-2012 (thru 6 Sep 2012)")
+ax[1].set_title("NWS Offices covering Iowa [DMX,DVN,ARX,FSD,OAX]")
+ax[0].set_ylabel("Count")
+ax[1].set_ylabel("Count")
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
47 scripts/feature/dwpt_bywarn.py
@@ -1,47 +0,0 @@
-# Extract Dwpt by warning
-
-import mx.DateTime
-from pyIEM import iemdb, stationTable
-i = iemdb.iemdb()
-postgis = i["postgis"]
-asos = i["asos"]
-st = stationTable.stationTable("/mesonet/TABLES/iowa.stns")
-
-output = open('month_dwpf.csv', 'w')
-
-for year in range(1986,2011):
- print year
- rs = postgis.query("""select x(c) as lon, y(c) as lat,
- extract(month from issue) as month, issue, phenomena from
- (select st_centroid(st_union(geom)) as c, wfo, phenomena, eventid,
- issue from warnings_%s WHERE ugc ~* '^IAC' and
- phenomena in ('TO','SV') and significance ='W' and gtype = 'C'
- GROUP by issue, phenomena, eventid, wfo) as foo""" % (year,)
- ).dictresult()
- for i in range(len(rs)):
- lat = rs[i]['lat']
- lon = rs[i]['lon']
- month = rs[i]['month']
- phenomena = rs[i]['phenomena']
- ts = mx.DateTime.strptime(rs[i]['issue'][:16], '%Y-%m-%d %H:%M')
-
- obs = asos.query("""SELECT * from t%s where valid BETWEEN '%s' and '%s'
- and station in %s and dwpf > 0 and tmpf > 0 ORDER by valid ASC""" % (year,
- (ts - mx.DateTime.RelativeDateTime(minutes=60)).strftime("%Y-%m-%d %H:%M"),
- (ts + mx.DateTime.RelativeDateTime(minutes=10)).strftime("%Y-%m-%d %H:%M"),
- `tuple(st.ids)`)).dictresult()
- dist = 1000
- val = None
- for j in range(len(obs)):
- slat = st.sts[obs[j]['station']]['lat']
- slon = st.sts[obs[j]['station']]['lon']
- sdist = ((slat - lat)**2 + (slon - lon)**2)**.5
- if slon < lon and slat > lat and sdist > 0.2: # Nothing NW please
- continue
- if sdist < dist:
- dist = sdist
- val = obs[j]
- if val:
- output.write("%s,%s,%s,%.3f,%s,%s\n" % (month, val['dwpf'], phenomena, dist, val['tmpf'] - val['dwpf'], ts.strftime("%Y-%m-%d %H:%M")))
-
-output.close()
View
125 scripts/feature/fair_wx.py
@@ -140,8 +140,9 @@
import numpy
import numpy.ma
-#avgT = numpy.ma.zeros( (2012-1880) )
-precip = numpy.ma.zeros( (2012-1880) )
+avgH = numpy.ma.zeros( (2013-1880) )
+avgL = numpy.ma.zeros( (2013-1880) )
+precip = numpy.ma.zeros( (2013-1880) )
precip[:] = -0.000000001
#heatcnt = numpy.ma.zeros( (2012-1880) )
#heatcnt[:] = -1
@@ -157,97 +158,77 @@
#if sts.year < 1933:
# continue
ccursor.execute("""
- SELECT day, precip from alldata where stationid = 'ia2203' and
- day >= '%s' and day <= '%s' ORDER by day ASC
+ SELECT avg(high), avg(low) from alldata_ia where station = 'IA2203' and
+ day >= '%s' and day <= '%s'
""" % (sts.strftime("%Y-%m-%d"), ets.strftime("%Y-%m-%d")))
i = 0
- for row in ccursor:
- if row[1] > 0.04:
- hits[i] += 1
- if row[1] > 0.00:
- total[i] += row[1]
- rains[i] += 1
- precip[sts.year-1880] += row[1]
- else:
- precip[sts.year-1880] += 0.001
- cnts[i] += 1
- i += 1
- #if row[0] > 50:
- # avgT[sts.year-1880] = row[0]
+ row = ccursor.fetchone()
+ avgH[sts.year-1880] = row[0]
+ avgL[sts.year-1880] = row[1]
# Heat index
- if sts.year < 1933:
- continue
- acursor.execute("""
- SELECT date(valid) as d, max(dwpf), min(dwpf) from t%s
- WHERE station = 'DSM' and tmpf > 0 and dwpf > 0 and
- valid BETWEEN '%s 00:00' and '%s 23:59' GROUP by d ORDER by d ASC
- """ % (sts.year, sts.strftime("%Y-%m-%d"), ets.strftime("%Y-%m-%d")))
- keys = {}
- cnt = 0
- data = ['M']*12
- data2 = ['M']*12
- cnt = 0
- for row in acursor:
- data[cnt] = str(row[1])
- data2[cnt] = str(row[2])
- cnt += 1
+# if sts.year < 1933:
+# continue
+# acursor.execute("""
+# SELECT date(valid) as d, max(dwpf), min(dwpf) from t%s
+# WHERE station = 'DSM' and tmpf > 0 and dwpf > 0 and
+# valid BETWEEN '%s 00:00' and '%s 23:59' GROUP by d ORDER by d ASC
+# """ % (sts.year, sts.strftime("%Y-%m-%d"), ets.strftime("%Y-%m-%d")))
+# keys = {}
+# cnt = 0
+# data = ['M']*12
+# data2 = ['M']*12
+# cnt = 0
+# for row in acursor:
+# data[cnt] = str(row[1])
+# data2[cnt] = str(row[2])
+# cnt += 1
#h = mesonet.heatidx(row[1], mesonet.relh(row[1], row[2]))
#if h >= 90:
# keys[ row[0].strftime("%Y%m%d%H") ] = 1
- dwpfs.write("%s,%s,%s\n" % (`sts.year`, ",".join(data), ",".join(data2)))
+# dwpfs.write("%s,%s,%s\n" % (`sts.year`, ",".join(data), ",".join(data2)))
#if cnt > 50:
# heatcnt[sts.year-1880] = len(keys.keys())
-dwpfs.close()
+#dwpfs.close()
-decades = numpy.ma.zeros( (2012-1880) )
+#decades = numpy.ma.zeros( (2012-1880) )
-for decade in range(1880,2020,10):
- avg = numpy.ma.average( precip[decade-1880:decade-1880+10] )
- decades[decade-1880:decade-1880+10] = avg
+#for decade in range(1880,2020,10):
+# avg = numpy.ma.average( precip[decade-1880:decade-1880+10] )
+# decades[decade-1880:decade-1880+10] = avg
#avgT[2011-1880] = 88
-#avgT.mask = numpy.where(avgT == 0, True, False)
+avgH.mask = numpy.where(avgH == 0, True, False)
+avgL.mask = numpy.where(avgL == 0, True, False)
#heatcnt.mask = numpy.where(heatcnt == -1, True, False)
-avgP = numpy.ma.average(precip)
-precip = numpy.where(precip < 0, 9, precip)
+avgHH = numpy.ma.average(avgH)
+avgLL = numpy.ma.average(avgL)
+#precip = numpy.where(precip < 0, 9, precip)
import matplotlib.pyplot as plt
-fig = plt.figure()
-ax = fig.add_subplot(211)
-bars = ax.bar(numpy.arange(1880,2012)-0.4, precip, edgecolor='r', facecolor='r')
+(fig, ax) = plt.subplots(2,1, sharex=True)
+bars = ax[0].bar(numpy.arange(1880,2013)-0.4, avgH, edgecolor='r', facecolor='r')
for bar in bars:
- if bar.get_height() == 9:
- bar.set_facecolor("#EEEEEE")
- bar.set_edgecolor("#EEEEEE")
- elif bar.get_height() > avgP:
- bar.set_facecolor("b")
- bar.set_edgecolor("b")
-ax.plot(numpy.arange(1880,2012),decades, color='black', label='Decade Average')
-bars[-1].set_facecolor('blue')
-bars[-1].set_edgecolor('blue')
-ax.set_xlim(1892.5,2011.5)
-#ax.set_ylim(numpy.ma.min(avgT) - 3, numpy.ma.max(avgT) + 3)
-ax.grid(True)
-ax.legend(loc=2)
-ax.set_title("Iowa State Fair Precipitation")
-ax.set_ylabel("Total Precipitation [inch]")
-ax.set_xlabel("*2011 data thru 17 August")
+ if bar.get_height() < avgHH:
+ bar.set_facecolor('b')
+ bar.set_edgecolor('b')
+ax[0].set_xlim(1879.5,2012.5)
+ax[0].set_ylim(70,100)
+ax[0].grid(True)
+ax[0].set_title("Iowa State Fair Average Daily Temps (Des Moines wx site)")
+ax[0].set_ylabel("High Temperature $^{\circ}\mathrm{F}$")
-ax2 = fig.add_subplot(212)
-ax2.bar(numpy.arange(1,13)-0.4, hits / cnts * 100. )
-ax2.set_ylabel("Frequency of 0.05+ inches [%]")
-#ax2.bar(numpy.arange(1880,2012)-0.4, heatcnt, edgecolor='r', facecolor='r')
-ax2.set_xlim(0.5,11.5)
-ax2.set_xticks( numpy.arange(1,12) )
-ax2.set_xticklabels( numpy.arange(1,12) )
-#ax2.set_yticks( numpy.arange(0,5*24,24))
-ax2.grid(True)
-#ax2.set_ylabel("Hours with\n Heat Index over 90")
-ax2.set_xlabel("Day of Fair")
+bars = ax[1].bar(numpy.arange(1880,2013)-0.4, avgL, edgecolor='r', facecolor='r')
+for bar in bars:
+ if bar.get_height() < avgLL:
+ bar.set_facecolor('b')
+ bar.set_edgecolor('b')
+ax[1].set_ylim(50,75)
+ax[1].grid(True)
+ax[1].set_ylabel("Low Temperature $^{\circ}\mathrm{F}$")
fig.savefig('test.ps')
import iemplot
View
20 scripts/feature/gdd_departure_since.py
@@ -21,6 +21,7 @@
jday = int(mx.DateTime.now().strftime("%j")) - 1
x = []
y = []
+maxy = []
xticks = []
xticklabels = []
for i in range(0,jday+1):
@@ -34,13 +35,21 @@
xticks.append( i )
xticklabels.append( ts.strftime("%-d %b"))
+ # Compute max accumulation
+ ccursor.execute("""SELECT year, sum(gdd50(high,low)) from alldata_ia
+ where station = 'IA0200' and sday BETWEEN '%s' and '0905'
+ GROUP by year ORDER by sum DESC LIMIT 1
+ """ % (ts.strftime("%m%d"),))
+ row = ccursor.fetchone()
+ maxy.append( float(row[1]) - climate_total)
+
import numpy
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1)
ax2 = ax.twinx()
-ax2.set_ylabel("GDD Departure [Climatology of 15 May days]")
+ax2.set_ylabel("GDD Departure [Climatology of 15 Jul days]")
def Tc(Tf):
- return Tf / 11.6
+ return Tf / 23.1
def update_ax2(ax1):
y1, y2 = ax1.get_ylim()
@@ -51,13 +60,16 @@ def update_ax2(ax1):
ax.callbacks.connect("ylim_changed", update_ax2)
ax.bar(numpy.array(x)-0.5,y, fc='r', ec='r')
+ax.plot(numpy.array(x),maxy, color='k', label='Maximum Departure')
+
ax.set_xticks(xticks)
ax.set_xlim(-0.5,max(x)+0.5)
ax.set_xticklabels( xticklabels )
-ax.set_xlabel("Period from 14 May back to date")
+ax.set_xlabel("Period from 5 Sep back to date")
ax.set_ylabel("Growing Degree Days Departure (base=50)")
-ax.set_title("2012 Ames Growing Degree Day Departure (base=50)\nfrom climatology over periods prior to today")
+ax.set_title("2012 Ames Growing Degree Day Departure (base=50)\nfrom climatology over periods prior to today, max 1893-2012")
ax.grid(True)
+ax.legend(loc='best')
fig.savefig('test.ps')
import iemplot
View
72 scripts/feature/heat_index_cnt.py
@@ -1,46 +1,36 @@
-import matplotlib.pyplot as plt
-fig = plt.figure()
-ax = fig.add_subplot(111)
-
-
-from pyIEM import iemdb, mesonet
-i = iemdb.iemdb()
-asos = i['asos']
+import mesonet
+import iemdb
+ASOS = iemdb.connect('iem', bypass=True)
+acursor = ASOS.cursor()
import mx.DateTime
import numpy
+import network
+nt = network.Table(("IA_ASOS", 'AWOS'))
+
+acursor.execute("""SELECT t.id, valid, tmpf, dwpf from current_log c JOIN stations t on
+ (t.iemid = c.iemid) WHERE t.network in ('IA_ASOS','AWOS') and valid
+ BETWEEN '2012-09-04' and '2012-09-05'
+ and tmpf > 0 and dwpf > 0 """)
-data = {}
-cnts = numpy.zeros( (2012-1971), 'f')
-for yr in range(1971,2012):
- rs = asos.query("""SELECT valid, tmpf, dwpf from t%s
- WHERE station = 'DSM' and tmpf > 0 and dwpf > 0 and extract(month from valid) = 7
- ORDER by valid ASC""" % (yr,)).dictresult()
- for i in range(len(rs)):
- ts = mx.DateTime.strptime(rs[i]['valid'][:16], '%Y-%m-%d %H:%M')
- key = ts.strftime("%Y%m%d%H")
- if data.has_key(key):
- continue
- data[key] = 0
- h = mesonet.heatidx(rs[i]['tmpf'], mesonet.relh(rs[i]['tmpf'], rs[i]['dwpf']))
- if h >= 100:
- cnts[yr-1971] += 1
+maxes = {}
+for row in acursor:
+ h = mesonet.heatidx(row[2], mesonet.relh(row[2], row[3]))
+ maxes[row[0]] = max(h, maxes.get(row[0], 0))
-avgV = numpy.average(cnts)
-bars = ax.bar( numpy.arange(1971,2012) -0.4, cnts , fc='r', ec='r')
-ax.plot( [1971,2012], [avgV,avgV] )
-for bar in bars:
- if bar.get_height() < avgV:
- bar.set_edgecolor('b')
- bar.set_facecolor('b')
-#bars[-1].set_edgecolor("b")
-#bars[-1].set_facecolor("b")
-ax.set_ylabel("Hours over 100 F")
-ax.set_xlim(1970.5,2011.5)
-ax.set_ylim(0,5*24)
-ax.set_yticks( numpy.arange(0,5*24,24) )
-#ax.set_xlabel("*2011 Total thru 15 July")
-ax.set_title("Des Moines *July* Hours over 100 F Heat Index [1971-2011]")
-ax.grid(True)
+vals = []
+lons = []
+lats = []
+for sid in maxes.keys():
+ vals.append( maxes[sid])
+ lons.append( nt.sts[sid]['lon'] )
+ lats.append( nt.sts[sid]['lat'] )
+
+cfg = {'wkColorMap': 'WhiteYellowOrangeRed',
+ '_title': '4 Sep 2012 ASOS/AWOS Maximum Heat Index',
+ 'lbTitleString': 'F',
+ '_showvalues': True,
+ '_format': '%.0f'
+ }
import iemplot
-fig.savefig('test.ps')
-iemplot.makefeature('test')
+tmpfp = iemplot.simple_contour(lons, lats, vals, cfg)
+iemplot.makefeature(tmpfp)
View
8 scripts/feature/isuag_soilt.py
@@ -17,7 +17,7 @@
icursor.execute("""
select valid, c30
from daily WHERE station = 'A130209'
- and valid >= '%s-01-01' and valid < '%s-05-01' ORDER by valid ASC
+ and valid >= '%s-01-01' and valid < '%s-09-23' ORDER by valid ASC
""" % (yr,yr) )
for row in icursor:
x.append( int(row[0].strftime("%j")) )
@@ -31,13 +31,13 @@
color = 'r'
ax.plot(x, y, color=color)
-ax.set_title("ISU AgClimate Ames Site 4 inch Soil Temperature\n1 Jan - 30 Apr Yearly Timeseries [1988-2012]")
-ax.set_xlabel("Red Line is 1 Jan - 23 Apr 2012")
+ax.set_title("ISU AgClimate Ames Site 4 inch Soil Temperature\n1 Jan - 22 Sep Yearly Timeseries [1988-2012]")
+ax.set_xlabel("Red Line is 1 Jan - 22 Sep 2012")
ax.grid(True)
ax.set_ylabel('Daily Avg Temp $^{\circ}\mathrm{F}$')
ax.set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
ax.set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
-ax.set_xlim(0,121)
+ax.set_xlim(0,267)
fig.savefig('test.ps')
import iemplot
iemplot.makefeature('test')
View
42 scripts/feature/lamkey.py
@@ -0,0 +1,42 @@
+import iemdb
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+
+import network
+nt = network.Table("IACLIMATE")
+
+ccursor.execute("""
+--- select station, avg(sum) from (select station, year, sum(precip) from alldata_ia where station in (select distinct station from alldata_ia where year = 1931 and precip > 0) and year between 1981 and 2010 GROUP by station, year) as foo GROUP by station
+
+select station, sum(precip) from climate71 GROUP by station
+""")
+previous = {}
+for row in ccursor:
+ previous[row[0]] = row[1]
+
+ccursor.execute("""
+select station, avg(sum) from (select station, year, sum(precip) from alldata_ia where station in (select distinct station from alldata_ia where year = 1931 and precip > 0) and year between 1931 and 1960 GROUP by station, year) as foo GROUP by station
+""")
+
+lats = []
+lons = []
+vals = []
+
+for row in ccursor:
+ if not nt.sts.has_key(row[0]):
+ continue
+ vals.append( previous[row[0]] )
+ lats.append( nt.sts[row[0]]['lat'] )
+ lons.append( nt.sts[row[0]]['lon'] )
+
+import iemplot
+
+cfg = {
+ '_title': '1981-2010 Yearly Average Precipitation',
+ 'lbTitleString': 'in',
+ '_showvalues': True,
+ '_format': '%.2f',
+ }
+
+tmpfp = iemplot.simple_contour(lons, lats, vals, cfg)
+iemplot.makefeature(tmpfp)
View
39 scripts/feature/last_twoday.py
@@ -0,0 +1,39 @@
+import iemdb
+import mx.DateTime
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+
+ccursor.execute("select high, max(sday) from (select day, sday, case when foo.high > foo2.high then foo2.high else foo.high end from (SELECT day, high from alldata_ia where station = 'IA2203') as foo, (SELECT day + '1 day'::interval as day2, sday, high from alldata_ia where station = 'IA2203') as foo2 WHERE foo.day = foo2.day2 ) as f GROUP by high ORDER by high DESC")
+
+threshold = []
+doy = []
+for row in ccursor:
+ if len(threshold) > 0 and row[0] + 1 != threshold[-1]:
+ threshold.append( row[0] + 1)
+ doy.append( doy[-1])
+ threshold.append( row[0] )
+ ts = mx.DateTime.strptime("2000"+row[1], '%Y%m%d')
+ jday = int(ts.strftime("%j"))
+ if len(doy) > 0 and jday < doy[-1]:
+ jday = doy[-1]
+ doy.append( jday )
+
+import matplotlib.pyplot as plt
+import numpy
+
+(fig, ax) = plt.subplots(1,1)
+
+ax.barh( numpy.array(threshold)-0.4, doy, fc='pink', ec='pink')
+ax.set_ylim(60,120)
+ax.grid(True)
+ax.set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax.set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax.set_xlim(min(doy),366)
+ax.set_title("Des Moines 1880-2012 Latest Date for\nTwo Days Above High Temperature")
+ax.text(242, 96, 'x', color='r', fontsize=15, va='center', ha='center')
+ax.set_ylabel("High Temperature $^{\circ}\mathrm{F}$")
+ax.set_xlabel("Red X is 29-30 August 2012 at 96$^{\circ}\mathrm{F}$")
+
+fig.savefig('test.ps')
+import iemplot
+iemplot.makefeature('test')
View
38 scripts/feature/longterm.py
@@ -0,0 +1,38 @@
+import iemdb, mx.DateTime
+ASOS = iemdb.connect('asos', bypass=True)
+acursor = ASOS.cursor()
+
+acursor.execute("""
+ SELECT extract(doy from one.d), two.avg - one.avg from
+ (SELECT date(valid) as d, avg(tmpf) from t2010 WHERE station = 'LWD'
+ GROUP by d) as one,
+ (SELECT date(valid) as d, avg(tmpf) from t2010 WHERE station = 'CSQ'
+ GROUP by d) as two WHERE one.d = two.d ORDER by one.d ASC
+ """)
+times = []
+vals = []
+for row in acursor:
+ times.append( row[0] )
+ vals.append( row[1] )
+
+import matplotlib.pyplot as plt
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(times, vals)
+#ax.legend()
+xticks = []
+xticklabels = []
+for i in range(0,366):
+ ts = mx.DateTime.DateTime(2000,1,1) + mx.DateTime.RelativeDateTime(days=i)
+ if ts.day == 1:
+ xticks.append(i)
+ xticklabels.append( ts.strftime("%b") )
+ax.set_xticks(xticks)
+ax.set_xticklabels(xticklabels)
+ax.set_xlim(0,366)
+ax.set_title("2011 Average Daily Temperature Difference\nAnkeny KIKV (AWOS) minus Des Moines KDSM (ASOS)")
+ax.set_ylabel("Temperature Difference [F]")
+
+ax.grid(True)
+fig.savefig('test.png')
View
35 scripts/feature/max_daily_precip.py
@@ -0,0 +1,35 @@
+import iemdb
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+import network
+nt = network.Table("IACLIMATE")
+
+ccursor.execute("""
+ SELECT one.station, max(one.precip + two.precip) from
+ (SELECT station, day as d1, precip from alldata_ia where year = 2012
+ and day >= '2012-05-01') as one,
+ (SELECT station, day + '1 day'::interval as d2, precip from alldata_ia where year = 2012
+ and day >= '2012-05-01') as two WHERE one.station = two.station
+ and one.d1 = two.d2 GROUP by one.station
+
+""")
+
+lats = []
+lons = []
+vals = []
+for row in ccursor:
+ if row[0][2] == 'C' or row[0][2:] == '0000':
+ continue
+ lats.append( nt.sts[row[0]]['lat'] )
+ lons.append( nt.sts[row[0]]['lon'] )
+ vals.append( row[1] )
+
+cfg = {'lbTitleString': 'inch',
+ '_title': "Max"
+
+ }
+
+import iemplot
+tmpfp = iemplot.simple_contour(lons, lats, vals, cfg)
+
+iemplot.makefeature(tmpfp)
View
51 scripts/feature/max_days_over.py
@@ -0,0 +1,51 @@
+import iemdb, iemplot
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+
+ts = []
+d2012s = []
+ms = []
+ys = []
+for thres in range(70,106):
+ ccursor.execute("""
+ select year, count(*) from alldata_ia where station = 'IA2203'
+ and high >= %s GROUP by year
+ ORDER by count DESC
+ """ % (thres,))
+ m = None
+ years = []
+ d2012 = None
+ for row in ccursor:
+ if m is None:
+ m = row[1]
+ if m == row[1]:
+ years.append( `row[0]` )
+ if row[0] == 2012:
+ d2012 = row[1]
+ print '%s %s %s %s' % (thres, d2012, m, ",".join(years))
+ ts.append( thres )
+ d2012s.append( d2012 )
+ ms.append( m )
+ ys.append( ",".join(years) )
+
+import numpy
+import matplotlib.pyplot as plt
+(fig, ax) = plt.subplots(1,1, figsize=(9,9))
+
+ax.barh(numpy.arange(70,106)-0.4, ms, height=0.8, fc='#efbd47', ec='#efbd47', label='Max')
+ax.barh(numpy.arange(70,106)-0.3, d2012s, height=0.6, fc='#c60c30', ec='#c60c30' , label='2012')
+for (yr, val, m) in zip(ys, numpy.arange(70,106), ms):
+ c = '#000000'
+ if yr.find("2012") > -1:
+ c = "#0000ff"
+ ax.text(m+1, val, "%s - %s" % (m, yr), va='center', size=10, color=c)
+ax.set_ylim(69,106)
+ax.set_xlim(0,225)
+ax.set_xlabel("Days per Year, thru 25 Sep 2012")
+ax.set_ylabel("Daily High Temperature $^{\circ}\mathrm{F}$")
+ax.set_title("Des Moines Maximum Number of Days\nAt or Above given High Temperature (1880-2012)")
+ax.legend()
+ax.grid()
+
+fig.savefig('test.ps')
+iemplot.makefeature('test')
View
43 scripts/feature/max_days_over2.py
@@ -0,0 +1,43 @@
+import iemdb, iemplot
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+import numpy
+
+data = numpy.zeros( (2013-1880,56), 'f')
+
+for thres in range(50,106):
+ ccursor.execute("""
+ select year, count(*) from alldata_ia where station = 'IA2203'
+ and high >= %s and year > 1879 and sday < '0912' GROUP by year
+ """ % (thres,))
+ for row in ccursor:
+ data[row[0]-1893,thres-50] = row[1]
+
+
+import matplotlib.pyplot as plt
+(fig, ax) = plt.subplots(1,1)
+
+cfg = {
+ 2012: {'c': '#FF0000', 'z': 3},
+ 1936: {'c': '#0000FF', 'z': 3},
+ 1934: {'c': '#FF00FF', 'z': 3},
+ 1931: {'c': '#000000', 'z': 3},
+ }
+
+for yr in range(1880,2013):
+ ax.plot(numpy.arange(50,106), data[yr-1893,:],
+ color=cfg.get(yr, {'c':'#bceebc'}).get('c'),
+ zorder=cfg.get(yr, {'z': 1}).get('z'),
+ label=cfg.get(yr, {'l': None}).get('l', `yr`),
+ linewidth=cfg.get(yr, {'z': 1}).get('z'))
+
+#ax.set_ylim(49,106)
+#ax.set_xlim(0,360)
+ax.set_ylabel("Days per Year")
+ax.set_xlabel("Daily High Temperature $^{\circ}\mathrm{F}$")
+ax.set_title("1 Jan - 11 Sep: Des Moines Maximum Number of Days\nAt or Above given High Temperature (1880-2012)")
+ax.legend()
+ax.grid()
+
+fig.savefig('test.ps')
+iemplot.makefeature('test')
View
8 scripts/feature/month_depatures.py
@@ -24,7 +24,7 @@
diff = []
ccursor.execute("""
SELECT year, month, avg((high+low)/2.0) from alldata_ia where
- station = 'IA0000' and year > 2006
+ station = 'IA0000' and year > 2006 and day < '2012-08-16'
GROUP by year, month ORDER by year, month ASC
""")
for row in ccursor:
@@ -50,11 +50,11 @@
xticklabels.append( ts.strftime(fmt) )
xticks.append( i )
-bars = ax.bar(numpy.arange(0, len(diff))-0.4, diff, fc='b', ec='b')
+bars = ax.bar(numpy.arange(0, len(diff))-0.4, diff, fc='r', ec='r')
for bar in bars:
if bar.get_xy()[1] < 0:
- bar.set_facecolor('r')
- bar.set_edgecolor('r')
+ bar.set_facecolor('b')
+ bar.set_edgecolor('b')
ax.set_ylabel("Departure $^{\circ}\mathrm{F}$")
ax.set_xlabel("* Aug 2012 total thru 15 Aug")
ax.grid(True)
View
39 scripts/feature/month_total.py
@@ -0,0 +1,39 @@
+import iemdb
+import network
+nt = network.Table( ("IA_COOP", "MN_COOP", 'NE_COOP', 'MO_COOP') )
+IEM = iemdb.connect('iem', bypass=True)
+icursor = IEM.cursor()
+
+icursor.execute("""
+ SELECT id, count(id) as cnt,
+ sum(CASE WHEN pday >= 0 THEN pday ELSE 0 END) as prectot
+ from summary_2012 s JOIN stations t on (t.iemid = s.iemid)
+ WHERE day >= '2012-05-01' and day < '2012-08-01'
+ and pday >= 0 and t.network in ('IA_COOP', 'MN_COOP', 'NE_COOP', 'MO_COOP')
+ GROUP by id ORDER by prectot ASC
+""")
+
+lats = []
+lons = []
+vals = []
+mask = []
+for row in icursor:
+ if row[1] < 80:
+ continue
+ mask.append( row[0][-2:] == 'I4' )
+ lats.append( nt.sts[row[0]]['lat'] )
+ lons.append( nt.sts[row[0]]['lon'] )
+ vals.append( row[2] )
+
+cfg = {'wkColorMap': 'WhiteBlueGreenYellowRed',
+ '_showvalues': True,
+ '_format': '%.2f',
+ 'lbTitleString': '[inch]',
+ '_valuemask': mask,
+ '_title' : 'May-July 2012 NWS COOP Rainfall Totals',
+ }
+
+import iemplot
+tmpfp = iemplot.simple_contour(lons, lats, vals, cfg)
+
+iemplot.makefeature(tmpfp)
View
99 scripts/feature/nhc_forecast_dprog_dt.py
@@ -0,0 +1,99 @@
+from mpl_toolkits.basemap import Basemap, cm
+import matplotlib
+import matplotlib.pyplot as plt
+from matplotlib.colors import rgb2hex
+from matplotlib.patches import Polygon
+from matplotlib.collections import PatchCollection
+import iemplot
+
+#fig = plt.figure(num=None, figsize=(11.3,7.00))
+fig = plt.figure(num=None, figsize=(3.85,3.1))
+# Set the axes instance
+ax = plt.axes([0.05,0.05,0.9,0.9], axisbg=(0.4471,0.6235,0.8117))
+
+# Set up the map
+#m = Basemap(projection='lcc', urcrnrlat=47.7, llcrnrlat=23.08, urcrnrlon=-62.5,
+# llcrnrlon=-120, lon_0=-98.7, lat_0=39, lat_1=33, lat_2=45,
+# resolution='l')
+m = Basemap(projection='lcc',
+ urcrnrlat=34.7, llcrnrlat=16.08, urcrnrlon=-37.5, llcrnrlon=-95,
+ lon_0=-88.7, lat_0=29, lat_1=23, lat_2=35,
+ resolution='l', ax=ax)
+m.fillcontinents(color='0.7',zorder=0)
+m.drawcountries(zorder=1)
+m.drawstates(zorder=1)
+m.drawcoastlines(zorder=1)
+
+import iemdb
+AFOS = iemdb.connect('afos', bypass=True)
+acursor = AFOS.cursor()
+
+acursor.execute("""SELECT data, entered from products_2012_0712 WHERE pil = 'TCDAT4'
+ ORDER by entered ASC""")
+
+verify_x = []
+verify_y = []
+verify_speed = []
+
+def get_color(speed):
+ if speed < 40:
+ return '#FFFFFF',1
+ if speed < 75:
+ return '#0000FF',2
+ if speed < 96:
+ return '#00FF00',3
+ if speed < 111:
+ return '#FF0000',4
+ if speed < 130:
+ return '#0000FF',5
+ if speed < 157:
+ return '#00FFFF',6
+
+for row in acursor:
+ print row[1]
+ track_x = []
+ track_y = []
+ colors = []
+ tokens = row[0].split("FORECAST POSITIONS AND MAX WINDS")
+ for line in tokens[1].split("\n"):
+ parts = line.strip().split()
+ if len(parts) != 8:
+ continue
+ lat = float(parts[2].replace("N",""))
+ lon = 0 - float(parts[3].replace("W",""))
+ x,y = m(lon,lat)
+ speed = float(parts[6])
+ #ax.text(x,y, parts[6])
+ if parts[0] == "INIT":
+ verify_x.append( lon )
+ verify_y.append( lat )
+ verify_speed.append( speed )
+ track_x.append( lon )
+ track_y.append( lat )
+ c, z = get_color(speed)
+ if parts[0] != "INIT":
+ if track_x[-2] - track_x[-1] > 5:
+ continue
+ x,y = m(track_x[-2:], track_y[-2:])
+ ax.plot(x,y, color=c, linestyle='-', linewidth=2, zorder=z)
+
+for i in range(1, len(verify_x)):
+ x,y = m(verify_x[i-1:i+1], verify_y[i-1:i+1])
+ c,z = get_color(verify_speed[i])
+ ax.plot(x,y, color='#FFFFFF', linewidth=4, linestyle="-",zorder=7)
+ ax.plot(x,y, color=c, linewidth=2, linestyle="-",zorder=8)
+
+from matplotlib.lines import Line2D
+l2 = Line2D([], [], linewidth=3, color=get_color(74)[0])
+l3 = Line2D([], [], linewidth=3, color=get_color(95)[0])
+l4 = Line2D([], [], linewidth=3, color=get_color(110)[0])
+
+import matplotlib.font_manager
+prop = matplotlib.font_manager.FontProperties(size=10)
+
+ax.set_title("National Hurricane Center\n Forecast Positions & Max Wind (TCD Product)\n for Isaac (2012)", size=10)
+ax.set_xlabel("Forecasts made between 5 AM 21 August 2012\n and 11 AM 27 August 2012, outlined is observation", size=10)
+ax.legend([l2, l3, l4], ["Tropical Storm", "Category 1", "Category 2"], prop=prop)
+
+fig.savefig('test.png')
+# iemplot.makefeature('test')
View
49 scripts/feature/percentile_monrain.py
@@ -1,49 +0,0 @@
-# Month percentile
-
-import sys, os, random
-sys.path.append("../lib/")
-import iemplot
-
-import mx.DateTime
-now = mx.DateTime.now()
-
-from pyIEM import iemdb, stationTable
-st = stationTable.stationTable("/mesonet/TABLES/coopClimate.stns")
-i = iemdb.iemdb()
-coop = i['coop']
-wepp = i['wepp']
-
-# Extract normals
-lats = []
-lons = []
-vals = []
-rs = coop.query("""SELECT stationid, sum(precip) as rain from alldata
- WHERE month in (6,7) and year = 2010 and
- stationid != 'ia0000' GROUP by stationid""").dictresult()
-for i in range(len(rs)):
- stid = rs[i]['stationid']
- lats.append( st.sts[ rs[i]['stationid'].upper() ]['lat'] )
- lons.append( st.sts[ rs[i]['stationid'].upper() ]['lon'] )
- obs = rs[i]['rain']
- # Find obs
- rs2 = coop.query("""SELECT year, sum(precip) as rain from alldata where
- stationid = '%s' and month IN (6,7) GROUP by year ORDER by rain ASC""" % (stid,)
- ).dictresult()
- for j in range(len(rs2)):
- if rs2[j]['rain'] > obs:
- break
- vals.append( (j+1) / float(len(rs2)) * 100.0 )
-
-cfg = {
- 'wkColorMap': 'BlAqGrYeOrRe',
- 'nglSpreadColorStart': 2,
- 'nglSpreadColorEnd' : -1,
- '_title' : "Iowa Jun/Jul Rainfall Percentile",
- '_valid' : "1 Jun - 31 Jul 2010",
- 'lbTitleString' : "[%]",
- '_showvalues' : True,
- '_format' : '%.0f',
-}
-# Generates tmp.ps
-fp = iemplot.simple_contour(lons, lats, vals, cfg)
-iemplot.postprocess(fp, "")
View
14 scripts/feature/percentile_springhigh.py
@@ -8,13 +8,13 @@
select foo2.station, max(foo2.min) as maxprev, min(foo.d2012) as thisyear,
sum(case when foo.d2012 > foo2.min then 1 else 0 end), count(*) from
- (select station, extract(year from day + '3 months'::interval) as yr,
- min(low) from alldata_ia where month in (10,11,12,1)
- and day < '2011-09-01' and (sday > '0901' or sday < '0111') and station = 'IA2203'
+ (select station, extract(year from day + '0 months'::interval) as yr,
+ avg((high+low)/2.0) as min from alldata_ia where month in (6,7,8)
GROUP by station, yr) as foo2,
- (select station, min(low) as d2012 from alldata_ia where day > '2011-09-01'
- GROUP by station) as foo
+ (select station, avg((high+low)/2.0) as d2012 from alldata_ia where month in (6,7,8)
+ and station != 'IA0000' and substring(station,2,1) != 'C'
+ and year = 2012 GROUP by station) as foo
where foo.station = foo2.station GROUP by foo2.station
""")
@@ -33,8 +33,8 @@
'wkColorMap': 'BlAqGrYeOrRe',
'nglSpreadColorStart': 2,
'nglSpreadColorEnd' : -1,
- '_title' : "2011-2012 Coldest Winter Temperature Percentile",
- '_valid' : "1 Nov 2011 - 9 Jan 2012, comparing 1893-2010 [100% Warmest]",
+ '_title' : "2012 Summer Average Temperature Percentile",
+ '_valid' : "1 Jun 2012 - 31 Aug 2012, comparing 1893-2011 [100% Warmest]",
'lbTitleString' : "[%]",
'_showvalues' : True,
'_format' : '%.0f',
View
61 scripts/feature/percentile_yrhigh.py
@@ -1,61 +0,0 @@
-# Month percentile
-
-import sys, os, random
-sys.path.append("../lib/")
-import iemplot
-
-import mx.DateTime
-now = mx.DateTime.now()
-
-from pyIEM import iemdb, stationTable
-st = stationTable.stationTable("/mesonet/TABLES/coopClimate.stns")
-i = iemdb.iemdb()
-coop = i['coop']
-iem = i['iem']
-mesosite = i['mesosite']
-
-xref = {}
-rs = mesosite.query("SELECT id, climate_site, x(geom) as lon, y(geom) as lat from stations WHERE network in ('IA_ASOS', 'AWOS')").dictresult()
-for i in range(len(rs)):
- xref[ rs[i]['id'] ] = rs[i]
-
-
-# Extract normals
-lats = []
-lons = []
-vals = []
-rs = iem.query("""SELECT station, max(max_tmpf) as high from summary_2010
- WHERE day > '2010-01-01' and network in ('IA_ASOS', 'AWOS') and
- max_tmpf > -30 GROUP by station
- ORDER by high ASC""").dictresult()
-for i in range(len(rs)):
- stid = rs[i]['station']
- lats.append( xref[ rs[i]['station'].upper() ]['lat'] )
- lons.append( xref[ rs[i]['station'].upper() ]['lon'] )
- ob = rs[i]['high']
- cid = xref[ stid ]['climate_site']
- # Find obs
- rs2 = coop.query("""SELECT year, max(high) as highw from alldata where
- stationid = '%s' and sday < '0525' GROUP by year
- ORDER by highw ASC
- """ % (cid.lower(),)
- ).dictresult()
- for j in range(len(rs2)):
- if rs2[j]['highw'] > ob:
- break
- print "[%s] 2010 high: %s rank: %s" % (stid, ob, j)
- vals.append( (j+1) / float(len(rs2)) * 100.0 )
-
-cfg = {
- 'wkColorMap': 'BlAqGrYeOrRe',
- 'nglSpreadColorStart': 2,
- 'nglSpreadColorEnd' : -1,
- '_title' : "2010 Iowa Maximum Temperature Percentile",
- '_valid' : "between 1 Jan - 24 May [100% Warmest]",
- 'lbTitleString' : "[%]",
- '_showvalues' : True,
- '_format' : '%.0f',
-}
-# Generates tmp.ps
-fp = iemplot.simple_contour(lons, lats, vals, cfg)
-iemplot.postprocess(fp, "")
View
52 scripts/feature/percentile_yrrain.py
@@ -1,52 +0,0 @@
-# Month percentile
-
-import sys, os, random
-sys.path.append("../lib/")
-import iemplot
-
-import mx.DateTime
-now = mx.DateTime.now()
-
-from pyIEM import iemdb, stationTable
-st = stationTable.stationTable("/mesonet/TABLES/coopClimate.stns")
-i = iemdb.iemdb()
-coop = i['coop']
-wepp = i['wepp']
-
-# Extract normals
-lats = []
-lons = []
-vals = []
-rs = coop.query("""SELECT stationid, sum(precip) as rain from alldata
- WHERE
- stationid != 'ia0000' and year = 2010 GROUP by stationid ORDER by stationid ASC""").dictresult()
-for i in range(len(rs)):
- stid = rs[i]['stationid']
- if not st.sts.has_key( rs[i]['stationid'].upper() ):
- continue
- lats.append( st.sts[ rs[i]['stationid'].upper() ]['lat'] )
- lons.append( st.sts[ rs[i]['stationid'].upper() ]['lon'] )
- obs = rs[i]['rain']
- # Find obs
- rs2 = coop.query("""SELECT year, sum(precip) as rain from alldata where
- stationid = '%s' and sday < '0818' GROUP by year ORDER by rain ASC""" % (stid,)
- ).dictresult()
- for j in range(len(rs2)):
- if rs2[j]['rain'] > obs:
- break
- vals.append( (j+1) / float(len(rs2)) * 100.0 )
- print rs2[j]['rain'], obs, rs[i]['stationid'], (j+1) / float(len(rs2)) * 100.0
-
-cfg = {
- 'wkColorMap': 'BlAqGrYeOrRe',
- 'nglSpreadColorStart': 2,
- 'nglSpreadColorEnd' : -1,
- '_title' : "Iowa Rainfall Percentile",
- '_valid' : "1 Jan - 18 Aug 2010",
- 'lbTitleString' : "[%]",
- '_showvalues' : True,
- '_format' : '%.0f',
-}
-# Generates tmp.ps
-tmpfp = iemplot.simple_contour(lons, lats, vals, cfg)
-iemplot.postprocess(tmpfp, "")
View
10 scripts/feature/plot_nam.py
@@ -16,7 +16,7 @@
# tot += grib.variables['APCP_P8_L1_GLC0_acc3h'][:]
# grib.close()
-grib = Nio.open_file('nam.t00z.conusnest.hiresf21.tm00.grib2', 'r')
+grib = Nio.open_file('nam.t00z.conusnest.hiresf60.tm00.grib2', 'r')
lats = grib.variables['gridlat_0'][:]
lons = grib.variables['gridlon_0'][:]
#CAPE = grib.variables['CAPE_P0_L1_GLC0'][:]
@@ -33,14 +33,14 @@
cfg = {
'cnLevelSelectionMode' : 'ManualLevels',
- 'cnLevelSpacingF' : 2.,
- 'cnMinLevelValF' : 92,
- 'cnMaxLevelValF' : 112.0,
+ 'cnLevelSpacingF' : 1.0,
+ 'cnMinLevelValF' : 18,
+ 'cnMaxLevelValF' : 45,
'nglSpreadColorEnd': -1,
'nglSpreadColorStart': 2,
# 'cnFillMode' : 'CellFill',
'wkColorMap': 'WhViBlGrYeOrRe',
- '_title' : "NAM 2 meter AGL Air Temperature",
+ '_title' : "NAM 2 meter AGL Air Temperature for Sunday Morning",
'_valid' : "%s Forecast for %s" % (
ts0.strftime("%-I %p %-d %B %Y"), ts.strftime("%-I %p %-d %B %Y")),
'lbTitleString' : "F",
View
61 scripts/feature/q-q.py
@@ -0,0 +1,61 @@
+import iemdb
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+
+years = {}
+ccursor.execute("""
+ select year, rank() over (ORDER by sum ASC), sum,
+ rank() over (ORDER by avg DESC), avg from
+ (select year, sum(precip), avg((high+low)/2.0) from alldata_ia
+ where station = 'IA0000' and month in (6,7,8) and year < 2012 GROUP by year) as foo
+""")
+for row in ccursor:
+ years[row[0]] = {'prank': row[1], 'precip': row[2],
+ 'trank': row[3], 'avgt': row[4]}
+
+ccursor.execute("""
+ select year, rank() over (ORDER by sum ASC), sum,
+ rank() over (ORDER by avg DESC), avg from
+ (select year, sum(precip), avg((high+low)/2.0) from alldata_ia
+ where station = 'IA0000' and month in (9,10,11) and year < 2012 GROUP by year) as foo
+""")
+for row in ccursor:
+ years[row[0]]['fall_prank'] = row[1]
+ years[row[0]]['fall_precip'] = row[2]
+ years[row[0]]['fall_trank'] = row[3]
+ years[row[0]]['fall_avgt'] = row[4]
+
+ccursor.execute("""
+ select year2, rank() over (ORDER by sum ASC), sum,
+ rank() over (ORDER by avg DESC), avg from
+ (select extract(year from day - '3 months'::interval) as year2, sum(precip), avg((high+low)/2.0) from alldata_ia
+ where station = 'IA0000' and month in (12,1,2) and day >= '1893-11-01' GROUP by year2) as foo
+""")
+for row in ccursor:
+ years[row[0]]['winter_prank'] = row[1]
+ years[row[0]]['winter_precip'] = row[2]
+ years[row[0]]['winter_trank'] = row[3]
+ years[row[0]]['winter_avgt'] = row[4]
+
+one = []
+two = []
+for year in range(1893,2012):
+ print '%s,%s,%.2f,%s,%.1f,%s,%.2f,%s,%.1f,%s,%.2f,%s,%.1f' % (year,
+ years[year]['prank'], years[year]['precip'],
+ years[year]['trank'], years[year]['avgt'],
+ years[year]['fall_prank'], years[year]['fall_precip'],
+ years[year]['fall_trank'], years[year]['fall_avgt'],
+ years[year]['winter_prank'], years[year]['winter_precip'],
+ years[year]['winter_trank'], years[year]['winter_avgt']
+ )
+
+ one.append(years[year]['prank'] )
+ two.append(years[year]['fall_trank'] )
+
+import matplotlib.pyplot as plt
+
+(fig, ax) = plt.subplots(1,1)
+
+ax.scatter(one, two)
+
+fig.savefig('test.png')
View
29 scripts/feature/weather_feels_like.py
@@ -0,0 +1,29 @@
+import iemdb
+COOP = iemdb.connect('coop', bypass=True)
+ccursor = COOP.cursor()
+import mx.DateTime
+
+sts = mx.DateTime.DateTime(2012,1,1)
+ets = mx.DateTime.now()
+interval = mx.DateTime.RelativeDateTime(days=14)
+now = sts
+
+while now < ets:
+ ccursor.execute("""select station,
+ SQRT( sum((dhigh - chigh)^2 + (dlow - clow)^2)/count(*) ) as diff
+ from (select c.station, c.valid, d.high as dhigh, c.high as chigh,
+ d.low as dlow, c.low as clow from ncdc_climate71 c JOIN alldata_ia d
+ on (d.sday = to_char(c.valid, 'mmdd')) WHERE d.station = 'IA2203'
+ and d.day BETWEEN '%s'::date and '%s'::date + '14 days'::interval) as foo
+ GROUP by station ORDER by diff ASC LIMIT 1""" % (now.strftime("%Y-%m-%d"),
+ now.strftime("%Y-%m-%d")))
+ row = ccursor.fetchone()
+ ccursor.execute("""SELECT name, state from stations where
+ id = '%s'""" % (row[0].upper(),))
+ row2 = ccursor.fetchone()
+
+ print '%s-%s %s, %s' % (now.strftime("%d %b"),
+ (now + mx.DateTime.RelativeDateTime(days=14)).strftime("%d %b"),
+ row2[0].capitalize(), row2[1])
+
+ now += interval
View
44 scripts/feature/xy2.py
@@ -2,25 +2,27 @@
from scipy import stats
import iemdb
-COOP = iemdb.connect('coop', bypass=True)
+COOP = iemdb.connect('asos', bypass=True)
ccursor = COOP.cursor()
#select foo.year, foo.april, foo2.march from (select year, avg(high) as april from alldata where stationid = 'ia2203' and month = 4 and sday < '0408' GROUP by year) as foo, (select year, avg(high) as march from alldata where stationid = 'ia2203' and sday > '0323' and sday < '0401' GROUP by year) as foo2 where foo.year = foo2.year
ccursor.execute("""
-select year, sum( case when sday < '0523' and high > 69 then 1 else 0 end), sum(case when sday > '0523' and high > 89 then 1 else 0 end ) from alldata where stationid = 'ia0200' and year < 2011 GROUP by year ORDER by year DESC
+ select extract(year from valid) as yr,
+ max(case when extract(month from valid) = 8 then dwpf else 0 end) -
+ max(case when extract(month from valid) = 9 then dwpf else 0 end) as d,
+ max(case when extract(month from valid) = 8 then dwpf else 0 end) as a,
+ max(case when extract(month from valid) = 9 then dwpf else 0 end) from alldata
+ where station = 'DSM' GROUP by yr ORDER by yr ASC
""")
x = []
y = []
-y2 = []
-x2 = []
+years = []
for row in ccursor:
- if row[0] < 1990:
- x2.append( float(row[1]) )
- y2.append( float(row[2]) )
- else:
- x.append( float(row[1]) )
- y.append( float(row[2]) )
+ if row[2] > 10:
+ x.append( float(row[2]) )
+ y.append( float(row[3]) )
+ years.append( row[0] )
from matplotlib import pyplot as plt
@@ -30,25 +32,27 @@
fig = plt.figure()
ax = fig.add_subplot(111)
-ax.scatter(x, y, label="After 1990")
-ax.scatter(x2, y2, color='r', label="Prior 1990")
+ax.scatter(x, y)
+ax.scatter(x[-1],y[-1], facecolor='r', label='2012')
+for i in range(len(years)):
+ if y[i] > x[i]:
+ ax.text(x[i]-0.4, y[i], "%.0f" % (years[i],), ha='right')
#ax.scatter(18, -4.5, color='g', marker='+', s=100)
-ax.text(15, 55, "2011", ha='center', va='center', color='b')
-ax.plot( [15,15], [0,70], color='b')
+ax.plot( [60,90], [60,90], color='k')
#ax.text(26, 56, "Average", ha='center', va='center', color='g')
#ax.plot( [20,80], [58,58], color='g')
#ax.plot( [0,31], [intercept, intercept + 31 * h_slope], color='r',
# label=r"Fit = $\frac{%.2f ^{\circ}\mathrm{F}}{day}, R^2 = %.2f$" % (
# h_slope, h_r_value ** 2))
-ax.set_title("Ames Daily High Temperatures [1893-2010]")
-ax.set_xlabel("Days prior to 23 May above 70 F")
-ax.set_ylabel("Days after 23 May above 90 F")
-ax.set_xlim(-0.5,50)
-ax.set_ylim(-0.5,70)
+ax.set_title("Des Moines Monthly Maximum Dew Point [1933-3 Sep 2012]")
+ax.set_xlabel("Maximum August Dew Point $^{\circ}\mathrm{F}$")
+ax.set_ylabel("Maximum September Dew Point $^{\circ}\mathrm{F}$")
+ax.set_xlim(60,90)
+ax.set_ylim(60,90)
#ax.set_yticks( range(0,30,6))
#ax.set_xticks( range(0,60,6))
-ax.legend()
+#ax.legend(loc='best')
ax.grid(True)
import iemplot
fig.savefig('test.ps')
View
39 scripts/feature/yearly_precip_range.py
@@ -1,19 +1,19 @@
import iemdb, network
import numpy
-COOP = iemdb.connect('postgis', bypass=True)
+COOP = iemdb.connect('coop', bypass=True)
ccursor = COOP.cursor()
ccursor.execute("""
-select extract(year from date) as yr, sum(case when extract(doy from date) < 232 then 1 else 0 end), sum(case when extract(doy from date) > 231 then 1 else 0 end) from (select distinct date(issue - '7 hours'::interval) from warnings where phenomena = 'TO' and ugc ~* 'IAC' and significance = 'W') as foo GROUP by yr ORDER by yr ASC
+select year, min(extract(doy from day)) from alldata_ia where station = 'IA2203' and high < 70 and month > 6 GROUP by year ORDER by year
+
""")
years = []
-prior = []
-total = []
+doy = []
for row in ccursor:
years.append( row[0] )
- prior.append( float(row[1]) )
- total.append( float(row[1] + row[2]) )
+ doy.append( float(row[1]) )
+# total.append( float(row[1] + row[2]) )
years = numpy.array(years)
@@ -23,19 +23,22 @@
prop = matplotlib.font_manager.FontProperties(size=12)
fig, ax = plt.subplots(1,1)
-bars = ax.bar( years - 0.4, total,
- facecolor='g', ec='g', zorder=1, label='After 19 Aug')
-bars = ax.bar( years - 0.4, prior,
- facecolor='lightblue', ec='lightblue', zorder=2, label='Prior 19 Aug')
-for i in range(len(total)):
- ax.text(1986+i, total[i]+0.5, "%.0f" % (total[i],), ha='center')
-
-ax.set_title("Days with Tornado Warning Issued in Iowa (1986-2012)")
+bars = ax.bar( years - 0.4, doy,
+ facecolor='purple', ec='purple', zorder=1, label='After 19 Aug')
+for bar in bars:
+ if bar.get_height() > bars[-1].get_height():
+ bar.set_facecolor('r')
+ bar.set_edgecolor('r')
+
+ax.set_title("First Day with High Temperature below 70$^{\circ}\mathrm{F}$\nfor Des Moines after 1 July (1880-2012)")
ax.grid(True)
-ax.set_xlabel("thru 19 August 2012")
-ax.set_ylabel('Number of Days (7 AM - 7 AM)')
-ax.set_xlim(1985.5,2013.5)
-ax.legend(loc='best')
+ax.set_ylabel('First Date')
+ax.set_xlim(1889.5,2013.5)
+ax.set_xlabel("years in red were later than 2012")
+ax.set_yticks( (1,32,60,91,121,152,182,188,195,202,213,220,227,233,244,251,258,265,274,305,335,365) )
+ax.set_yticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul 1','Jul 8', 'Jul 15', 'Jul 22', 'Aug 1','Aug 8', 'Aug 15', 'Aug 22', 'Sep 1','Sep 8', 'Sep 15', 'Sep 22', 'Oct','Nov','Dec') )
+ax.set_ylim(182,275)
+#ax.legend(loc='best')
fig.savefig('test.ps')
View
62 scripts/feature/yearly_range.py
@@ -1,43 +1,57 @@
import iemdb, network
import numpy
import numpy.ma
+import mx.DateTime
-COOP = iemdb.connect('asos', bypass=True)
+COOP = iemdb.connect('coop', bypass=True)
ccursor = COOP.cursor()
#select extract(year from o.day + '2 months'::interval) as yr, sum(case when o.high > c.high then 1 else 0 end), count(*) from alldata_ia o , climate51 c WHERE c.station = o.station and o.station = 'IA0200' and extract(month from c.valid) = extract(month from o.day) and extract(day from c.valid) = extract(day from o.day) and (o.month in (11,12,1) or (o.month = 2 and extract(day from o.day) < 9)) and o.day > '1893-06-01' GROUP by yr ORDER by yr ASC
ccursor.execute("""
-select doy, year, count(*) from (select distinct extract(year from valid) as year, extract(hour from valid) as hour, extract(doy from valid) as doy from alldata where station = 'DSM' and p01i > 0) as foo GROUP by doy, year
+ select year, max(case when high > 89 then day else '1880-01-01'::date end),
+ avg((high+low)/2.0) from alldata_ia
+ where station = 'IA2203' and year < 2012 and year > 1879 GROUP by year ORDER by year ASC
""")
-doy = []
-count = []
-data = numpy.ma.zeros((24,366), 'f')
-
+years = []
+last = []
+total = []
for row in ccursor:
- data[:int(row[2])-1,int(row[0])-1] += 1
+ ts = mx.DateTime.strptime("%s%s%s" % (row[1].year, row[1].month, row[1].day), "%Y%m%d")
+ years.append( row[0] )
+ last.append( int(ts.strftime("%j")) )
+ total.append( float(row[2]) )
-data.mask = numpy.where(data==0,True,False)
+years = numpy.array( years )
+last = numpy.array( last )
import matplotlib.pyplot as plt
import iemplot
-fig, ax = plt.subplots(1, 1)
-
-#ax.scatter( doy, count)
-res= ax.imshow(data,aspect='auto', interpolation='nearest', extent=[1,366,0.5,24.5],
- origin='lower')
-fig.colorbar(res)
-
-ax.set_yticks( numpy.arange(1,25,4))
-ax.set_ylabel("Hours")
-
-ax.set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
-ax.set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
-ax.set_xlim(0, 366)
-ax.set_title("Des Moines Hours with Precipitation Reported [1933-2012]\nFrequency of Number of Hours per Day of Year")
-ax.set_ylim(0.5,24.5)
-ax.grid(True)
+fig, ax = plt.subplots(2, 1)
+
+bars = ax[0].bar(years-0.5, last, fc='pink', ec='pink')
+for bar in bars:
+ if bar.get_height() < numpy.average(last):
+ bar.set_edgecolor('skyblue')
+ bar.set_facecolor('skyblue')
+ax[0].set_yticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax[0].set_yticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax[0].set_ylim(150,300)
+ax[0].set_ylabel("Date of last 90+$^{\circ}\mathrm{F}$")
+ax[0].set_title("1880-2011 Des Moines Last Date Above 90$^{\circ}\mathrm{F}$")
+ax[0].set_xlim(1879.5,2012)
+ax[0].grid(True)
+ax[0].plot([1880,2011], [numpy.average(last), numpy.average(last)], color='k')
+
+ax[1].scatter(total, last)
+#ax[1].set_ylim(30,80)
+ax[1].set_yticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax[1].set_yticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax[1].set_ylim(150,300)
+ax[1].set_xlabel("Average Temperature for Year $^{\circ}\mathrm{F}$")
+ax[1].set_ylabel("Date of last 90+$^{\circ}\mathrm{F}$")
+ax[1].grid(True)
fig.savefig('test.ps')
View
62 scripts/feature/ytd_lines.py
@@ -14,40 +14,56 @@
# Save 1931 (record year)
obs1931 = obs[1931-1880,:].copy()
-
-# Now, we replace each year's 1 Jan - 18 Jul with 2012
-obs[:,0:200] = obs[-1,0:200]
-
-ytd = numpy.zeros((2013-1880, 366), 'f')
+obs2 = obs.copy()
ytd1931 = numpy.zeros((366,), 'f')
for day in range(366):
- ytd[:,day] = numpy.average(obs[:,:(day+1)],1)
ytd1931[day] = numpy.average(obs1931[:(day+1)])
-count = 0
-total = 2012-1880
-for year in range(1880,2012):
- if ytd[year-1880,-2] > ytd1931[-2]:
- count += 1
+
+freq = numpy.zeros((366), 'f')
+for dy in range(60,263):
+ # Now, we replace each year's 1 Jan - 18 Jul with 2012
+ obs2[:,0:dy] = obs[-1,0:dy]
+
+ ytd = numpy.zeros((2013-1880, 366), 'f')
+
+ for day in range(366):
+ ytd[:,day] = numpy.average(obs2[:,:(day+1)],1)
+
+ count = 0
+ total = 2012-1880
+ for year in range(1880,2012):
+ if ytd[year-1880,-2] > ytd1931[-2]:
+ count += 1
+ freq[dy] = float(count) / float(total) * 100.0
import matplotlib.pyplot as plt
-(fig, ax) = plt.subplots(1,1)
+(fig, ax) = plt.subplots(2,1, sharex=True)
-for year in range(1880,2012):
- ax.plot(range(366), ytd[year-1880,:], c='tan')
+for year in range(1880,2013):
+ if year == 2012:
+ ax[0].plot(range(263), ytd[year-1880,:263], c='red', label='2012')
+ else:
+ ax[0].plot(range(366), ytd[year-1880,:], c='tan')
-ax.plot(range(366), ytd1931[:], c='r', label='1931')
-ax.set_ylabel("YTD Average Temperature $^{\circ}\mathrm{F}$")
-ax.set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
-ax.set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
-ax.grid(True)
-ax.set_title("Chance of 2012 being Warmest Year on Record for Des Moines\nScenarios: append previous year time series to this year's YTD value\n%s/%s (%.1f%%) years (1880-2011) finish warmer than 1931" % (count, total,
+ax[0].plot(range(366), ytd1931[:], c='b', label='1931')
+ax[0].set_ylabel("YTD Average Temperature $^{\circ}\mathrm{F}$")
+ax[0].set_xticks( (1,32,60,91,121,152,182,213,244,274,305,335,365) )
+ax[0].set_xticklabels( ('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )
+ax[0].grid(True)
+ax[0].set_title("Chance of 2012 being Warmest Year on Record for Des Moines\nScenarios: append previous year time series to this year's YTD value\n%s/%s (%.1f%%) years (1880-2011) finish warmer than 1931" % (count, total,
float(count) / float(total) * 100.0))
-ax.set_xlim(200,364)
-ax.set_ylim(50,65)
-ax.legend()
+ax[0].set_xlim(60,364)
+ax[0].set_ylim(20,65)
+ax[0].legend(loc=4)
+
+ax[1].plot(range(263), freq[:263], c='r', label='1931')
+ax[1].grid(True)
+ax[1].set_ylabel("Scenario Probability [%]")
+ax[1].set_ylim(0,100)
+
plt.tight_layout()
fig.savefig('test.ps')
import iemplot

0 comments on commit fd51a61

Please sign in to comment.