This notebook showcases a basic usage of pyspec for computing 2D spectrum and its associated isotropic spectrum. Other featrures such as bin average in log space and confidence limit estimation are also shown.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.signal as signal
import seawater as sw
from pyspec import spectrum as spec
fni = "data/lg0703_nb150.npz"
data = np.load(fni)
fmax = 135
lon,lat = data['lon1'][:fmax],data['lat1'][:fmax]
u, v = data['u1'][0,:fmax],data['v1'][0,:fmax]
urot, vrot = data['u1_rot'][0,:fmax],data['v1_rot'][0,:fmax]
d,ang = sw.dist(lat,lon,units='km')
dist = np.append(0.,np.cumsum(d))
uspec = spec.Spectrum(urot[:].copy(),dt=5)
vspec = spec.Spectrum(vrot[:].copy(),dt=5)
ksp, Pu_sp = signal.welch(urot[:], fs=1/5., window='hanning', nperseg=urot[:].size, noverlap=False,
detrend='constant', scaling='density', axis=-1)
_, Pv_sp = signal.welch(vrot[:], fs=1/5., window='hanning', nperseg=urot[:].size, noverlap=False,
detrend='constant', scaling='density', axis=-1)
fig = plt.figure(figsize=(9,7))
ax = fig.add_subplot(111)
ax.loglog(uspec.f,uspec.spec,'m')
ax.loglog(ksp,2*Pu_sp,'m--')
ax.loglog(vspec.f,vspec.spec)
ax.loglog(ksp,2*Pv_sp,'b--')
ax.set_xlabel('Wavenumber [cpkm]')
ax.set_ylabel(r'KE spectral density [m$^2$ s$^{-1}$/cpkm]')