In [42]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#import seaborn as sns
#sns.set(style="darkgrid")
#sns.set(style="whitegrid")
from pyspec import spectrum
In [43]:
plt.rcParams.update({'font.size': 12
, 'legend.markerscale': 1., 'axes.titlesize': 12, 'axes.labelsize' : 12,
'legend.fontsize' : 10,'legend.handlelength': 3})
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
In [44]:
color2 = '#6495ed'
color1 = '#ff6347'
color5 = '#8470ff'
color3 = '#3cb371'
color4 = '#ffd700'
color6 = '#ba55d3'
lw=3
aph=.7
In [45]:
data_path = './outputs/'
slab1=np.load(data_path+'spectra_'+'slab1_bcf.npz')
slab2=np.load(data_path+'spectra_'+'slab2_bcf.npz')
slab3=np.load(data_path+'spectra_'+'slab3_bcf.npz')
sampled = np.load(data_path+"spec_0m_sampling.npz")
In [46]:
## -2 and -3 slopes in the loglog space
ks = np.array([1.e-3,1])
Es2 = .2e-4*(ks**(-2))
Es3 = .35e-6*(ks**(-3))
rd1 = 22.64 # [km]
Enoise = np.ones(2)*2.*1.e-4
def add_second_axis(ax1):
""" Add a x-axis at the top of the spectra figures """
ax2 = ax1.twiny()
ax2.set_xscale('log')
ax2.set_xlim(ax1.axis()[0], ax1.axis()[1])
kp = 1./np.array([500.,200.,100.,40.,20.,10.,5.])
lp=np.array([500,200,100,40,20,10,5])
ax2.set_xticks(kp)
ax2.set_xticklabels(lp)
plt.xlabel('Wavelength [km]')
def plt_spec_error(x=0.034,y=0.15, sn=200):
Erl, Eru = spectrum.spec_error(y,sn=sn)
err=np.array([y-Erl,Eru-y])[...,np.newaxis]
plt.errorbar(x, y, yerr=err, color='k',fmt='.')
plt.text(x+0.005,y-0.025,r'95 $\%$',fontsize=10)
def plt_adcp_spectrum(slab,vlevel=1,lw=3):
""" Plots ADCP spectrum in the given vertical level
slab is a dictionary contaning the spectra """
if vlevel==1:
ltit = r'26-50 m, 400 DOF'
fig_num = 'a'
elif vlevel==2:
ltit=r'58-98 m, 400 DOF'
fig_num = 'b'
elif vlevel==3:
ltit=r'106-202 m, 400 DOF'
fig_num = 'c'
fig = plt.figure(facecolor='w', figsize=(12.,10.))
ax1 = fig.add_subplot(111)
#ax1.fill_between(slab['k'],slab['Eul']/2.,slab['Euu']/2., color=color1, alpha=0.35)
#ax1.fill_between(slab['k'],slab['Evl']/2.,slab['Evu']/2., color=color2, alpha=0.35)
#ax1.fill_between(slab['k'],slab['Eufl']/2.,slab['Eufu']/2., color=color1, alpha=0.35)
#ax1.fill_between(slab['k'],slab['Evfl']/2.,slab['Evfu']/2., color=color2, alpha=0.35)
#ax1.set_xscale('log'); ax1.set_yscale('log')
ax1.loglog(slab['k'],slab['Eu']/2.,color=color1,linewidth=lw,
label=r'$\hat{C}^u$: across-track')
ax1.loglog(slab['k'],slab['Ev']/2.,color=color2,linewidth=lw,
label=r'$\hat{C}^v$: along-track')
ax1.loglog(slab['k'],slab['Euf']/2.,'--',color=color1,linewidth=lw)
ax1.loglog(slab['k'],slab['Evf']/2.,'--',color=color2,linewidth=lw)
ax1.loglog(slab['kK'],slab['Kpsi']/2,color=color3,linewidth=lw,
label='$\hat{C}^\psi$: rotational')
ax1.loglog(slab['kK'],slab['Kphi']/2,color=color4,linewidth=lw,
label='$\hat{C}^\phi$: divergent')
plt_spec_error()
ax1.loglog(slab['kK'],slab['Kpsif']/2,'--',color=color3,linewidth=lw)
ax1.loglog(slab['kK'],slab['Kphif']/2,'--',color=color4,linewidth=lw)
ax1.axis((1./(1000),1./4,.4e-5,10))
ax1.loglog(ks,Es2,'-', color='0.5',linewidth=2.)
ax1.loglog(ks,Es3,'-', color='0.5',linewidth=2.)
ax1.axis((1./(1000),1./4,.4e-5,10))
plt.text(0.002, 1.941,u'-2')
plt.text(0.006, 1.951,u'-3')
plt.xlabel('Along-track wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [ m$^{2}$ s$^{-2}$/ cpkm]')
lg = plt.legend(loc=3,title=ltit, numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,5e0))
plt.text(1./20, 3., "llc 4320 simulation", size=25, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
plt.text(1./5, 3., fig_num, size=35, rotation=0.)
add_second_axis(ax1)
plt.savefig('figs/spec_model_slab'+str(vlevel),bbox_inches='tight')
In [ ]:
In [47]:
## 58-98 m
plt_adcp_spectrum(slab2,vlevel=2)
In [48]:
## 106-202 m
plt_adcp_spectrum(slab3,vlevel=3)
In [49]:
aph=.35
Eerr = 1.e-2
std_E = (1/np.sqrt(400))
Eerrl = Eerr/(1 + 2*std_E)
Eerru = Eerr/(1 - 2*std_E)
win = np.hanning(sampled['k'].size)
fac = (sampled['k'].size/(win**2).sum())
In [50]:
lw=1
fig = plt.figure(figsize=(8.27/2-.25,11.69/3-.25))
ax1 = fig.add_subplot(111)
ax1.loglog(sampled['k'],fac*sampled['Eu']/2.,color=color1,linewidth=lw,
label=r'snapshots')
ax1.loglog(sampled['k'],fac*sampled['Eus']/2.,color=color3,linewidth=lw,
label=r'sampled')
ax1.loglog(sampled['ksb'],fac*sampled['Eusb']/2.,color='b',linewidth=lw,
label=r'sampled+block-averaged',alpha=aph)
# error bar reference
#ax1.loglog([.075,.075],[Eerrl,Eerru],'k',linewidth=2.)
#ax1.text (.0785,0.0088,r'$95\%$',fontsize=20)
plt_spec_error()
ax1.loglog(ks,fac*Es2,'-', color='0.5',linewidth=1.)
ax1.loglog(ks,fac*Es3,'-', color='0.5',linewidth=1.)
plt.text(0.00185, 4.25,u'-2')
plt.text(0.00375, 4.25,u'-3')
plt.xlabel('Along-track wavenumber [cpkm]')
plt.ylabel(u'KE density [m$^{2}$ s$^{-2}$/ cpkm]')
lg = plt.legend(loc=3,title='across-track spectrum', numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1.,1./1.e6,1.e1))
lg = plt.legend(loc=3,title="", numpoints=1,ncol=1)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,10e0))
plt.text(1./20, 3., "llc 4320 simulation", size=10, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
add_second_axis(ax1)
plt.savefig('figs/spec_model_0m_sampled',bbox_inches='tight')
plt.savefig('figs/spec_model_0m_sampled.eps',bbox_inches='tight')
plt.savefig('figs/spec_model_0m_sampled.pdf',bbox_inches='tight')
In [51]:
#compare with adcp
In [52]:
slab2_adcp=np.load('../adcp/outputs/adcp_spec_slab2.npz')
slab2_model=np.load('1d_spec_100m.npz')
In [57]:
eiso = np.load('Eiso_uv_small.npz')
eiso.keys()
Out[57]:
In [60]:
fig = plt.figure(figsize=(8.27/2-.25,11.69/3-.25))
ax1 = fig.add_subplot(111)
ax1.fill_between(slab2['k'],fac*slab2['Eul']/2.,fac*slab2['Euu']/2., color=color1, alpha=0.35)
ax1.fill_between(slab2['k'],fac*slab2['Evl']/2.,fac*slab2['Evu']/2., color=color2, alpha=0.35)
ax1.fill_between(slab2_adcp['k'],fac*slab2_adcp['Eul']/2,fac*slab2_adcp['Euu']/2, color=color1, alpha=0.35)
ax1.fill_between(slab2_adcp['k'],fac*slab2_adcp['Evl']/2,fac*slab2_adcp['Evu']/2, color=color2, alpha=0.35)
ax1.set_xscale('log'); ax1.set_yscale('log')
ax1.loglog(slab2_adcp['k'],fac*slab2_adcp['Eu']/2.,color=color1,linewidth=lw,
label=r'$\hat{C}^u$: across-track, adcp')
ax1.loglog(slab2_adcp['k'],fac*slab2_adcp['Ev']/2.,color=color2,linewidth=lw,
label=r'$\hat{C}^v$: along-track, adcp')
ax1.loglog(slab2['k'],fac*slab2['Eu']/2.,'--',color=color1,linewidth=lw,
label=r'$\hat{C}^u$: across-track, model')
ax1.loglog(slab2['k'],fac*slab2['Ev']/2.,'--',color=color2,linewidth=lw,
label=r'$\hat{C}^v$: along-track, model')
# error bar reference
#ax1.loglog([.075,.075],[Eerrl,Eerru],'k',linewidth=2.)
#ax1.text (.0785,0.0088,r'$95\%$',fontsize=20)
#plt_spec_error()
ax1.loglog(ks,fac*Es2,'-', color='0.5',linewidth=1.)
ax1.loglog(ks,fac*Es3,'-', color='0.5',linewidth=1.)
plt.text(0.00185, 4.25,u'-2')
plt.text(0.00375, 4.25,u'-3')
plt.xlabel('Along-track wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [m$^{2}$ s$^{-2}$/ cpkm]')
lg = plt.legend(loc=3,title='', numpoints=1,ncol=1)
lg.draw_frame(False)
plt.text(1./20, 3., "llc 4320 simulation \n vs. ADCP", size=10, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
plt.axis((1./1.e3,1.,1./1.e6,1.e1))
plt.axis((1./1.e3,1./4.,.5/1.e4,1e1))
add_second_axis(ax1)
plt.savefig('figs/spec_model_slab2_comparison',bbox_inches='tight')
plt.savefig('figs/spec_model_slab2_comparison.eps',bbox_inches='tight')
plt.savefig('figs/spec_model_slab2_comparison.pdf',bbox_inches='tight')
In [54]:
fig = plt.figure(facecolor='w', figsize=(12.,10.))
ax1 = fig.add_subplot(111)
#ax1.fill_between(slab1['k'],slab1['Eufl']/2.,slab1['Eufu']/2., color=color1, alpha=0.35)
#ax1.fill_between(slab1['k'],slab1['Evfl']/2.,slab1['Evfu']/2., color=color2, alpha=0.35)
#ax1.fill_between(slab2['k'],slab2['Eufl']/2.,slab2['Eufu']/2., color=color1, alpha=0.35)
#ax1.fill_between(slab2['k'],slab2['Evfl']/2.,slab2['Evfu']/2., color=color2, alpha=0.35)
#ax1.fill_between(slab3['k'],slab3['Eufl']/2.,slab3['Eufu']/2., color=color1, alpha=0.35)
#ax1.fill_between(slab3['k'],slab3['Evfl']/2.,slab3['Evfu']/2., color=color2, alpha=0.35)
ax1.loglog(slab1['k'],slab1['Euf']/2,'-',color=color1,alpha=.8,linewidth=lw,label='26-50 m')
ax1.loglog(slab2['k'],slab2['Euf']/2,'-',color=color3,alpha=.8,linewidth=lw,label='58-98 m')
ax1.loglog(slab3['k'],slab3['Euf']/2,'-',color=color2,alpha=.8,linewidth=lw,label='106-202 m')
ax1.loglog(slab1['k'],slab1['Evf']/2,'--',color=color1,alpha=.8,linewidth=lw)
ax1.loglog(slab2['k'],slab2['Evf']/2,'--',color=color3,alpha=.8,linewidth=lw)
ax1.loglog(slab3['k'],slab3['Evf']/2,'--',color=color2,alpha=.8,linewidth=lw)
plt_spec_error()
ax1.loglog(ks,Es2/3.8,'-', color='0.5',linewidth=2.)
ax1.loglog(ks,Es3/1.5,'-', color='0.5',linewidth=2.)
plt.text(1./20, 3., "llc 4320 simulation", size=25, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
plt.text(0.002, 1.941,u'-2')
plt.text(0.006, 1.951,u'-3')
plt.xlabel('Along-track wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [ m$^{2}$ s$^{-2}$/ cpkm]')
lg = plt.legend(loc=3,title='', numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,5e0))
add_second_axis(ax1)
plt.savefig('figs/spec_model_depth_dependence',bbox_inches='tight')
In [55]:
Eerrl,Eerru, 0.1
Out[55]:
In [56]:
Erl, Eru = spectrum.spec_error(0.1,sn=200)
Erl, Eru,0.1
Out[56]:
In [ ]: