text_html_output_1.png
from past.builtins import xrange

import pandas as pd

import seaborn as sns

webpth = 'http://www.files.benlaken.com/documents/'
monsoon = pd.read_csv('../input/Monsoon_data.csv', parse_dates=['Date'])
monsoon.index = monsoon.Date
monsoon = monsoon.drop('Date', 1)
olou = pd.read_csv('../input/Olou_counts.csv', parse_dates=['Date'])
olou.index = olou.Date
olou = olou.drop('Date', 1)

# Plot the simple time series 
my_ts = plt.figure()
my_ts.set_size_inches(10,5)         # Specify the output size
ax1 = my_ts.add_subplot(211)        # Add an axis frame object to the plot (i.e. a pannel)
ax2 = my_ts.add_subplot(212) 

ax1.step(monsoon.index.date,monsoon.Precip,lw=1.0)
ax1.set_title(r'Monthly Precipitation and NM counts')
ax1.set_ylabel(r'Precipitation (mm)')
ax1.grid(True)
#ax1.set_yscale('log')

ax2.plot(olou.index.date,olou.Counts/1000,'r.',ms=3.0)
ax2.set_ylabel(r'Olou NM (cnt./min.$\times10^{3}$)')
ax2.set_xlabel('Date')
ax2.grid(True)
plt.show(my_ts)
my_ts.savefig('Monthly_ts.png',dpi=300)

def return_stderr(data):
    """Calculate uncertainty of a np array as Standard Error of the Mean"""
    return np.nanstd(data) / np.sqrt(np.count_nonzero(data) - 1)
climo = {}
climo['means'] = [np.mean(monsoon.Precip[monsoon.index.month == mnth + 1]) for mnth in xrange(12)]
climo['error'] = [return_stderr(monsoon.Precip[monsoon.index.month == mnth + 1].values) for mnth in xrange(12)]

# -- Plot the climatology --
my_climo = plt.figure()
my_climo.set_size_inches(5,5)
ax1 = my_climo.add_subplot(111)
ax1.errorbar(x=range(12),y=climo['means'],yerr=climo['error'])
ax1.set_title(r'Precipitation climatology')
ax1.set_ylabel(r'Precipitation (mm)')
ax1.set_xlabel(r'Month')
ax1.set_xlim(0,11)
ax1.set_xticklabels(labels=['Jan','Mar','May','Jul','Sep','Nov'])
ax1.grid(True)
plt.show(my_climo)
my_climo.savefig('Monthly_climo.png',dpi=300)

delta = []
for date in monsoon.Precip.index:
    delta.append(monsoon.Precip[date] - climo['means'][date.month - 1])
dseries = pd.Series(delta, index=monsoon.index)

def lookup_index(yr):
    return (monsoon.index.year == yr) & (monsoon.index.month >= 5) & (monsoon.index.month <= 9)
mjjas = {}
mjjas['means'] = [np.mean(dseries[lookup_index(yr)]) for yr in xrange(1964, 2012, 1)]
mjjas['SEM'] = [return_stderr(dseries[lookup_index(yr)]) for yr in xrange(1964, 2012, 1)]
mjjas['sum'] = [np.sum(dseries[lookup_index(yr)]) for yr in xrange(1964, 2012, 1)]

sns.set(style='darkgrid')
yrange = xrange(1964, 2012, 1)
my_mjjas = plt.figure()
my_mjjas.set_size_inches(10, 5)
ax1 = my_mjjas.add_subplot(121)
ax2 = my_mjjas.add_subplot(122)
ax1.errorbar(x=yrange, y=mjjas['means'], yerr=mjjas['SEM'], fmt='.', ms=10)
ax1.set_xlim(min(yrange) - 1, max(yrange) + 1)
ax1.set_title('Mean MJJAS precipitation anomaly')
ax1.set_ylabel('$\\delta$ precipitation (mm/month)')
ax1.set_xlabel('Year')
ax1.grid(True)
sns.distplot(mjjas['means'], ax=ax2)
ax2.set_title('Distribution of MJJAS anomalies')
ax2.set_xlabel('$\\delta$ precipitation (mm/month)')
ax2.set_ylabel('Density')
plt.show(my_mjjas)
my_mjjas.savefig('delta_precip_pop.png', dpi=300)