#!/usr/bin/env python
# coding: utf-8

# # Namakan and Rainy Lake Water Levels 1970-2014

# [Rainy Lake](http://en.wikipedia.org/wiki/Rainy_Lake) is a relatively large lake on Minnesota/Ontario border that makes a portion of [Voyageurs National Park](http://en.wikipedia.org/wiki/Voyageurs_National_Park). In recent years Rainy Lake has experienced periods of high water levels and flooding. The purpose of these notes and calculations here is to quantify those episodes using publically available data sets, explore possible reasons for these events, and what can be done to reduce the risk of future high water and flooding events.

# ### Initialization

# Initilization of the graphics system and computational modules used in this IPython notebook.

# In[1]:


# Display graphics inline with the notebook
get_ipython().run_line_magic('matplotlib', 'inline')

# Standard Python modules
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import datetime

# Modules to display images and data tables
from IPython.display import Image
from IPython.core.display import display

import seaborn as sns
sns.set_context('talk')

# Data directory
dir = '../data/'
img = '../images/'


# ## Namakan Lake and Rainy Lake Levels

# In[71]:


plt.figure(figsize=(12,6))

RL = pd.read_pickle(dir+'RL.pkl')
NL = pd.read_pickle(dir+'NL.pkl')

RL.plot(lw=1.5)
NL.plot(lw=1.5)

ax1 = plt.gca()
x1,x2 = ax1.get_xlim()
y1,y2 = ax1.get_ylim()

def rcPlot(yr,mo,dy):
    t = datetime.datetime(yr,mo,dy).toordinal() - datetime.datetime(1969,12,31).toordinal()
    plt.plot([t,t],[336,y2],'r--',lw=1.5)
    plt.text(t,335.7,'{0:4d}'.format(yr),ha='center',va='top',size=14)

rcPlot(1940,10,3)
rcPlot(1949,6,8)
rcPlot(1957,10,1)
rcPlot(1970,7,29)
rcPlot(2000,1,5)

plt.legend([RL.name,NL.name]);
plt.title('Lake Levels' + ' ' 
         + str(RL.dropna().index[0].year) + '-' + str(RL.dropna().index[-1].year))
plt.ylabel('meters')

m2f = 3.28084

ax2 = plt.twinx()
ax2.set_xlim(x1,x2)
ax2.set_yticks([m2f*y for y in ax1.get_yticks()])
ax2.set_ylim(m2f*y1,m2f*y2)
ax2.set_ylabel('feet')
ax2.grid(None)

fname = img + 'RainyNamakanLakeLevels.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')
get_ipython().system('convert $fname -transparent white $fname')


# ## Rule Curve Performance 1970-1999

# Load the rule curves from another notebook.

# In[72]:


NL1970 = pd.read_pickle(dir+'NL1970.pkl')
RL1970 = pd.read_pickle(dir+'RL1970.pkl')

plt.figure(figsize=(10,7))

NL1970['LRC'].plot(color='y')
NL1970['URC'].plot(color='y')

RL1970['LRC'].plot(color='y')
RL1970['URC'].plot(color='y')

plt.fill_between(NL1970.index,
                 NL1970['LRC'].tolist(),
                 NL1970['URC'].tolist(),
                 color='y',alpha=1)

plt.fill_between(NL1970.index,
                 NL1970['LRC'].tolist(),
                 NL1970['URC'].tolist(),
                 color='y', alpha=1)

plt.fill_between(RL1970.index,
                 RL1970['LRC'].tolist(),
                 RL1970['URC'].tolist(),
                 color='y', alpha=1)

plt.ylabel('Lake Level [meters]')
plt.title('Rule Curves for Namakan and Rainy Lakes: 1971-1999')

e = pd.Series([])

for (yr,r) in RL['1971':'1999'].groupby(RL['1971':'1999'].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r = r.tshift(shift.days)
    r.plot()
    e = e.append((r - RL1970['URC']).tshift(-shift.days))

for (yr,r) in NL['1971':'1999'].groupby(NL['1971':'1999'].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r.tshift(shift.days).plot()
    
import matplotlib.dates as mdates

plt.gca().xaxis.set_major_locator(mdates.MonthLocator())
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b'))

fname = img + 'RuleCurvePerformance1970-1999.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')


# In[73]:


NL2000 = pd.read_pickle(dir+'NL2000.pkl')
RL2000 = pd.read_pickle(dir+'RL2000.pkl')

plt.figure(figsize=(10,7))

NL2000['LRC'].plot(color='b')
NL2000['URC'].plot(color='b')

RL2000['LRC'].plot(color='b')
RL2000['URC'].plot(color='b')

plt.fill_between(NL2000.index,
                 NL2000['LRC'].tolist(),
                 NL2000['URC'].tolist(),
                 color='b', alpha=1)

plt.fill_between(RL2000.index,
                 RL2000['LRC'].tolist(),
                 RL2000['URC'].tolist(),
                 color='b', alpha=1)

plt.ylabel('Lake Level [meters]')
plt.title('Rule Curves for Namakan and Rainy Lakes: 2000-2010')

for (yr,r) in RL['2000':'2010'].groupby(RL['2000':'2010'].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r = r.tshift(shift.days)
    r.plot()
    e = e.append((r - RL2000['URC']).tshift(-shift.days))
    
for (yr,r) in NL['2000':'2010'].groupby(NL['2000':'2010'].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r.tshift(shift.days).plot()

import matplotlib.dates as mdates

plt.gca().xaxis.set_major_locator(mdates.MonthLocator())
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b'))

fname = img + 'RuleCurvePerformance2000-2010.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')
get_ipython().system('convert $fname -transparent white $fname')


# ## Frequency and Distribution of High Water Events

# In[75]:


p1970 = RL['1971':'1999'].index
pA = p1970[(p1970.month >= 5) & (p1970.month <= 9)]

p2000 = RL['2000':'2010'].index
pB = p2000[(p2000.month >= 5) & (p2000.month <= 9)]

def highEvents(str,S,h):
    print("\n{:s}".format(str),end="")
    print("exceeded {:d} days out of {:d} days.".format(S[S>h].count(),len(S.index)))
    print("{:>20s}: {:6.2f}%".format("Frequency", 100.0*float(S[S>h].count())/len(S.index)))
    print("{:>20s}: {:7.3f} meters".format("Median Value", (S[S>h]-h).median()))
    print("{:>20s}: {:7.3f} meters".format("95th Percentile", (S[S>h]-h).quantile(0.95)))
    r = pd.Series([])
    r['Days in Period'] = len(S.index)
    r['Days exceeded'] = S[S>h].count()
    r['Frequency (%)'] = 100.0*float(S[S>h].count())/len(S.index)
    r['Median'] = (S[S>h]-h).median()
    r['95th Percentile'] = (S[S>h]-h).quantile(0.95)
    return r
    
def highWater(per):
    print("\n\nPeriod: ", per[0].date().strftime("%Y/%b/%d"), "-", per[-1].date().strftime("%Y/%b/%d"))
    r = highEvents("Rule Curve",e[per],0.0)
    r = highEvents("Emergency High Water",RL[per],337.75)
    r = highEvents("All Gates Open",RL[per],337.90)
    
highWater(pA)
highWater(pB)


# In[76]:


plt.subplot(2,1,1)
plt.hist(e[pA],bins=np.arange(0,1.0,.02),color='y',alpha=0.5,normed=True)
plt.xlim(0,1.0)
plt.ylim(0,12)
plt.title("{0} - {1}".format(pA[0].date().strftime("%Y/%b/%d"),pA[-1].date().strftime("%Y/%b/%d")))

plt.subplot(2,1,2)
plt.hist(e[pB],bins=np.arange(0,1.0,.02),color='b',alpha=0.5,normed=True)
plt.xlim(0,1.0)
plt.ylim(0,12)
plt.title("{0} - {1}".format(pB[0].date().strftime("%Y/%b/%d"),pB[-1].date().strftime("%Y/%b/%d")))

plt.tight_layout()


# In[79]:


plt.figure(figsize=(9,5))

pd.Series(e[pA]).hist(alpha=0.4,bins=np.arange(0,1,.05),color='y',normed=True)
pd.Series(e[pB]).hist(alpha=0.4,bins=np.arange(0,1,.05),color='b',normed=True)


# ## Stage Frequency for Rainy Lake Levels

# In[80]:


xlo = 336.3
xhi = 338.6
dx = 0.02

plt.figure(figsize=(10,4))
ts = RL['1971':'1999']
ts.hist(bins=np.arange(xlo,xhi+dx,dx))
ax = plt.axis()
plt.axis([xlo,xhi,ax[2],ax[3]])
plt.title('Distribution of Rainy Lake Levels, 1971-1999')

fname = img + 'RainyLakeLevelDistribution1970.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')
get_ipython().system('convert $fname -transparent white $fname')

plt.figure(figsize=(12,4))
ts = RL['2000':'2012']
ts.hist(bins=np.arange(xlo,xhi+dx,dx))
ax = plt.axis()
plt.axis([xlo,xhi,ax[2],ax[3]])
plt.xlabel('Lake Level [meters]')
plt.title('Distribution of Rainy Lake Levels, 2000-2012');

fname = img + 'RainyLakeLevelDistribution2000.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')


# In[83]:


plt.figure(figsize=(10,4))

RL['1970':'1999'].hist(cumulative=True,
        bins=np.arange(xlo,xhi+dx,dx),
        alpha=0.4,
        color = 'b',
        normed=True)

RL['2000':'2012'].hist(cumulative=True,
        bins=np.arange(xlo,xhi+dx,dx),
        alpha=0.4,
        color = 'y',
        normed=True)

plt.axis([xlo,xhi,0,1.1])
plt.xlabel('Lake Level [meters]')

plt.legend(['1970-1999','2000-2012'],loc = 'upper left')
plt.title('Cumulative Distribution of Daily Levels for Rainy Lake');


# In[84]:


tsa = RL['1970':'1999']
tsb = RL['2000':'2012']
mostr = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

plt.figure(figsize=(10,8))
for mo in range(1,13):
    plt.subplot(4,3,mo)
    tsa[tsa.index.month == mo].hist(cumulative=True,
                                    normed=True,
                                    bins=np.arange(xlo,xhi+3*dx,3*dx),
                                    color = 'b',
                                    alpha = 0.4)
    tsb[tsb.index.month == mo].hist(cumulative=True,
                                    normed=True,
                                    bins=np.arange(xlo,xhi+3*dx,3*dx),
                                    color = 'y',
                                    alpha = 0.4)
    plt.legend(['1970-1999','2000-2012'],loc='upper left')
    plt.axis([xlo,xhi,0,1.05])
    plt.title(mostr[mo-1])
    
plt.tight_layout()


# In[85]:


plt.figure()

hist,bins = np.histogram([t for t in RL['1970':'1999'] if pd.notnull(t)],bins=100)
hist = np.cumsum(hist)
plt.semilogx([1-float(h)/hist.max()  for h in hist],bins[:-1],color='b')

hist,bins = np.histogram([t for t in RL['2000':'2012'] if pd.notnull(t)],bins=100)
hist = np.cumsum(hist)
plt.semilogx([1-float(h)/hist.max()  for h in hist],bins[:-1],color='r')

bins.size
hist.size


# In[ ]:





# ## Rule Curve Performance 2000-2014

# Load the rule curves from another notebook.

# In[86]:


NL2000 = pd.read_pickle(dir+'NL2000.pkl')
RL2000 = pd.read_pickle(dir+'RL2000.pkl')

plt.figure(figsize=(12,7))

NL2000['LRC'].plot(color='y')
NL2000['URC'].plot(color='y')

RL2000['LRC'].plot(color='y')
RL2000['URC'].plot(color='y')

NL2000['ELW'].plot(color='b',alpha=0.5)
NL2000['EDL'].plot(color='b',alpha=0.5)
NL2000['EHW'].plot(color='b',alpha=0.5)
NL2000['AGO'].plot(color='b',alpha=0.5)

RL2000['ELW'].plot(color='r',alpha=0.5)
RL2000['EDL'].plot(color='r',alpha=0.5)
RL2000['EHW'].plot(color='r',alpha=0.5)
RL2000['AGO'].plot(color='r',alpha=0.5)

plt.fill_between(NL2000.index,
                 NL2000['LRC'].tolist(),
                 NL2000['URC'].tolist(),
                 color='b', alpha=0.5)

plt.fill_between(RL2000.index,
                 RL2000['LRC'].tolist(),
                 RL2000['URC'].tolist(),
                 color='r', alpha=0.5)

plt.ylabel('Lake Level [meters]')
plt.title('Rule Curves for Namakan and Rainy Lakes: ' 
        + '2000 - ' + str(RL.dropna().index[-1].year))

e = pd.Series([])

for (yr,r) in RL['2000':].groupby(RL['2000':].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r = r.tshift(shift.days)
    r.plot()
    e = e.append((r - RL2000['URC']).tshift(-shift.days))

for (yr,r) in NL['2000':].groupby(NL['2000':].index.year):
    shift = datetime.datetime(2014,1,1) - datetime.datetime(yr,1,1)
    r.tshift(shift.days).plot()
    
plt.ylim(336,341.5)

import matplotlib.dates as mdates

plt.gca().xaxis.set_major_locator(mdates.MonthLocator())
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b'))

fname = img + 'RuleCurvePerformance2000-2014.png'

plt.savefig(fname)
get_ipython().system('convert $fname -trim $fname')


# In[ ]: