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.
Initilization of the graphics system and computational modules used in this IPython notebook.
# 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/'
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)
!convert $fname -trim $fname
!convert $fname -transparent white $fname
Load the rule curves from another notebook.
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)
!convert $fname -trim $fname
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)
!convert $fname -trim $fname
!convert $fname -transparent white $fname
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)
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()
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)
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
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');
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()
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
Load the rule curves from another notebook.
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)
!convert $fname -trim $fname