In [1]:
%pylab inline
import os
import sys
import math
import datetime
import pandas as pd
import scipy.signal as signl
import scipy.linalg as splin
import scipy.optimize as sciopt
from matplotlib import dates as mdates
#--default matplotlib
rcParams['mathtext.default'] = 'regular'
rcParams['legend.fontsize'] = 6
rcParams['axes.labelsize'] = 6
rcParams['xtick.labelsize'] = 6
rcParams['ytick.labelsize'] = 6
rcParams['figure.subplot.wspace'] = 0.075
rcParams['figure.subplot.hspace'] = 0.075
rcParams['figure.subplot.left'] = 0.125
rcParams['figure.subplot.right'] = 0.9
rcParams['figure.subplot.bottom'] = 0.175
rcParams['figure.subplot.top'] = 0.9
days = mdates.DayLocator() # every day
hours = mdates.HourLocator(interval=3) # every hour
daysFmt = mdates.DateFormatter('%m/%d')
def set_daysFmt(ax):
ax.xaxis.set_major_formatter(daysFmt)
ax.xaxis.set_minor_locator(hours)
$\hspace{10pt}$ accessed from http://tidesandcurrents.noaa.gov/harcon.html?unit=0&timezone=1&id=8723214&name=Virginia+Key&state=FL
In [2]:
# [name, period in hours]
tide_spec = [['M2',12.421],
['S2',12.000],
['N2',12.658],
['K1',23.935],
['M4',6.210],
['O1',25.819],
['M6',4.140],
['MK3',8.177],
['S4',6.000],
['MN4',6.269],
['NU2',12.626],
['S6',4.000],
['MU2',12.872],
['2N2',12.905],
['OO1',22.306],
['LAM2',12.222],
['S1',24.000],
['M1',24.841],
['J1',23.099],
['MM',661.311],
['SSA',4383.076],
['SA',8766.153],
['MSF',354.367],
['MF',327.860],
['RHO',26.723],
['Q1',26.868],
['T2',12.016],
['R2',11.983],
['2Q1',28.006],
['P1',24.0659],
['2SM2',11.607],
['M3',8.280],
['L2',12.192],
['2MK3',8.386],
['K2',11.967],
['M8',3.105],
['MS4',6.103]]
nfreq = 8 #37
species = tide_spec[:nfreq]
In [3]:
def harmonic_analysis(station,pharmonics,time_sec,stage,print_results=False):
npfreq = len(pharmonics)
npsamples = stage.shape[0]
nf2 = 2 * npfreq
nf21, nf22 = nf2+1, nf2+2
aa = zeros((nf22, nf22), float)
ci = zeros(npfreq, float)
si = zeros(npfreq, float)
x = zeros(nf21, float)
trig = zeros(nf22, float)
amplitude = zeros(npfreq, float)
theta = zeros(npfreq, float)
trig[0] = 1.
for i in xrange(npsamples):
for j in xrange(1,nf2,2):
jh = (j-1) / 2
period = pharmonics[jh][1] * 60. * 60.
if period==0.0: print period
fr = 1. / period
omega = 2. * pi * fr
arc = omega * t[i]
j1 = j + 1
trig[j] = cos(arc)
trig[j1] = sin(arc)
trig[nf22-1] = stage[i]
# fill matrix
for i1 in xrange(nf21):
for j1 in xrange(i1,nf22):
aa[i1,j1] = aa[i1,j1] + trig[i1] * trig[j1]
aa[j1,i1] = aa[i1,j1]
# scale diagonal and off diagonals
for i in xrange(nf21):
diagonal = aa[i,i]
aa[i,i] = 1.0
for j in xrange(nf22):
if (j != i):
aa[i,j] = aa[i,j] / diagonal
# solve equations
determinant = splin.det(aa[:nf21,:nf21])
if abs(determinant) < 1.e-20:
print 'singular or ill-conditioned matrix'
print 'determinant = {0}'.format(determinant)
ainv = splin.inv(aa[:nf21,:nf21])
# calculate x
for i in xrange(nf21):
x[i] = 0.
for j in xrange(nf21):
x[i] += ainv[i,j] * aa[j,nf22-1]
# calculate amplitude and phase
harmonic_results = []
for i in xrange(npfreq):
i2 = (i * 2) + 1
ci[i] = x[i2]
si[i] = x[i2+1]
amplitude[i] = sqrt(ci[i] * ci[i] + si[i] * si[i])
ttheta = arctan2(si[i], ci[i]) #* 180. / pi #degrees
if ttheta < 0.:
ttheta += 2. * pi
theta[i] = rad2deg(ttheta)
harmonic_results.append([pharmonics[i][0], pharmonics[i][1], amplitude[i], theta[i]])
#--write the data
if print_results:
print_hr(station,harmonic_results)
return amplitude, theta
def print_hr(station,pharmonics,difference=False):
print station
if difference:
print ' {0:5s} {1:^16s} {2:^16s}'.format('', 'Tidal Effic.', 'Phase Diff.')
else:
print ' {0:5s} {1:^16s} {2:^16s}'.format('', 'Amplitude', 'Phase')
print ' {0:5s} {1:^16s} {2:^16s}'.format('Name', 'meters', 'degrees')
print ' {0:5s} {1:^16s} {2:^16s}'.format('-----', '----------------', '----------------')
for idx,[c,period,pa,ptheta] in enumerate(pharmonics):
print ' {0:5s} {1: 16.3f} {2: 16.1f}'.format(c, pa, ptheta)
print ''
In [4]:
#function to parse data and normalize it
def parse_data(fname):
data_dtype = dtype([('datetime','a20'),('head','f4')])
raw = loadtxt(fname, delimiter=',', dtype=data_dtype, comments='#')
dtime = []
stage = []
for [d,v] in raw:
dt = datetime.datetime.strptime(d, '%m/%d/%Y %H:%M:%S')
dtime.append(dt)
stage.append(v)
dtime = array(dtime)
stage = array(stage)
stagen = np.copy(stage)
stagen -= average(stage)
stagen = signl.detrend(stage)
return dtime, stagen, stage
In [5]:
data_dir = os.path.join('data')
tide_file = 'DE_Tide.csv'
cat_file = 'DE_Catalog.csv'
tide_datetime, tide_v, tide_vr = parse_data(os.path.join(data_dir,tide_file))
tide = ['tide', 0.0, tide_datetime, tide_v, tide_vr]
time_min = tide_datetime[0]
time_max = tide_datetime[-1]
dt_sec = (tide_datetime[1] - tide_datetime[0]).total_seconds()
fo = open(os.path.join(data_dir,cat_file), 'r')
gw_data = []
for line in fo:
if len(line) < 1:
break
if line[0] == '#':
continue
tc = line.strip().replace(' ','').split(',')
station = tc[0]
x = float(tc[2])
dtime, dv, dvr = parse_data(os.path.join(data_dir,tc[1]))
gw_data.append([station, x, dtime, dv, dvr])
time_min = min(time_min, dtime[0])
time_max = max(time_max, dtime[-1])
fo.close()
#--create time_all
tstep = datetime.timedelta(seconds=dt_sec)
dt = time_min
time_all = []
while dt <= time_max:
time_all.append(dt)
dt += tstep
time_all = array(time_all)
#--plot raw normalized data
fig1, ax_all = pyplot.subplots(nrows=1, ncols=1, figsize=(6., 2.5), dpi=300, facecolor='w', edgecolor='k')
ax = ax_all
dtime, dvr = tide[2], tide[4]
ax.plot(dtime, dvr, linewidth=1.0, color='black', label='Tide')
nstat = len(gw_data)
for istat,[station,x,gtime,gv,gvr] in enumerate(gw_data):
fcolor = float(istat) / float(nstat)
ax.plot(gtime, gvr, linewidth=1.0, color=cm.jet_r(fcolor), label=station)
ax.plot([gtime[0], gtime[-1]], [0., 0.], linewidth=0.5, color='0.5', label='_None')
set_daysFmt(ax)
l = ax.legend(loc='best',ncol=7)
l.draw_frame(False)
ax.set_ylim(-0.5,1.5)
fig1.savefig(os.path.join('Figures','DE_RawData.png'), dpi=300)
show()
In [6]:
#def pair_data(pdatetime_in,pv_in,pdatetime_base,pv_base):
# a = pd.DataFrame(pv_in,index=pdatetime_in)
# b = pd.DataFrame(pv_base,index=pdatetime_base)
# c = pd.concat([a,b], axis=1, join='inner')
# pdatetime_out = c.index.to_pydatetime()
# return pdatetime_out, c.values[:,0], c.values[:,1]
def pair_data(pdatetime_in,pv_in,pdatetime_base,pv_base):
dt_min = max(pdatetime_in[0], pdatetime_base[0])
dt_max = min(pdatetime_in[-1], pdatetime_base[-1])
dt_sec = datetime.timedelta(seconds=6. * 60.)
pdatetime_index = []
dt_on = dt_min
while dt_on <= dt_max:
pdatetime_index.append(dt_on)
dt_on += dt_sec
pdatetime_index = array(pdatetime_index)
a = pd.Series(pv_in,index=pdatetime_in)
b = pd.Series(pv_base,index=pdatetime_base)
c = a.reindex(set(a.index).union(pdatetime_index)).sort_index().interpolate('time').ix[pdatetime_index]
d = b.reindex(set(b.index).union(pdatetime_index)).sort_index().interpolate('time').ix[pdatetime_index]
return pdatetime_index, c.values, d.values
def calculate_t(datetime_in, datetime_start=None):
total_sec = (datetime_in[-1] - datetime_in[0]).total_seconds()
nsamples = datetime_in.shape[0]
t_out = linspace(0.,float(total_sec),nsamples)
if datetime_start is not None:
adj_total_sec = (datetime_in[0] - datetime_start).total_seconds()
t_out += adj_total_sec
return t_out
In [7]:
def create_signal(tconstituents,t_cs):
nsamples = t_cs.shape
v_out = zeros(nsamples, float)
for idx,[c,period,a,theta,te,phase] in enumerate(tconstituents):
period *= 60. * 60.
fr = 1. / period
omega = 2. * pi * fr
theta *= pi / 180.
phase *= pi / 180.
v_out += te * a * cos(omega*t_cs - theta - phase)
return v_out
def signal_error(tconstituents, datetime_se, obsv, datetime_start=None):
t_se = calculate_t(datetime_se, datetime_start=datetime_start)
simv = create_signal(tconstituents,t_se)
resid = subtract(simv,obsv)
me = average(resid)
rmse = std(resid)
return me, rmse
In [8]:
#--process tide
station = 'tide'
t = calculate_t(tide_datetime, datetime_start=time_min)
tide_alpha, tide_theta = harmonic_analysis(station,species,t,tide_v,print_results=True)
#--plot simulated and observed tide
pt = calculate_t(time_all)
harmonic_results = []
for idx,[c,period] in enumerate(species):
harmonic_results.append([c, period, tide_alpha[idx], tide_theta[idx], 1.0, 0.0])
#dsec = (time_min-dtime[0]).total_seconds()
#pt += dsec
vr = create_signal(harmonic_results,pt)
vr_paired = create_signal(harmonic_results,t)
me, rmse = signal_error(harmonic_results, tide_datetime, tide_v, datetime_start=time_min)
ctxt = '{0} - ME: {1:7.2g} meters, RMSE: {2:7.2g} meters'.format('Tide', me, rmse)
ax = subplot(1,1,1)
ax.plot(time_all, vr, linewidth=3.0, color='black', zorder=9)
ax.plot(tide_datetime, tide_v, linewidth=0.75, color='red', zorder=10)
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=8, transform=ax.transAxes)
set_daysFmt(ax)
xlabel('Date')
ylabel('Normalized elevation, m')
ax.set_ylim(-0.5,0.5)
show()
ax = subplot(1,1,1)
ax.plot(tide_v, (vr_paired-tide_v), linewidth=0.0, marker='o', markersize=2, color='black')
xlabel('Observed normalized elevation, m')
ylabel('Residual, m')
show()
In [9]:
nprow = 4
npcol = 2
heightp = float(nprow) * 1.25
#--process each groundwater file
fighc, ax_all = pyplot.subplots(nrows=nprow, ncols=npcol, figsize=(10., heightp), dpi=300, facecolor='w', edgecolor='k')
fighc.tight_layout()
pt = calculate_t(time_all)
vr0 = None
# calculate simulated tidal value
for irow in xrange(nprow):
for jcol in xrange(npcol):
ispecies = irow * npcol + jcol
harmonic_results = []
for jdx in xrange(ispecies+1):
c, period= species[jdx]
harmonic_results.append([c, period, tide_alpha[jdx], tide_theta[jdx], 1.0, 0.0])
vr = create_signal(harmonic_results,pt)
# plot the results
ctxt = '{0:2d} tidal constituent'.format(ispecies+1)
if ispecies > 0:
ctxt = '{0}s'.format(ctxt)
ctxt = '{0} ['.format(ctxt)
for [tc, tp, ta, tt, t1, t2] in harmonic_results:
ctxt = '{0} {1},'.format(ctxt, tc)
ctxt = '{0} ]'.format(ctxt)
ax = ax_all[irow, jcol]
if vr0 is not None:
ax.plot(time_all, vr0, linewidth=1.0, color='red')
ax.plot(time_all, vr, linewidth=1.0, color='blue')
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=6, transform=ax.transAxes)
set_daysFmt(ax)
ax.set_ylim(-0.5,0.5)
set_daysFmt(ax)
# save results
vr0 = np.copy(vr)
#if irow+1 < nprow:
# ax.get_xaxis().set_ticks([])
fighc.savefig(os.path.join('Figures','DE_HarmonicConstituents.png'), dpi=300)
show()
In [10]:
station_te = []
station_dtheta = []
for istat,[station, x, gd, gv, gvr] in enumerate(gw_data):
pdatetime, pval, pbase = pair_data(gd, gv, tide_datetime, tide_v)
t = calculate_t(pdatetime, datetime_start=time_min)
alpha, theta = harmonic_analysis('tide',species,t,pbase,print_results=False)
a2, t2 = harmonic_analysis(station,species,t,pval,print_results=False)
te = []
dtheta = []
# write the station data
harmonic_results = []
for idx,[c,period] in enumerate(species):
v1 = a2[idx]/alpha[idx]
v2 = t2[idx]
v3 = theta[idx]
if v2 < 0.0: v2 += 360.
if v2 > 360.: v2 -= 360.
if v3 < 0.0: v3 += 360.
if v3 > 360.: v3 -= 360.
v4 = v2 - v3
te.append(v1)
dtheta.append(v4)
harmonic_results.append([c, period, v1, v4])
station_te.append(te)
station_dtheta.append(dtheta)
print_hr(station,harmonic_results,difference=True)
In [11]:
npcol = 2
nprow = 8
# calculate simulated tidal value
t = calculate_t(time_all)
harmonic_results = []
for idx,[c,period] in enumerate(species):
harmonic_results.append([c, period, tide_alpha[idx], tide_theta[idx], 1.0, 0.0])
vr = create_signal(harmonic_results,t)
me, rmse = signal_error(harmonic_results, tide_datetime, tide_v, datetime_start=time_min)
ctxt = '{0} - ME: {1:7.2g} meters, RMSE: {2:7.2g} meters'.format('Tide', me, rmse)
ifig = 1
fig1, ax_all = pyplot.subplots(nrows=1, ncols=1, figsize=(6., 1.0), dpi=300, facecolor='w', edgecolor='k')
ax = ax_all
ax.plot(time_all, vr, linewidth=1.0)
ax.plot(tide_datetime, tide_v, linewidth=0.0, marker='x', markersize=2)
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=6, transform=ax.transAxes)
set_daysFmt(ax)
ax.set_ylim(-0.5,0.5)
fig1.savefig(os.path.join('Figures','DE_TidalFit.png'), dpi=300)
#--process each groundwater file
ifig = 0
fig2, ax_all = pyplot.subplots(nrows=nprow, ncols=npcol, figsize=(6., 6.5), dpi=300, facecolor='w', edgecolor='k')
print ax_all.shape
#fig2 = figure(num=None, figsize=(6., 6.5), dpi=300, facecolor='w', edgecolor='k')
fig2.tight_layout()
irow = 0
icol = 0
for istat,[station, x, gd, gv, gvr] in enumerate(gw_data):
ifig += 1
harmonic_results = []
for idx,[c,period] in enumerate(species):
harmonic_results.append([c, period, tide_alpha[idx], tide_theta[idx], station_te[istat][idx], station_dtheta[istat][idx]])
vr = create_signal(harmonic_results, t)
me, rmse = signal_error(harmonic_results, gd, gv, datetime_start=time_min)
ctxt = '{0} - ME: {1:7.2g} meters RMSE: {2:7.2g} meters'.format(station, me, rmse)
#ax = fig2.add_subplot(nprow,npcol,ifig)
ax = ax_all[irow,icol]
icol += 1
if icol > 1:
irow += 1
icol = 0
#print ax
ax.plot(time_all, vr, linewidth=1.0)
ax.plot(gd, gv, linewidth=0.0, marker='x', markersize=2)
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=6, transform=ax.transAxes)
set_daysFmt(ax)
ax.set_ylim(-0.5,0.5)
#for idx in xrange(1,ifig-1):
#ax = subplot(nprow,npcol,idx)
for idx in xrange(0,nprow-1):
for jdx in xrange(npcol):
ax = ax_all[idx, jdx]
ax.get_xaxis().set_ticks([])
fig2.savefig(os.path.join('Figures','DE_GWObsFit.png'), dpi=300)
show()
In [12]:
npcol = 2
nprow = 1 #+ (len(gw_data)) / npcol
nstat = len(gw_data)
#figure(num=None, figsize=(12, 12), dpi=300, facecolor='w', edgecolor='k')
fig, ax_all = pyplot.subplots(nrows=2, ncols=1, figsize=(3.75, 7.5), dpi=300, facecolor='w', edgecolor='k')
fig.tight_layout()
figte, axte = pyplot.subplots(nrows=1, ncols=1, figsize=(4.25, 4.25), dpi=300, facecolor='w', edgecolor='k')
figte.tight_layout()
for istat,[station, xpos, gd, gv, gvr] in enumerate(gw_data):
if xpos < 0: continue
fcolor = float(istat) / float(nstat)
te = array(station_te[istat])
dt = array(station_dtheta[istat])
Dalpha = []
Dtheta = []
for idx,[c,period] in enumerate(species):
#period *= 60. * 60. # period in seconds
period /= 24. # period in days
fr = 1. / period
omega = 2. * pi / fr
v = omega / (2.0 * power((log(te[idx])/xpos),2.))
Dalpha.append(log10(v))
dtrad = dt[idx] * pi / 180.
v = 1. / (2. * omega * power((dtrad/xpos), 2.0))
Dtheta.append(log10(v))
x = average(te)
xerr = 2.*std(te)
Dalpha = array(Dalpha)
y = average(Dalpha)
yerr = 2.*std(Dalpha)
#ax = subplot(2,2,1)
ax = ax_all[0]
ax.errorbar(x, y, yerr=yerr, xerr=xerr, fmt='o', color=cm.jet_r(fcolor), label=station)
#
x = average(Dalpha)
xerr = 2.*std(Dalpha)
Dtheta = array(Dtheta)
y = average(Dtheta)
yerr = 2.*std(Dtheta)
#ax = subplot(2,2,2)
ax = ax_all[1]
ax.errorbar(x, y, yerr=yerr, xerr=xerr, fmt='o', color=cm.jet_r(fcolor), label=station)
# te figure
r = np.ones((len(te)), np.float) * xpos
#print r
if '580A' not in station:
axte.semilogx(r, te, marker='o', linewidth=0., color=cm.jet_r(fcolor), label=station)
#rval = xpos - 0.2
#if istat > 0:
# rval = xpos + 2000
#if istat == 0 or istat+1==len(gw_data):
# for idx,[c,period] in enumerate(species):
# axte.text(rval, te[idx], c, ha='center', va='center', size=6)
D_min, D_max = -2, 10
#ax = subplot(2,2,1)
ax = ax_all[0]
l = ax.legend(loc='best',ncol=2,numpoints=1)
l.draw_frame(False)
ax.set_ylim(D_min, D_max)
ax.set_xlabel('Tidal efficiency, unitless')
ax.set_ylabel(r'log$_{10}$(Hydraulic diffusivity), $m^{2}/d$')
ax.xaxis.set_label_coords(0.5, -0.05)
ax.yaxis.set_label_coords(-0.05, 0.5)
#ax = subplot(2,2,2)
ax = ax_all[1]
ax.plot([D_min, D_max],[D_min, D_max], color='black', linewidth=1)
ax.set_xlim(D_min, D_max)
ax.set_ylim(D_min, D_max)
ax.set_xlabel(r'log$_{10}$(Hydraulic diffusivity (TE)), $m^{2}/d$')
ax.set_ylabel(r'log$_{10}$(Hydraulic diffusivity ($\Delta \theta$)), $m^{2}/d$')
ax.xaxis.set_label_coords(0.5, -0.05)
ax.yaxis.set_label_coords(-0.05, 0.5)
l = axte.legend(loc='best',ncol=2,numpoints=1)
l.draw_frame(False)
axte.set_ylim(0, 1)
axte.set_xlabel('Distance from Biscayne Bay, m')
axte.set_ylabel('Tidal efficiency, unitless')
axte.xaxis.set_label_coords(0.5, -0.045)
axte.yaxis.set_label_coords(-0.05, 0.5)
fig.savefig(os.path.join('Figures','DE_CalcDh.png'), dpi=300)
figte.savefig(os.path.join('Figures','DE_TidalEfficiency.png'), dpi=300)
show()
In [12]:
In [12]: