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