For results including the Northern Boundary, plot the Foreman et al (1995) model results against our model results, along a thalweg (the thalweg is defined and plotted by Nancy in computer_thalweg.ipynb)
# imports
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import netCDF4 as NC
#get thalwag co-ordinates (provided by Nancy)
lines = np.loadtxt('../../tools/analysis_tools/thalweg.txt', delimiter=" ", unpack=False)
lines = lines.astype(int)
#get Foreman data
#data in /ocean/klesouef/meopar/tools/compare_tides
path = '/ocean/klesouef/meopar/tools/compare_tides/'
consts = np.loadtxt(path+'foreman_m2_thalweg.txt')
m2_amp_foreman = consts[:,0]
m2_pha_foreman = consts[:,1]
m2_pha_foreman[m2_pha_foreman>180]=m2_pha_foreman[m2_pha_foreman>180]-360
consts = np.loadtxt(path+'foreman_k1_thalweg.txt')
k1_amp_foreman = consts[:,0]
k1_pha_foreman = consts[:,1]
k1_pha_foreman[k1_pha_foreman>180]=k1_pha_foreman[k1_pha_foreman>180]-360
#get model data
from salishsea_tools import tidetools
runname = '20dec25dec','26dec31dec', '1jan5jan', '6jan10jan', '11jan20jan', '21jan30jan', \
'31jan9feb', '10feb19feb','20feb1mar','2mar11mar','12mar21mar','22mar31mar', '1apr10apr', \
'11apr20apr', '21apr30apr', '1may10may', '11may20may', '21may30may', '31may9jun', \
'10jun19jun', '20jun29jun', '30jun9jul'
resultsloc = '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/'
#resultsloc = '/ocean/sallen/allen/research/Meopar/myResults/SalishSea/'
#runname = 'M2only'
mod_M2_amp, mod_K1_amp, mod_M2_pha, mod_K1_pha = tidetools.get_amp_phase_data(runname,resultsloc)
#rename our model results
m2_amp_model_thalweg = mod_M2_amp[lines[:,0],lines[:,1]]*100
m2_pha_model_thalweg = mod_M2_pha[lines[:,0],lines[:,1]]
k1_amp_model_thalweg = mod_K1_amp[lines[:,0],lines[:,1]]*100
k1_pha_model_thalweg = mod_K1_pha[lines[:,0],lines[:,1]]
#plot everything
thalnum = range(0,len(m2_amp_model_thalweg))
#M2
plt.figure(figsize=(12,5))
#the Foreman model doesn't quite cover our domain, so ignore some points
plt.plot(thalnum[60:620],m2_amp_foreman[60:620],'b',label='Foreman')
plt.plot(thalnum,m2_amp_model_thalweg,'m',label='model')
plt.legend(loc=5)
plt.xlabel('thalweg point')
plt.ylabel('m2 amp [cm]')
plt.ylim([30,100])
plt.xlim([0,max(thalnum)])
plt.title('Foreman model and our model: M2 amplitude along thalweg')
plt.figure(figsize=(12,5))
plt.plot(thalnum[60:620],m2_pha_foreman[60:620],'b',label='Foreman')
plt.plot(thalnum,m2_pha_model_thalweg,'m',label='model')
plt.ylim([-120,50])
plt.legend(loc=4)
plt.xlim([0,max(thalnum)])
plt.xlabel('thalweg point')
plt.ylabel('m2 pha [deg]')
plt.title('Foreman model and our model: M2 phase along thalweg')
#K1
plt.figure(figsize=(12,5))
#the Foreman model doesn't quite cover our domain, so ignore some points
plt.plot(thalnum[60:620],k1_amp_foreman[60:620],'b',label='Foreman')
plt.plot(thalnum,k1_amp_model_thalweg,'m',label='model')
plt.legend(loc=4)
plt.xlabel('thalweg point')
plt.ylabel('k1 amp [cm]')
plt.ylim([40,90])
plt.xlim([0,max(thalnum)])
plt.title('Foreman model and our model: K1 amplitude along thalweg')
plt.figure(figsize=(12,5))
plt.plot(thalnum[60:620],k1_pha_foreman[60:620],'b',label='Foreman')
plt.plot(thalnum,k1_pha_model_thalweg,'m',label='model')
plt.ylim([-120,-50])
plt.legend(loc=4)
plt.xlim([0,max(thalnum)])
plt.xlabel('thalweg point')
plt.ylabel('k1 pha [deg]')
plt.title('Foreman model and our model: K1 phase along thalweg')
<matplotlib.text.Text at 0x5f87c50>
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.plot(thalnum[60:620],m2_amp_model_thalweg[60:620]-m2_amp_foreman[60:620])
plt.title('M2 amp (cm)')
plt.subplot(2,2,2)
plt.plot(thalnum[60:620],m2_pha_model_thalweg[60:620]-m2_pha_foreman[60:620])
plt.title('M2 pha (deg)')
plt.ylim([-10,15])
plt.subplot(2,2,3)
plt.plot(thalnum[60:620],k1_amp_model_thalweg[60:620]-k1_amp_foreman[60:620])
plt.title('K1 amp (cm)')
plt.subplot(2,2,4)
plt.plot(thalnum[60:620],k1_pha_model_thalweg[60:620]-k1_pha_foreman[60:620])
plt.title('K1 pha (deg)')
<matplotlib.text.Text at 0x5de1410>
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.plot(m2_amp_model_thalweg[60:620],m2_amp_foreman[60:620],'o',np.arange(40,100),np.arange(40,100),'-')
plt.title('M2 amp (cm)')
plt.subplot(2,2,2)
plt.plot(m2_pha_model_thalweg[60:620],m2_pha_foreman[60:620],'o',np.arange(-120,60),np.arange(-120,60),'-')
plt.title('M2 pha (deg)')
plt.subplot(2,2,3)
plt.plot(k1_amp_model_thalweg[60:620],k1_amp_foreman[60:620],'o',np.arange(45,90),np.arange(45,90),'-')
plt.title('K1 amp (cm)')
plt.subplot(2,2,4)
plt.plot(k1_pha_model_thalweg[60:620],k1_pha_foreman[60:620],'o',np.arange(-110,-55),np.arange(-110,-55),'-')
plt.title('K1 pha (deg)')
<matplotlib.text.Text at 0x40c9c90>