Compute and plot spectra of stochastically-generated data


In [1]:
import numpy as np
import scipy.signal
import scipy as sp
import sys

import matplotlib.pyplot as plt
%matplotlib inline

sys.path.append('../src/')
from aux_func import *
import aux_func_3dfields as my


Vendor:  Continuum Analytics, Inc.
Package: mkl
Message: trial mode expires in 20 days

In [2]:
plt.rcParams.update({'font.size': 25, 'legend.handlelength'  : 1.5
    , 'legend.markerscale': 1.})
plt.rc('xtick', labelsize=25) 
plt.rc('ytick', labelsize=25)

In [3]:
color1 = '#ff6347'
color2 = '#6495ed'
color3 = '#ffd700'
color4 = '#8470ff'
color5 = '#ff69b4'
color6 = '#006400'

lw1=3
aph=.7

In [5]:
# parameters
L = 800.    # [km]
dx = 2.5
dk = 1./L/dx
kNy = 1./(2.*dx)

# create isotropic spectral ramp
k = np.arange(-kNy,kNy,dk)
k = k/k[-1] # non-dimensional form to reduce numerical error
ki,li = np.meshgrid(k,k)
K = np.sqrt((ki**2) + (li**2))

kt = 1.e-10  # low wavenumber cutoff to avoid numerical error 
Esyn2D = (1./K**3)/K
Esyn2D = np.ma.masked_array(Esyn2D,K<=kt)
Esyn2D2 = (1./K**2)/K
Esyn2D2 = np.ma.masked_array(Esyn2D2,K<=kt)

In [6]:
# plot 2D spectral ramp (K^{-3} rolloff)
fig = plt.figure(facecolor='w', figsize=(12.,10.))
plt.contourf(k,k,np.log10(Esyn2D), 25,cmap='Spectral_r')
cb = plt.colorbar()
cb.set_label(u'Spectral density  [(m$^2$s$^{-2}$)/(cycles/km)$^2$]')
plt.axis('equal')
plt.xlim(-.1,.1)
plt.ylim(-.1,.1)

plt.xlabel('Zonal wavenumber [cycles/km]')
plt.ylabel('Meridional wavenumber [cycles/km]')
plt.savefig('figs/2d_Esyn')



In [7]:
# plot 2D spectral ramp (K^{-2} rolloff)
fig = plt.figure(facecolor='w', figsize=(12.,10.))
#plt.contour(k,k,np.log10(Esyn2D2), 25)#,cmap='Spectral_r')
plt.contour(np.log10(Esyn2D2))
#cb = plt.colorbar()
#cb.set_label(u'Spectral density  [(m$^2$s$^{-2}$)/(cycles/km)$^2$]')
plt.axis('equal')
plt.xlim(-.1,.1)
plt.ylim(-.1,.1)

plt.xlabel('Zonal wavenumber [cycles/km]')
plt.ylabel('Meridional wavenumber [cycles/km]')
plt.savefig('figs/2d_Esyn2')



In [ ]:


In [8]:
# create random phase
n = 10
ix,jx=Esyn2D.shape
pha = 2*np.pi*(np.random.rand(ix*jx*n).reshape(ix,jx,n))  # assume phases are correlated
ii = np.complex(0,1)
pha = np.cos(pha) + ii*np.sin(pha)
Esyn2D = np.repeat(Esyn2D,n).reshape(ix,jx,n)
Esyn2D = np.sqrt(Esyn2D)*pha
Esyn2D2 = np.repeat(Esyn2D2,n).reshape(ix,jx,n)
Esyn2D2 = np.sqrt(Esyn2D2)*pha

# back to fourier coefs
an = np.fft.fftshift(Esyn2D)*((dk*dk)*((ix*jx)**2))
an2 = np.fft.fftshift(Esyn2D2)*((dk*dk)*((ix*jx)**2))

# back to physical space (i.e., create synthetic u and v)
U = np.fft.ifft2(an,axes=(0,1))
U = 2.5*(U/U.max())   # normalize to have the same level of energy as in llc4320
u = np.real(U)
v = np.imag(U)

U2 = np.fft.ifft2(an2,axes=(0,1))
U2 = 2.5*(U2/U2.max())
u2 = np.real(U2)
v2 = np.imag(U2)

up,vp,us,vs = ps(u,v,dx,dx)
up2,vp2,us2,vs2 = ps(u2,v2,dx,dx)

# filter divergent part of the flow (cut-off about 40 km)
nx = 40/dx
ny = 40/dx
x, y = np.mgrid[-nx/2:nx/2, -ny/2:ny/2]
    
rx = 40./dx
ry = 40./dx

g = np.exp( -  ( (x/rx)**2 + (y/ry)**2 ) )
g = g/g.sum()

ix,jx,kx = up.shape
upf = np.zeros(up.shape)
vpf = np.zeros(vp.shape)
upf2 = np.zeros(up2.shape)
vpf2 = np.zeros(vp2.shape)

In [9]:
for i in range(kx):
    upm = up[:,:,i].mean()
    upi = up[:,:,i]
    upf[:,:,i] = sp.signal.convolve2d(upi-upm,g, mode='same') + upm 
    vpm = vp[:,:,i].mean()
    vpi = vp[:,:,i]
    vpf[:,:,i] = sp.signal.convolve2d(vpi-vpm,g, mode='same') + vpm 
    del upi, upm

    upm = up2[:,:,i].mean()
    upi = up2[:,:,i]
    upf2[:,:,i] = sp.signal.convolve2d(upi-upm,g, mode='same') + upm 
    vpm = vp2[:,:,i].mean()
    vpi = vp2[:,:,i]
    vpf2[:,:,i] = sp.signal.convolve2d(vpi-vpm,g, mode='same') + vpm 
    del upi, upm

# the divergent part associated with 'small scales'
ud = up-upf
vd = vp-vpf
ud2 = up2-upf2
vd2 = vp2-vpf2

# make the flow slightly div.

nd = 1.  # div. to non-div variance ratio
uu = nd*ud + us
vv = nd*vd + vs

nd = .25*(ud.std()/ud2.std())  # signal to noise ratio in the <40 km subrange
uu2 = nd*ud2 + us
vv2 = nd*vd2 + vs

EUn2,l2d,k2d,_,_,_,_=spec_est2(uu2,dx,dx,win=True)
EVn2,_,_,_,_,_,_=spec_est2(vv2,dx,dx,win=True)
E2D_n2=(EUn2+EVn2)/2.

EUn,l2d,k2d,_,_,_,_=spec_est2(uu,dx,dx,win=True)
EVn,_,_,_,_,_,_=spec_est2(vv,dx,dx,win=True)
E2D_n=(EUn+EVn)/2.

# add white noise and divergent flow at small scales
nn= 0.0  # noise-to-signal ratio
Au = nn*us.std()
Av = nn*vs.std()
ix,jx,kx = us.shape
nu = Au*(np.random.randn(ix*jx*kx)).reshape(ix,jx,kx)
nv = Av*(np.random.randn(ix*jx*kx)).reshape(ix,jx,kx)

# spectral window
ix,jx,kx=u.shape
window = np.repeat(np.hanning(ix),jx).reshape(ix,jx)
window=np.repeat(window,kx).reshape(ix,jx,kx)

# total flow
Eut,kut,dku,kuNy = my.spec_est_meridional(u*window,dx)
Evt,kvt,dkv,kvNy = my.spec_est_meridional(v*window,dx)
Evt=Evt.mean(axis=1)
Eut=Eut.mean(axis=1)

# horizontally non-divergent
Eu,ku,dku,kuNy = my.spec_est_meridional(us*window,dx)
Ev,kv,dkv,kvNy = my.spec_est_meridional(vs*window,dx)
Ev=Ev.mean(axis=1)
Eu=Eu.mean(axis=1)

# add divergent flow at small scales and random noise
Eun,_,_,_ = my.spec_est_meridional((uu+nu)*window,dx)
Evn,_,_,_ = my.spec_est_meridional((vv+nv)*window,dx)
Evn=Evn.mean(axis=1)
Eun=Eun.mean(axis=1)
Eun2,_,_,_ = my.spec_est_meridional((uu2+nu)*window,dx)
Evn2,_,_,_ = my.spec_est_meridional((vv2+nv)*window,dx)
Evn2=Evn2.mean(axis=1)
Eun2=Eun2.mean(axis=1)

Ek = (np.sum(Esyn2D,axis=1)*dk)[k.size/2:]
El = (np.sum(Esyn2D,axis=0)*dk)[k.size/2:]
k = k[k.size/2:]

# add random noise
Enoise=.1*((Ev[240]+Eu[240])/2.)*np.ones(Ev.size)
Eunn,Evnn=Eu+Enoise,Ev+Enoise

# mask very low and very high wavenumbers
L = 1./ku
fm = ((L<=5)|(L>=800))
Eu = np.ma.masked_array(Eu,fm)
Ev = np.ma.masked_array(Ev,fm)
Eun = np.ma.masked_array(Eun,fm)
Evn = np.ma.masked_array(Evn,fm)
Eun2 = np.ma.masked_array(Eun2,fm)
Evn2 = np.ma.masked_array(Evn2,fm)

# compute ratios in an arbitrary range
f = ((L>=40)&(L<=200))
rn = ((Eun/Evn)[f]).mean()
r = ((Eu/Ev)[f]).mean()
rn2 = ((Eun2/Evn2)[f]).mean()

# BCF decomposition of synthetic data
Kpsi_nd,Kphi_nd,kK=bcf(kv,Eu,Ev)
Kpsi_d1,Kphi_d1,_=bcf(kv,Eun,Evn)
Kpsi_d2,Kphi_d2,_=bcf(kv,Eun2,Evn2)


aux_func.py:96: ComplexWarning: Casting complex values to real discards the imaginary part
  EU[i+1:]*np.cosh(s[i]-s[i+1:]))*ds).sum()
aux_func.py:98: ComplexWarning: Casting complex values to real discards the imaginary part
  EU[i+1:]*np.sinh(s[i]-s[i+1:]))*ds).sum()

In [10]:
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]')

In [11]:
# -2 and -3 slopes in the loglog space
ks2 = np.array([1.e-3,1.])
Es2 = .5e-6*(ks2**(-2))

ks3 = np.array([1.e-3,1.])
Es3 = .4e-7*(ks3**(-3))

# nd + d
fig = plt.figure(facecolor='w', figsize=(11.,12.))
plt.loglog(kut,Eut,color=color1,label='across-track',linewidth=4.,alpha=.5)
plt.loglog(kvt,Evt,color=color2,label='along-track',linewidth=4.,alpha=.5)
plt.loglog(ks3,Es3,'--',color='k',linewidth=4.,alpha=.5)
plt.text(0.0041159850623393082, 0.29427271762092821,u'$\kappa^{-3}$')
plt.axis((1./(400),1./2.5,1.e-6,1.))
plt.ylabel('Spectral density  [m$^2$ s$^{-2}$/cpkm]')
plt.xlabel('Along-track wavenumber  [cpkm]')
lg = plt.legend(loc=1,title= u'', prop={'size':22}, numpoints=1)
lg.draw_frame(False)
my.leg_width(lg,5.)
figtit = 'figs/EuEv_synthetic_total.png'
plt.savefig(figtit,format='png', bbox_inches='tight')


/Users/crocha/anaconda/lib/python2.7/site-packages/numpy/core/numeric.py:462: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

In [18]:
# nd + high-pass d
fig = plt.figure(facecolor='w', figsize=(11.,12.))
ax1 = fig.add_subplot(111)
p1,= ax1.loglog(ku,Eu,color=color1,label=u'$\hat{C}^u$: across-track, rotational',linewidth=3.)
ax1.loglog(kv,Ev,color=color2,label=u'$\hat{C}^v$: along-track, rotational',linewidth=3.)
p2, =ax1.loglog(ku,Eun,color=color3,label=u'rot. + high-pass div. (k$^{-3}$), across-track',linewidth=3.)
plt.loglog(kv,Evn,color=color6,label=u'rot. + high-pass div. (k$^{-3}$), along-track',linewidth=3.)
p3, = ax1.loglog(ku,Eun2,'--',color=color3,label=u'rot. + 0.25 x high-pass div. (k$^{-2}$), across-track',linewidth=3.)
ax1.loglog(kv,Evn2,'--',color=color6,label=u'rot. + 0.25 x high-pass div. (k$^{-2}$), along-track',linewidth=3.)
ax1.axis((1./(1000),1./5.,1.e-6,10.))

ax1.loglog(ks2,Es2,'--', color='k',linewidth=2.,alpha=.5)
ax1.loglog(ks3,Es3,'--', color='k',linewidth=2.,alpha=.5)

ax1.text(0.0023277454363699311, 4.2296271323591785,u'k$^{-3}$')
ax1.text(0.0011367633953758597, 0.35651953657755464,u'k$^{-2}$')
ax1.set_ylabel('KE spectral density  [m$^2$ s$^{-2}$/cpkm]')
ax1.set_xlabel('Along-track wavenumber  [cpkm]')
lg = plt.legend(loc=(.24,.69),title= u'', prop={'size':18,}, numpoints=1)
lg.draw_frame(False)
my.leg_width(lg,5.)

ax1.axis((1./1.e3,1.,1./1.e6,1.e1))

ax1.text(0.75, 5., "a", fontsize=32)

add_second_axis(ax1)

ax1.text(.5/1.e2, 5./1.e6, "Synthetic data", size=25, rotation=0.,
         ha="center", va="center",
         bbox = dict(boxstyle="round",ec='k',fc='w'))

figtit = 'figs/EuEv_synthetic.png'
plt.savefig(figtit,format='png', bbox_inches='tight')



In [13]:
# nd + high-pass d: BCF decomposition
fig = plt.figure(facecolor='w', figsize=(11.,12.))
ax1 = fig.add_subplot(111)
p1,= ax1.loglog(kK,Kpsi_nd,color=color1,label=u'rotational',linewidth=3.)
ax1.loglog(kK,Kphi_nd,color=color2,label=u'divergent',linewidth=3.)
p2, =ax1.loglog(kK,Kpsi_d1,color=color3,label=u'rotational',linewidth=3.)
ax1.loglog(kK,Kphi_d1,color=color6,label=u'divergent',linewidth=3.)
p3, = ax1.loglog(kK,Kpsi_d2,'--',color=color3,label=u'rotational',linewidth=3.)
ax1.loglog(kK,Kphi_d2,'--',color=color6,label=u'divergent',linewidth=3.)
ax1.axis((1./(1000),1./5.,1.e-6,10.))

ax1.loglog(ks2,Es2,'--', color='k',linewidth=2.,alpha=.5)
ax1.loglog(ks3,Es3,'--', color='k',linewidth=2.,alpha=.5)

ax1.text(0.0023277454363699311, 4.2296271323591785,u'k$^{-3}$')
ax1.text(0.0011367633953758597, 0.35651953657755464,u'k$^{-2}$')
ax1.set_ylabel('Spectral density  [m$^2$ s$^{-2}$/cpkm]')
ax1.set_xlabel('Along-track wavenumber  [cpkm]')
#lg = plt.legend(loc=(.6x,.65),title= u'', numpoints=1)
#lg.draw_frame(False)
#my.leg_width(lg,5.)

ax1.text(0.75, 5., "b", fontsize=32)

ax1.axis((1./1.e3,1.,1./1.e6,1.e1))

ax1.text(.5/1.e2, 5./1.e6, "Synthetic data", size=25, rotation=0.,
         ha="center", va="center",
         bbox = dict(boxstyle="round",ec='k',fc='w'))

add_second_axis(ax1)

figtit = 'figs/EuEv_synthetic_bcf.png'
plt.savefig(figtit,format='png', bbox_inches='tight')



In [12]:
fig = plt.figure(facecolor='w', figsize=(12.,12.))
plt.loglog(ku,Eu,color=color1,label=u'across-track',linewidth=3.)
plt.loglog(kv,Ev,color=color2,label=u'along-track',linewidth=3.)

plt.loglog(ku,Eunn,'--',color=color1,linewidth=3.)
plt.loglog(kv,Evnn,'--',color=color2,linewidth=3.)

plt.axis((1./(1000),1./5.,1.e-6,10.))

plt.loglog(ks2,Es2,'--', color='k',linewidth=2.,alpha=.5)
plt.loglog(ks3,Es3,'--', color='k',linewidth=2.,alpha=.5)

plt.text(0.0023277454363699311, 4.2296271323591785,u'k$^{-3}$')
plt.text(0.0011367633953758597, 0.35651953657755464,u'k$^{-2}$')
plt.ylabel('KE spectral density  [m$^2$ s$^{-2}$/ cpkm]')
plt.xlabel('Along-track wavenumber  [cpkm]')
lg = plt.legend(loc=1,title= u'rotational component', numpoints=1)
lg.draw_frame(False)

plt.axis((1./1.e3,1.,1./1.e6,1.e1))

plt.text(.5/1.e2, 5./1.e6, "Synthetic data", size=25, rotation=0.,
         ha="center", va="center",
         bbox = dict(boxstyle="round",ec='k',fc='w'))

figtit = 'figs/EuEv_noise.png'
plt.savefig(figtit,format='png', bbox_inches='tight')



In [ ]: