# Namakan and Rainy Lake Water Levels 1970-2014¶

Rainy Lake is a relatively large lake on Minnesota/Ontario border that makes a portion of 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
%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.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)
!convert $fname -trim$fname
!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')

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)
!convert $fname -trim$fname

In [73]:
NL2000 = pd.read_pickle(dir+'NL2000.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)
!convert $fname -trim$fname
!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)


Period:  1971/May/01 - 1999/Sep/30

Rule Curveexceeded 661 days out of 4437 days.
Frequency:  14.90%
Median Value:   0.068 meters
95th Percentile:   0.383 meters

Emergency High Waterexceeded 346 days out of 4437 days.
Frequency:   7.80%
Median Value:   0.048 meters
95th Percentile:   0.356 meters

All Gates Openexceeded 84 days out of 4437 days.
Frequency:   1.89%
Median Value:   0.124 meters
95th Percentile:   0.288 meters

Period:  2000/May/01 - 2010/Sep/30

Rule Curveexceeded 303 days out of 1683 days.
Frequency:  18.00%
Median Value:   0.230 meters
95th Percentile:   0.704 meters

Emergency High Waterexceeded 235 days out of 1683 days.
Frequency:  13.96%
Median Value:   0.187 meters
95th Percentile:   0.710 meters

All Gates Openexceeded 147 days out of 1683 days.
Frequency:   8.73%
Median Value:   0.171 meters
95th Percentile:   0.621 meters

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)

Out[79]:
<matplotlib.axes._subplots.AxesSubplot at 0x118ceaf28>

## 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)
!convert $fname -trim$fname
!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)
!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

Out[85]:
100

## Rule Curve Performance 2000-2014¶

Load the rule curves from another notebook.

In [86]:
NL2000 = pd.read_pickle(dir+'NL2000.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)
!convert $fname -trim$fname