danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added wave_height_steepness method to TimeSeries class

Browse Source
master
per.andreas.brodtkorb 14 years ago
parent 843ddb90c4
commit e2161e8696
7 changed files with 343 additions and 224 deletions
Show Stats Download Patch File Download Diff File

39
pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
Unescape Escape View File

@ -1,28 +1,6 @@
from scipy import * from scipy import *
from pylab import * from pylab import *
#import wafo.spectrum.models as sm
#Sj = sm.Jonswap()
#S = Sj.tospecdata()
#S.data[0:40] = 0.0
#S.data[100:-1] = 0.0
#Nt = len(S.data)-1
#acf = S.tocovdata(nr=0, nt=Nt)
#S2 = acf.tospecdata()
#S.plot('r')
#S2.plot('b:')
#
#show()
#
#import wafo
#import wafo.objects as wo
#xn = wafo.data.sea()
#ts = wo.mat2timeseries(xn)
#Sest = ts.tospecdata(method='cov')
#Sest.setplotter('semilogy')
#Sest.plot()
#show()
# pyreport -o chapter1.html chapter1.py # pyreport -o chapter1.html chapter1.py
#! CHAPTER1 demonstrates some applications of WAFO #! CHAPTER1 demonstrates some applications of WAFO
@ -42,8 +20,8 @@ from pylab import *
#! Simulation of the sea surface from spectrum #! Simulation of the sea surface from spectrum
#! The following code generates 200 seconds of data sampled with 10Hz from #! The following code generates 200 seconds of data sampled with 10Hz from
#! the Torsethaugen spectrum #! the Torsethaugen spectrum
import wafo.spectrum.models as sm import wafo.spectrum.models as wsm
S = sm.Torsethaugen(Hm0=6, Tp=8); S = wsm.Torsethaugen(Hm0=6, Tp=8);
S1 = S.tospecdata() S1 = S.tospecdata()
S1.plot() S1.plot()
show() show()
@ -103,8 +81,8 @@ Nt = 101; # number of angles
th0 = pi / 2; # primary direction of waves th0 = pi / 2; # primary direction of waves
Sp = 15; # spreading parameter Sp = 15; # spreading parameter
D1 = sm.Spreading(type='cos', theta0=th0, method=None) # frequency independent D1 = wsm.Spreading(type='cos', theta0=th0, method=None) # frequency independent
D12 = sm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent D12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent
SD1 = D1.tospecdata2d(S1) SD1 = D1.tospecdata2d(S1)
SD12 = D12.tospecdata2d(S1) SD12 = D12.tospecdata2d(S1)
@ -165,11 +143,12 @@ show()
#! which damage calculations and fatigue life predictions can be #! which damage calculations and fatigue life predictions can be
#! performed. #! performed.
#! #!
#! The matlab commands below computes the intensity of rainflowcycles for #! The commands below computes the intensity of rainflowcycles for
#! the Gaussian model with spectrum S1 using the Markow approximation. #! the Gaussian model with spectrum S1 using the Markov approximation.
#! The rainflow cycles found in the simulated load signal are shown in the #! The rainflow cycles found in the simulated load signal are shown in the
# figure. #! figure.
#clf
#clf()
#paramu = [-6 6 61]; #paramu = [-6 6 61];
#frfc = spec2cmat(S1, [], 'rfc', [], paramu); #frfc = spec2cmat(S1, [], 'rfc', [], paramu);
#pdfplot(frfc); #pdfplot(frfc);

54
pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
Unescape Escape View File

@ -2,7 +2,7 @@ import numpy as np
from scipy import * from scipy import *
from pylab import * from pylab import *
# pyreport -o chapter1.html chapter1.py # pyreport -o chapter2.html chapter2.py
#! CHAPTER2 Modelling random loads and stochastic waves #! CHAPTER2 Modelling random loads and stochastic waves
#!======================================================= #!=======================================================
@ -18,21 +18,7 @@ from pylab import *
#! process. #! process.
#! #!
#! Some of the commands are edited for fast computation. #! Some of the commands are edited for fast computation.
#! Each set of commands is followed by a 'pause' command.
#! #!
#!
#! Tested on Matlab 5.3, 7.01
#! History
#! Revised by Georg Lindgren sept 2009 for WAFO ver 2.5 on Matlab 7.1
#! Revised pab sept2005
#! Added sections -> easier to evaluate using cellmode evaluation.
#! Revised pab Dec2004
#! Created by GL July 13, 2000
#! from commands used in Chapter 2
#!
#! Section 2.1 Introduction and preliminary analysis #! Section 2.1 Introduction and preliminary analysis
#!==================================================== #!====================================================
#! Example 1: Sea data #! Example 1: Sea data
@ -48,7 +34,6 @@ xx[:,1] -= me
ts = wo.mat2timeseries(xx) ts = wo.mat2timeseries(xx)
tp = ts.turning_points() tp = ts.turning_points()
cc = tp.cycle_pairs() cc = tp.cycle_pairs()
lc = cc.level_crossings() lc = cc.level_crossings()
lc.plot() lc.plot()
@ -61,12 +46,15 @@ show()
#! crossing intensity curve, as follows. #! crossing intensity curve, as follows.
T = xx[:, 0].max() - xx[:, 0].min() T = xx[:, 0].max() - xx[:, 0].min()
f0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency f0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency
print('f0 = %g' % f0)
#! Turningpoints and irregularity factor #! Turningpoints and irregularity factor
#!---------------------------------------- #!----------------------------------------
fm = len(tp.data)/(2*T) #! frequency of maxima fm = len(tp.data) / (2 * T) # frequency of maxima
alfa = f0/fm #! approx Tm24/Tm02 alfa = f0 / fm # approx Tm24/Tm02
print('fm = %g, alpha = %g, ' % (fm, alfa))
#! Visually examine data #! Visually examine data
#!------------------------ #!------------------------
@ -97,7 +85,7 @@ dt = ts.sampling_period()
dcrit = 5 * dt dcrit = 5 * dt
ddcrit = 9.81 / 2 * dt * dt ddcrit = 9.81 / 2 * dt * dt
zcrit = 0 zcrit = 0
inds, indg = wm.findoutliers(xx[:,1],zcrit,dcrit,ddcrit, verbose=True) inds, indg = wm.findoutliers(ts.data, zcrit, dcrit, ddcrit, verbose=True)
#! Section 2.2 Frequency Modeling of Load Histories #! Section 2.2 Frequency Modeling of Load Histories
#!---------------------------------------------------- #!----------------------------------------------------
@ -113,12 +101,12 @@ show()
#! Calculate moments #! Calculate moments
#!------------------- #!-------------------
mom, text = S.moment(nr=4) mom, text = S.moment(nr=4)
[sa, sqrt(mom[0])] print('sigma = %g, m0 = %g' % (sa, sqrt(mom[0])))
#! Section 2.2.1 Random functions in Spectral Domain - Gaussian processes #! Section 2.2.1 Random functions in Spectral Domain - Gaussian processes
#!-------------------------------------------------------------------------- #!--------------------------------------------------------------------------
#! Smoothing of spectral estimate #! Smoothing of spectral estimate
#1---------------------------------- #!----------------------------------
#! By decreasing Lmax the spectrum estimate becomes smoother. #! By decreasing Lmax the spectrum estimate becomes smoother.
clf() clf()
@ -132,7 +120,7 @@ show()
#! Estimated autocovariance #! Estimated autocovariance
#!---------------------------- #!----------------------------
#! Obviously knowing the spectrum one can compute the covariance #! Obviously knowing the spectrum one can compute the covariance
#! function. The following matlab code will compute the covariance for the #! function. The following code will compute the covariance for the
#! unimodal spectral density S1 and compare it with estimated #! unimodal spectral density S1 and compare it with estimated
#! covariance of the signal xx. #! covariance of the signal xx.
clf() clf()
@ -171,24 +159,21 @@ sk, ku=S1.stats_nl(moments='sk' )
clf() clf()
import wafo.transform.models as wtm import wafo.transform.models as wtm
gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata() gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata()
g = wtm.TrLinear(mean=me, sigma=sa ).trdata() g = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation
#! Linear transformation
glc, gemp = lc.trdata(mean=me, sigma=sa) glc, gemp = lc.trdata(mean=me, sigma=sa)
g.plot('r')
glc.plot('b-') #! Transf. estimated from level-crossings glc.plot('b-') #! Transf. estimated from level-crossings
gh.plot('b-.') #! Hermite Transf. estimated from moments gh.plot('b-.') #! Hermite Transf. estimated from moments
g.plot('r')
grid('on') grid('on')
show() show()
#! Test Gaussianity of a stochastic process. #! Test Gaussianity of a stochastic process
#!--------------------------------------------- #!------------------------------------------
#! TESTGAUSSIAN simulates e(g(u)-u) = int (g(u)-u)^2 du for Gaussian processes #! TESTGAUSSIAN simulates e(g(u)-u) = int (g(u)-u)^2 du for Gaussian processes
#! given the spectral density, S. The result is plotted if test0 is given. #! given the spectral density, S. The result is plotted if test0 is given.
#! This is useful for testing if the process X(t) is Gaussian. #! This is useful for testing if the process X(t) is Gaussian.
#! If 95#! of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5#! level. #! If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
#! #!
#! As we see from the figure below: none of the simulated values of test1 is #! As we see from the figure below: none of the simulated values of test1 is
#! above 1.00. Thus the data significantly departs from a Gaussian distribution. #! above 1.00. Thus the data significantly departs from a Gaussian distribution.
@ -196,8 +181,9 @@ clf()
test0 = glc.dist2gauss() test0 = glc.dist2gauss()
#! the following test takes time #! the following test takes time
N = len(xx) N = len(xx)
test1 = S1.testgaussian(ns=N,cases=50,t0=test0) test1 = S1.testgaussian(ns=N, cases=50, test0=test0)
sum(test1>test0)<5 is_gaussian = sum(test1 > test0) < 5
print(is_gaussian)
show() show()
#! Normalplot of data xx #! Normalplot of data xx
@ -223,7 +209,7 @@ show()
#! Directional spectrum #! Directional spectrum
clf() clf()
D = wsm.Spreading('cos2s') D = wsm.Spreading('cos2s')
Sd = D.tospecdata2d(S) Sd = D.tospecdata2d(spec)
Sd.plot() Sd.plot()
show() show()

199
pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
Unescape Escape View File

@ -1,40 +1,39 @@
%% CHAPTER3 Demonstrates distributions of wave characteristics import numpy as np
% from scipy import *
% Chapter3 contains the commands used in Chapter3 in the tutorial. from pylab import *
%
% Some of the commands are edited for fast computation. #! CHAPTER3 Demonstrates distributions of wave characteristics
% Each set of commands is followed by a 'pause' command. #!=============================================================
% #!
#! Chapter3 contains the commands used in Chapter3 in the tutorial.
% Tested on Matlab 5.3 #!
% History #! Some of the commands are edited for fast computation.
% Revised pab sept2005 #!
% Added sections -> easier to evaluate using cellmode evaluation. #! Section 3.2 Estimation of wave characteristics from data
% Revised by pab Feb 2005 #!----------------------------------------------------------
% -updated calls to kdetools+spec2XXpdf programs #! Example 1
% Created by GL July 12, 2000 #!~~~~~~~~~~
% from commands used in Chapter 3, written by IR
%
%% Section 3.2 Estimation of wave characteristics from data
%% Example 1
pstate = 'off';
speed = 'fast'
%speed = 'slow'
xx = load('sea.dat');
xx(:,2) = detrend(xx(:,2));
rate = 8;
Tcrcr = dat2wa(xx,0,'c2c','tw',rate);
Tc = dat2wa(xx,0,'u2d','tw',rate);
disp('Block = 1'), pause(pstate)
%% Histogram of crestperiod compared to the kernel density estimate speed = 'fast'
clf #speed = 'slow'
mean(Tc)
max(Tc) import wafo.data as wd
import wafo.misc as wm
import wafo.objects as wo
xx = wd.sea()
xx[:,1] = wm.detrend(xx[:,1])
ts = wo.mat2timeseries(xx)
Tcrcr, ix = ts.wave_periods(vh=0, pdef='c2c', wdef='tw', rate=8)
Tc, ixc = ts.wave_periods(vh=0, pdef='u2d', wdef='tw', rate=8)
#! Histogram of crestperiod compared to the kernel density estimate
#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import wafo.kdetools as wk
clf()
print(Tc.mean())
print(Tc.max())
t = linspace(0.01,8,200); t = linspace(0.01,8,200);
kopt = kdeoptset('L2',0); kopt = kdeoptset('L2',0);
tic tic
ftc1 = kde(Tc,kopt,t); ftc1 = kde(Tc,kopt,t);
@ -44,7 +43,7 @@ histgrm(Tc,[],[],1)
axis([0 8 0 0.5]) axis([0 8 0 0.5])
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 2'), pause(pstate) disp('Block = 2'), pause(pstate)
%% #!#!
tic tic
kopt.inc = 128; kopt.inc = 128;
ftc2 = kdebin(Tc,kopt); ftc2 = kdebin(Tc,kopt);
@ -54,7 +53,7 @@ title('Kernel Density Estimates')
hold off hold off
disp('Block = 3'), pause(pstate) disp('Block = 3'), pause(pstate)
%% Extreme waves - model check: the highest and steepest wave #!#! Extreme waves - model check: the highest and steepest wave
clf clf
method = 0; method = 0;
rate = 8; rate = 8;
@ -69,9 +68,9 @@ spwaveplot(yn,[indA indS],'k.')
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 5'), pause(pstate) disp('Block = 5'), pause(pstate)
%% Does the highest wave contradict a transformed Gaussian model? #!#! Does the highest wave contradict a transformed Gaussian model?
clf clf
inds1 = (5965:5974)'; % points to remove inds1 = (5965:5974)'; #! points to remove
Nsim = 10; Nsim = 10;
[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ... [y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ...
reconstruct(xx,inds1,Nsim); reconstruct(xx,inds1,Nsim);
@ -87,28 +86,28 @@ wafostamp([],'(ER)')
disp('Block = 6'), disp('Block = 6'),
pause(pstate) pause(pstate)
%% Expected value (solid) compared to data removed #!#! Expected value (solid) compared to data removed
clf clf
plot(xx(inds1,1),xx(inds1,2),'+'), hold on plot(xx(inds1,1),xx(inds1,2),'+'), hold on
mu = tranproc(mu1o,fliplr(grec1)); mu = tranproc(mu1o,fliplr(grec1));
plot(y1(inds1,1), mu), hold off plot(y1(inds1,1), mu), hold off
disp('Block = 7'), pause(pstate) disp('Block = 7'), pause(pstate)
%% Crest height PDF #!#! Crest height PDF
% Transform data so that kde works better #! Transform data so that kde works better
clf clf
L2 = 0.6; L2 = 0.6;
plotnorm(Ac.^L2) plotnorm(Ac.^L2)
%% #!#!
% #!
fac = kde(Ac,{'L2',L2},linspace(0.01,3,200)); fac = kde(Ac,{'L2',L2},linspace(0.01,3,200));
pdfplot(fac) pdfplot(fac)
wafostamp([],'(ER)') wafostamp([],'(ER)')
simpson(fac.x{1},fac.f) simpson(fac.x{1},fac.f)
disp('Block = 8'), pause(pstate) disp('Block = 8'), pause(pstate)
%% Empirical crest height CDF #!#! Empirical crest height CDF
clf clf
Fac = flipud(cumtrapz(fac.x{1},flipud(fac.f))); Fac = flipud(cumtrapz(fac.x{1},flipud(fac.f)));
Fac = [fac.x{1} 1-Fac]; Fac = [fac.x{1} 1-Fac];
@ -117,7 +116,7 @@ axis([0 2 0 1])
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 9'), pause(pstate) disp('Block = 9'), pause(pstate)
%% Empirical crest height CDF compared to a Transformed Rayleigh approximation #!#! Empirical crest height CDF compared to a Transformed Rayleigh approximation
facr = trraylpdf(fac.x{1},'Ac',grec1); facr = trraylpdf(fac.x{1},'Ac',grec1);
Facr = cumtrapz(facr.x{1},facr.f); Facr = cumtrapz(facr.x{1},facr.f);
hold on hold on
@ -126,7 +125,7 @@ axis([1.25 2.25 0.95 1])
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 10'), pause(pstate) disp('Block = 10'), pause(pstate)
%% Joint pdf of crest period and crest amplitude #!#! Joint pdf of crest period and crest amplitude
clf clf
kopt2 = kdeoptset('L2',0.5,'inc',256); kopt2 = kdeoptset('L2',0.5,'inc',256);
Tc = Tcf+Tcb; Tc = Tcf+Tcb;
@ -139,7 +138,7 @@ hold off
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 11'), pause(pstate) disp('Block = 11'), pause(pstate)
%% Example 4: Simple wave characteristics obtained from Jonswap spectrum #!#! Example 4: Simple wave characteristics obtained from Jonswap spectrum
clf clf
S = jonswap([],[5 10]); S = jonswap([],[5 10]);
[m, mt]= spec2mom(S,4,[],0); [m, mt]= spec2mom(S,4,[],0);
@ -150,15 +149,15 @@ spec2bw(S)
[ch Sa2] = spec2char(S,[1 3]) [ch Sa2] = spec2char(S,[1 3])
disp('Block = 13'), pause(pstate) disp('Block = 13'), pause(pstate)
%% Section 3.3.2 Explicit form approximations of wave characteristic densities #!#! Section 3.3.2 Explicit form approximations of wave characteristic densities
%% Longuett-Higgins model for Tc and Ac #!#! Longuett-Higgins model for Tc and Ac
clf clf
t = linspace(0,15,100); t = linspace(0,15,100);
h = linspace(0,6,100); h = linspace(0,6,100);
flh = lh83pdf(t,h,[m(1),m(2),m(3)]); flh = lh83pdf(t,h,[m(1),m(2),m(3)]);
disp('Block = 14'), pause(pstate) disp('Block = 14'), pause(pstate)
%% Transformed Longuett-Higgins model for Tc and Ac #!#! Transformed Longuett-Higgins model for Tc and Ac
clf clf
[sk, ku ]=spec2skew(S); [sk, ku ]=spec2skew(S);
sa = sqrt(m(1)); sa = sqrt(m(1));
@ -166,14 +165,14 @@ gh = hermitetr([],[sa sk ku 0]);
flhg = lh83pdf(t,h,[m(1),m(2),m(3)],gh); flhg = lh83pdf(t,h,[m(1),m(2),m(3)],gh);
disp('Block = 15'), pause(pstate) disp('Block = 15'), pause(pstate)
%% Cavanie model for Tc and Ac #!#! Cavanie model for Tc and Ac
clf clf
t = linspace(0,10,100); t = linspace(0,10,100);
h = linspace(0,7,100); h = linspace(0,7,100);
fcav = cav76pdf(t,h,[m(1) m(2) m(3) m(5)],[]); fcav = cav76pdf(t,h,[m(1) m(2) m(3) m(5)],[]);
disp('Block = 16'), pause(pstate) disp('Block = 16'), pause(pstate)
%% Example 5 Transformed Rayleigh approximation of crest- vs trough- amplitude #!#! Example 5 Transformed Rayleigh approximation of crest- vs trough- amplitude
clf clf
xx = load('sea.dat'); xx = load('sea.dat');
x = xx; x = xx;
@ -194,7 +193,7 @@ hold off
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 17'), pause(pstate) disp('Block = 17'), pause(pstate)
%% #!#!
clf clf
TC = dat2tc(xx, me); TC = dat2tc(xx, me);
tc = tp2mm(TC); tc = tp2mm(TC);
@ -203,7 +202,7 @@ At = -tc(:,1);
AcAt = Ac+At; AcAt = Ac+At;
disp('Block = 18'), pause(pstate) disp('Block = 18'), pause(pstate)
%% #!#!
clf clf
Fac_h = [fac_h.x{1} cumtrapz(fac_h.x{1},fac_h.f)]; Fac_h = [fac_h.x{1} cumtrapz(fac_h.x{1},fac_h.f)];
subplot(3,1,1) subplot(3,1,1)
@ -232,8 +231,8 @@ title('At+Ac CDF')
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 19'), pause(pstate) disp('Block = 19'), pause(pstate)
%% Section 3.4 Exact wave distributions in transformed Gaussian Sea #!#! Section 3.4 Exact wave distributions in transformed Gaussian Sea
%% Section 3.4.1 Density of crest period, crest length or encountered crest period #!#! Section 3.4.1 Density of crest period, crest length or encountered crest period
clf clf
S1 = torsethaugen([],[6 8],1); S1 = torsethaugen([],[6 8],1);
D1 = spreading(101,'cos',pi/2,[15],[],0); D1 = spreading(101,'cos',pi/2,[15],[],0);
@ -242,7 +241,7 @@ SD1 = mkdspec(S1,D1);
SD12 = mkdspec(S1,D12); SD12 = mkdspec(S1,D12);
disp('Block = 20'), pause(pstate) disp('Block = 20'), pause(pstate)
%% Crest period #!#! Crest period
clf clf
tic tic
f_tc = spec2tpdf(S1,[],'Tc',[0 11 56],[],4); f_tc = spec2tpdf(S1,[],'Tc',[0 11 56],[],4);
@ -252,13 +251,13 @@ wafostamp([],'(ER)')
simpson(f_tc.x{1},f_tc.f) simpson(f_tc.x{1},f_tc.f)
disp('Block = 21'), pause(pstate) disp('Block = 21'), pause(pstate)
%% Crest length #!#! Crest length
if strncmpi(speed,'slow',1) if strncmpi(speed,'slow',1)
opt1 = rindoptset('speed',5,'method',3); opt1 = rindoptset('speed',5,'method',3);
opt2 = rindoptset('speed',5,'nit',2,'method',0); opt2 = rindoptset('speed',5,'nit',2,'method',0);
else else
% fast #! fast
opt1 = rindoptset('speed',7,'method',3); opt1 = rindoptset('speed',7,'method',3);
opt2 = rindoptset('speed',7,'nit',2,'method',0); opt2 = rindoptset('speed',7,'nit',2,'method',0);
end end
@ -271,7 +270,7 @@ else
disp('NIT=5 may take time, running with NIT=3 in the following') disp('NIT=5 may take time, running with NIT=3 in the following')
NITa = 3; NITa = 3;
end end
%f_Lc = spec2tpdf2(S1,[],'Lc',[0 200 81],opt1); % Faster and more accurate #!f_Lc = spec2tpdf2(S1,[],'Lc',[0 200 81],opt1); #! Faster and more accurate
f_Lc = spec2tpdf(S1,[],'Lc',[0 200 81],[],NITa); f_Lc = spec2tpdf(S1,[],'Lc',[0 200 81],[],NITa);
pdfplot(f_Lc,'-.') pdfplot(f_Lc,'-.')
wafostamp([],'(ER)') wafostamp([],'(ER)')
@ -279,27 +278,27 @@ disp('Block = 22'), pause(pstate)
f_Lc_1 = spec2tpdf(S1,[],'Lc',[0 200 81],1.5,NITa); f_Lc_1 = spec2tpdf(S1,[],'Lc',[0 200 81],1.5,NITa);
%f_Lc_1 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1); #!f_Lc_1 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1);
hold on hold on
pdfplot(f_Lc_1) pdfplot(f_Lc_1)
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 23'), pause(pstate) disp('Block = 23'), pause(pstate)
%% #!#!
clf clf
simpson(f_Lc.x{1},f_Lc.f) simpson(f_Lc.x{1},f_Lc.f)
simpson(f_Lc_1.x{1},f_Lc_1.f) simpson(f_Lc_1.x{1},f_Lc_1.f)
disp('Block = 24'), pause(pstate) disp('Block = 24'), pause(pstate)
%% #!#!
clf clf
tic tic
f_Lc_d1 = spec2tpdf(rotspec(SD1,pi/2),[],'Lc',[0 300 121],[],NITa); f_Lc_d1 = spec2tpdf(rotspec(SD1,pi/2),[],'Lc',[0 300 121],[],NITa);
f_Lc_d12 = spec2tpdf(SD12,[],'Lc',[0 200 81],[],NITa); f_Lc_d12 = spec2tpdf(SD12,[],'Lc',[0 200 81],[],NITa);
% f_Lc_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1); #! f_Lc_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1);
% f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1); #! f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1);
toc toc
pdfplot(f_Lc_d1,'-.'), hold on pdfplot(f_Lc_d1,'-.'), hold on
pdfplot(f_Lc_d12), hold off pdfplot(f_Lc_d12), hold off
@ -307,7 +306,7 @@ wafostamp([],'(ER)')
disp('Block = 25'), pause(pstate) disp('Block = 25'), pause(pstate)
%% #!#!
clf clf
@ -317,26 +316,26 @@ if strncmpi(speed,'slow',1)
f_Lc_d1_5 = spec2tpdf(SD1r,[], 'Lc',[0 300 121],[],5); f_Lc_d1_5 = spec2tpdf(SD1r,[], 'Lc',[0 300 121],[],5);
pdfplot(f_Lc_d1_5), hold on pdfplot(f_Lc_d1_5), hold on
else else
% fast #! fast
disp('Run the following example only if you want a check on computing time') disp('Run the following example only if you want a check on computing time')
disp('Edit the command file and remove %') disp('Edit the command file and remove #!')
end end
f_Lc_d1_3 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],3); f_Lc_d1_3 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],3);
f_Lc_d1_2 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],2); f_Lc_d1_2 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],2);
f_Lc_d1_0 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],0); f_Lc_d1_0 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],0);
%f_Lc_d1_n4 = spec2tpdf2(SD1r,[],'Lc',[0 400 161],opt1); #!f_Lc_d1_n4 = spec2tpdf2(SD1r,[],'Lc',[0 400 161],opt1);
pdfplot(f_Lc_d1_3), hold on pdfplot(f_Lc_d1_3), hold on
pdfplot(f_Lc_d1_2) pdfplot(f_Lc_d1_2)
pdfplot(f_Lc_d1_0) pdfplot(f_Lc_d1_0)
%pdfplot(f_Lc_d1_n4) #!pdfplot(f_Lc_d1_n4)
%simpson(f_Lc_d1_n4.x{1},f_Lc_d1_n4.f) #!simpson(f_Lc_d1_n4.x{1},f_Lc_d1_n4.f)
disp('Block = 26'), pause(pstate) disp('Block = 26'), pause(pstate)
%% Section 3.4.2 Density of wave period, wave length or encountered wave period #!#! Section 3.4.2 Density of wave period, wave length or encountered wave period
%% Example 7: Crest period and high crest waves #!#! Example 7: Crest period and high crest waves
clf clf
tic tic
xx = load('sea.dat'); xx = load('sea.dat');
@ -354,7 +353,7 @@ t = linspace(0.01,8,200);
ftc1 = kde(Tc,{'L2',0},t); ftc1 = kde(Tc,{'L2',0},t);
pdfplot(ftc1) pdfplot(ftc1)
hold on hold on
% f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,4); #! f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,4);
f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,2); f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,2);
simpson(f_t.x{1},f_t.f) simpson(f_t.x{1},f_t.f)
pdfplot(f_t,'-.') pdfplot(f_t,'-.')
@ -363,7 +362,7 @@ wafostamp([],'(ER)')
toc toc
disp('Block = 27'), pause(pstate) disp('Block = 27'), pause(pstate)
%% #!#!
clf clf
tic tic
@ -372,7 +371,7 @@ if strncmpi(speed,'slow',1)
else else
NIT = 2; NIT = 2;
end end
% f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4); #! f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4);
tic tic
f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],Hs/2,NIT); f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],Hs/2,NIT);
toc toc
@ -390,14 +389,14 @@ wafostamp([],'(ER)')
toc toc
disp('Block = 28'), pause(pstate) disp('Block = 28'), pause(pstate)
%% Example 8: Wave period for high crest waves #!#! Example 8: Wave period for high crest waves
% clf #! clf
tic tic
f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1); f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1);
toc toc
simpson(f_tcc2.x{1},f_tcc2.f) simpson(f_tcc2.x{1},f_tcc2.f)
f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],3,5); f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],3,5);
% f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5); #! f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5);
simpson(f_tcc3.x{1},f_tcc3.f) simpson(f_tcc3.x{1},f_tcc3.f)
pdfplot(f_tcc2,'-.') pdfplot(f_tcc2,'-.')
hold on hold on
@ -406,7 +405,7 @@ toc
toc toc
disp('Block = 29'), pause(pstate) disp('Block = 29'), pause(pstate)
%% #!#!
clf clf
[TC tc_ind v_ind] = dat2tc(yn,[],'dw'); [TC tc_ind v_ind] = dat2tc(yn,[],'dw');
N = length(tc_ind); N = length(tc_ind);
@ -418,7 +417,7 @@ spwaveplot(yn,ind(2:4))
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 30'), pause(pstate) disp('Block = 30'), pause(pstate)
%% #!#!
clf clf
Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1); Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1);
t = linspace(0.01,14,200); t = linspace(0.01,14,200);
@ -441,10 +440,10 @@ disp('Block = 32'), pause(pstate)
disp('The rest of this chapter deals with joint densities.') disp('The rest of this chapter deals with joint densities.')
disp('Some calculations may take some time.') disp('Some calculations may take some time.')
disp('You could experiment with other NIT.') disp('You could experiment with other NIT.')
%return #!return
%% Section 3.4.3 Joint density of crest period and crest height #!#! Section 3.4.3 Joint density of crest period and crest height
%% Example 9. Some preliminary analysis of the data #!#! Example 9. Some preliminary analysis of the data
clf clf
tic tic
yy = load('gfaksr89.dat'); yy = load('gfaksr89.dat');
@ -457,7 +456,7 @@ v = v(2)
toc toc
disp('Block = 33'), pause(pstate) disp('Block = 33'), pause(pstate)
%% #!#!
clf clf
tic tic
[TC, tc_ind, v_ind] = dat2tc(yy,v,'dw'); [TC, tc_ind, v_ind] = dat2tc(yy,v,'dw');
@ -477,7 +476,7 @@ At = v-yy(t_ind,2);
toc toc
disp('Block = 34'), pause(pstate) disp('Block = 34'), pause(pstate)
%% #!#!
clf clf
tic tic
t = linspace(0.01,15,200); t = linspace(0.01,15,200);
@ -497,8 +496,8 @@ wafostamp([],'(ER)')
toc toc
disp('Block = 35'), pause(pstate) disp('Block = 35'), pause(pstate)
%% Example 10: Joint characteristics of a half wave: #!#! Example 10: Joint characteristics of a half wave:
%% position and height of a crest for a wave with given period #!#! position and height of a crest for a wave with given period
clf clf
tic tic
ind = find(4.4<Tc & Tc<4.6); ind = find(4.4<Tc & Tc<4.6);
@ -510,7 +509,7 @@ wafostamp([],'(ER)')
toc toc
disp('Block = 36'), pause(pstate) disp('Block = 36'), pause(pstate)
%% #!#!
clf clf
tic tic
opt1 = rindoptset('speed',5,'method',3); opt1 = rindoptset('speed',5,'method',3);
@ -532,7 +531,7 @@ toc
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 37'), pause(pstate) disp('Block = 37'), pause(pstate)
%% #!#!
clf clf
f_tcac_s = spec2thpdf(SS,[],'TcAc',[0 12 81],[Hs/2:0.1:2*Hs],opt1); f_tcac_s = spec2thpdf(SS,[],'TcAc',[0 12 81],[Hs/2:0.1:2*Hs],opt1);
disp('Block = 38'), pause(pstate) disp('Block = 38'), pause(pstate)
@ -553,16 +552,16 @@ toc
wafostamp([],'(ER)') wafostamp([],'(ER)')
disp('Block = 39'), pause(pstate) disp('Block = 39'), pause(pstate)
%% #!#!
clf clf
% f_tcac = spec2thpdf(SS,[],'TcAc',[0 12 81],[0:0.2:8],opt1); #! f_tcac = spec2thpdf(SS,[],'TcAc',[0 12 81],[0:0.2:8],opt1);
% pdfplot(f_tcac) #! pdfplot(f_tcac)
disp('Block = 40'), pause(pstate) disp('Block = 40'), pause(pstate)
%% Section 3.4.4 Joint density of crest and trough height #!#! Section 3.4.4 Joint density of crest and trough height
%% Section 3.4.5 Min-to-max distributions – Markov method #!#! Section 3.4.5 Min-to-max distributions <20> Markov method
%% Example 11. (min-max problems with Gullfaks data) #!#! Example 11. (min-max problems with Gullfaks data)
%% Joint density of maximum and the following minimum #!#! Joint density of maximum and the following minimum
clf clf
tic tic
tp = dat2tp(yy); tp = dat2tp(yy);
@ -578,7 +577,7 @@ wafostamp([],'(ER)')
toc toc
disp('Block = 41'), pause(pstate) disp('Block = 41'), pause(pstate)
%% The joint density of ”still water separated” maxima and minima. #!#! The joint density of <20>still water separated<65> maxima and minima.
clf clf
tic tic
ind = find(Mm(:,1)>v & Mm(:,2)<v); ind = find(Mm(:,1)>v & Mm(:,2)<v);
@ -595,7 +594,7 @@ toc
disp('Block = 42'), pause(pstate) disp('Block = 42'), pause(pstate)
%% #!#!
clf clf
tic tic
facat = kde([Ac At]); facat = kde([Ac At]);

2
pywafo/src/wafo/misc.py
Unescape Escape View File

@ -199,7 +199,7 @@ def detrendma(x, L):
y[n - L::] = x1[n - L::] - trend[-1] y[n - L::] = x1[n - L::] - trend[-1]
return y return y
def ecross(t, f, ind, v): def ecross(t, f, ind, v=0):
''' '''
Extracts exact level v crossings Extracts exact level v crossings

229
pywafo/src/wafo/objects.py
Unescape Escape View File

@ -16,6 +16,9 @@ from __future__ import division
from wafo.transform.core import TrData from wafo.transform.core import TrData
from wafo.transform.models import TrHermite, TrOchi, TrLinear from wafo.transform.models import TrHermite, TrOchi, TrLinear
from wafo.stats import edf from wafo.stats import edf
from wafo.misc import (nextpow2, findtp, findtc, findcross,
ecross, JITImport, DotDict, gravity)
from wafodata import WafoData
from wafo.interpolate import SmoothSpline from wafo.interpolate import SmoothSpline
from scipy.interpolate.interpolate import interp1d from scipy.interpolate.interpolate import interp1d
from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport
@ -26,7 +29,7 @@ import numpy as np
from numpy import (inf, pi, zeros, ones, sqrt, where, log, exp, sin, arcsin, mod, finfo, interp, #@UnresolvedImport from numpy import (inf, pi, zeros, ones, sqrt, where, log, exp, sin, arcsin, mod, finfo, interp, #@UnresolvedImport
linspace, arange, sort, all, abs, vstack, hstack, atleast_1d, #@UnresolvedImport linspace, arange, sort, all, abs, vstack, hstack, atleast_1d, #@UnresolvedImport
finfo, polyfit, r_, nonzero, cumsum, ravel, size, isnan, nan, floor, ceil, diff, array) #@UnresolvedImport finfo, polyfit, r_, nonzero, cumsum, ravel, size, isnan, nan, ceil, diff, array) #@UnresolvedImport
from numpy.fft import fft from numpy.fft import fft
from numpy.random import randn from numpy.random import randn
from scipy.integrate import trapz from scipy.integrate import trapz
@ -34,15 +37,14 @@ from pylab import stineman_interp
from matplotlib.mlab import psd, detrend_mean from matplotlib.mlab import psd, detrend_mean
import scipy.signal import scipy.signal
from wafo.misc import (nextpow2, findtp, findtc, findcross,
ecross, JITImport, DotDict)
from wafodata import WafoData
from plotbackend import plotbackend from plotbackend import plotbackend
import matplotlib import matplotlib
from scipy.stats.stats import skew, kurtosis from scipy.stats.stats import skew, kurtosis
from scipy.signal.windows import parzen from scipy.signal.windows import parzen
from scipy import special from scipy import special
floatinfo = finfo(float) floatinfo = finfo(float)
matplotlib.interactive(True) matplotlib.interactive(True)
_wafocov = JITImport('wafo.covariance') _wafocov = JITImport('wafo.covariance')
@ -64,6 +66,7 @@ class LevelCrossings(WafoData):
number of upcrossings number of upcrossings
args : array-like args : array-like
crossing levels crossing levels
Examples Examples
-------- --------
>>> import wafo.data >>> import wafo.data
@ -80,29 +83,30 @@ class LevelCrossings(WafoData):
def __init__(self, *args, **kwds): def __init__(self, *args, **kwds):
options = dict(title='Level crossing spectrum', options = dict(title='Level crossing spectrum',
xlab='Levels', ylab='Count', xlab='Levels', ylab='Count',
plotmethod='semilogy',
plot_args=['b'], plot_args=['b'],
plot_args_children=['r--'],) plot_args_children=['r--'],)
options.update(**kwds) options.update(**kwds)
super(LevelCrossings, self).__init__(*args, **options) super(LevelCrossings, self).__init__(*args, **options)
self.stdev = kwds.get('stdev', None) self.sigma = kwds.get('sigma', None)
self.mean = kwds.get('mean', None) self.mean = kwds.get('mean', None)
self.setplotter(plotmethod='step') #self.setplotter(plotmethod='step')
icmax = self.data.argmax() icmax = self.data.argmax()
if self.data != None: if self.data != None:
if self.stdev is None or self.mean is None: if self.sigma is None or self.mean is None:
logcros = where(self.data == 0.0, inf, -log(self.data)) logcros = where(self.data == 0.0, inf, -log(self.data))
logcmin = logcros[icmax] logcmin = logcros[icmax]
logcros = sqrt(2 * abs(logcros - logcmin)) logcros = sqrt(2 * abs(logcros - logcmin))
logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1] logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1]
p = polyfit(self.args[10:-9], logcros[10:-9], 1) #least square fit p = polyfit(self.args[10:-9], logcros[10:-9], 1) #least square fit
if self.stdev is None: if self.sigma is None:
self.stdev = 1.0 / p[0] #estimated standard deviation of x self.sigma = 1.0 / p[0] #estimated standard deviation of x
if self.mean is None: if self.mean is None:
self.mean = -p[1] / p[0] #self.args[icmax] self.mean = -p[1] / p[0] #self.args[icmax]
cmax = self.data[icmax] cmax = self.data[icmax]
x = (self.args - self.mean) / self.stdev x = (self.args - self.mean) / self.sigma
y = cmax * exp(-x ** 2 / 2.0) y = cmax * exp(-x ** 2 / 2.0)
self.children = [WafoData(y, self.args)] self.children = [WafoData(y, self.args)]
@ -212,7 +216,7 @@ class LevelCrossings(WafoData):
scr = trapz(lc1 ** 2 * lc, lc1) scr = trapz(lc1 ** 2 * lc, lc1)
scr = sqrt(scr - mcr ** 2) scr = sqrt(scr - mcr ** 2)
lc2 = LevelCrossings(lc, lc1, mean=mcr, stdev=scr) lc2 = LevelCrossings(lc, lc1, mean=mcr, sigma=scr)
g = lc2.trdata() g = lc2.trdata()
@ -350,7 +354,7 @@ class LevelCrossings(WafoData):
if mean is None: if mean is None:
mean = self.mean mean = self.mean
if sigma is None: if sigma is None:
sigma = self.stdev sigma = self.sigma
opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05, opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05,
param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000, ne=7, gvar=1) param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000, ne=7, gvar=1)
@ -401,7 +405,7 @@ class LevelCrossings(WafoData):
inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over
slc22 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu) slc22 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu)
g = TrData(slc22.copy(), g1.copy(), mean=mean, sigma=sigma) #*sa; #multiply with stdev g = TrData(slc22.copy(), g1.copy(), mean=mean, sigma=sigma)
if opt.chkder: if opt.chkder:
for ix in range(5): for ix in range(5):
@ -446,7 +450,7 @@ class CyclePairs(WafoData):
''' '''
def __init__(self, *args, **kwds): def __init__(self, *args, **kwds):
self.kind = kwds.pop('kind', 'max2min') self.kind = kwds.pop('kind', 'max2min')
self.stdev = kwds.pop('stdev', None) self.sigma = kwds.pop('sigma', None)
self.mean = kwds.pop('mean', None) self.mean = kwds.pop('mean', None)
options = dict(title=self.kind + ' cycle pairs', options = dict(title=self.kind + ' cycle pairs',
@ -590,7 +594,7 @@ class CyclePairs(WafoData):
dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx] dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx]
elif defnr == 3: ## This are upcrossings + minima + maxima elif defnr == 3: ## This are upcrossings + minima + maxima
dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx] dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx]
return LevelCrossings(dcount, levels, mean=self.mean, stdev=self.stdev) return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma)
class TurningPoints(WafoData): class TurningPoints(WafoData):
''' '''
@ -613,7 +617,7 @@ class TurningPoints(WafoData):
''' '''
def __init__(self, *args, **kwds): def __init__(self, *args, **kwds):
self.name_ = kwds.pop('name', 'WAFO TurningPoints Object') self.name_ = kwds.pop('name', 'WAFO TurningPoints Object')
self.stdev = kwds.pop('stdev', None) self.sigma = kwds.pop('sigma', None)
self.mean = kwds.pop('mean', None) self.mean = kwds.pop('mean', None)
options = dict(title='Turning points') options = dict(title='Turning points')
@ -672,7 +676,7 @@ class TurningPoints(WafoData):
kind = 'max2min' kind = 'max2min'
M = self.data[iM:-1:2] M = self.data[iM:-1:2]
m = self.data[iM + 1::2] m = self.data[iM + 1::2]
return CyclePairs(M, m, kind=kind, mean=self.mean, stdev=self.stdev) return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma)
def mat2timeseries(x): def mat2timeseries(x):
""" """
@ -772,7 +776,7 @@ class TimeSeries(WafoData):
with attributes: with attributes:
data : ACF vector length L+1 data : ACF vector length L+1
args : time lags length L+1 args : time lags length L+1
stdev : estimated large lag standard deviation of the estimate sigma : estimated large lag standard deviation of the estimate
assuming x is a Gaussian process: assuming x is a Gaussian process:
if R(k)=0 for all lags k>q then an approximation if R(k)=0 for all lags k>q then an approximation
of the variance for large samples due to Bartlett of the variance for large samples due to Bartlett
@ -828,14 +832,15 @@ class TimeSeries(WafoData):
t = linspace(0, lag * dt, lag + 1) t = linspace(0, lag * dt, lag + 1)
#cumsum = np.cumsum #cumsum = np.cumsum
acf = _wafocov.CovData1D(R[lags], t) acf = _wafocov.CovData1D(R[lags], t)
acf.stdev = sqrt(r_[ 0, r0**2 , r0**2 + 2 * cumsum(R[1:] ** 2)] / Ncens) acf.sigma = sqrt(r_[ 0, r0**2 , r0**2 + 2 * cumsum(R[1:] ** 2)] / Ncens)
acf.children = [WafoData(-2. * acf.stdev[lags], t), WafoData(2. * acf.stdev[lags], t)] acf.children = [WafoData(-2. * acf.sigma[lags], t), WafoData(2. * acf.sigma[lags], t)]
acf.plot_args_children = ['r:'] acf.plot_args_children = ['r:']
acf.norm = norm acf.norm = norm
return acf return acf
def specdata(self, L=None, tr=None, method='cov', detrend=detrend_mean, window=parzen, noverlap=0, pad_to=None): def _specdata(self, L=None, tr=None, method='cov', detrend=detrend_mean, window=parzen, noverlap=0, pad_to=None):
""" """
Obsolete: Delete?
Return power spectral density by Welches average periodogram method. Return power spectral density by Welches average periodogram method.
Parameters Parameters
@ -963,7 +968,7 @@ class TimeSeries(WafoData):
R = tsy.tocovdata() R = tsy.tocovdata()
if estimate_L: if estimate_L:
#finding where ACF is less than 2 st. deviations. #finding where ACF is less than 2 st. deviations.
L = max_L + 2 - (np.abs(R.data[max_L::-1]) > 2 * R.stdev[max_L::-1]).argmax() # a better L value L = max_L + 2 - (np.abs(R.data[max_L::-1]) > 2 * R.sigma[max_L::-1]).argmax() # a better L value
if wdef == 1: # modify L so that hanning and Parzen give appr. the same result if wdef == 1: # modify L so that hanning and Parzen give appr. the same result
L = min(int(4 * L / 3), n - 2) L = min(int(4 * L / 3), n - 2)
print('The default L is set to %d' % L) print('The default L is set to %d' % L)
@ -1025,6 +1030,157 @@ class TimeSeries(WafoData):
# S.S = zeros(nf+1,m-1); # S.S = zeros(nf+1,m-1);
return spec return spec
def wave_height_steepness(self, rate=1, method=1, g=None):
'''
Returns waveheights and steepnesses from data.
Parameters
----------
rate : scalar integer
interpolation rate. Interpolates with spline if greater than one.
method : scalar integer
0 max(Vcf, Vcb) and corresponding wave height Hd or Hu in H
1 crest front (rise) speed (Vcf) in S and wave height Hd in H. (default)
-1 crest back (fall) speed (Vcb) in S and waveheight Hu in H.
2 crest front steepness in S and the wave height Hd in H.
-2 crest back steepness in S and the wave height Hu in H.
3 total wave steepness in S and the wave height Hd in H
for zero-downcrossing waves.
-3 total wave steepness in S and the wave height Hu in H.
for zero-upcrossing waves.
Returns
-------
S, H = Steepness and the corresponding wave height according to method
The parameters are calculated as follows:
Crest front speed (velocity) = Vcf = Ac/Tcf
Crest back speed (velocity) = Vcb = Ac/Tcb
Crest front steepness = 2*pi*Ac./Td/Tcf/g
Crest back steepness = 2*pi*Ac./Tu/Tcb/g
Total wave steepness (zero-downcrossing wave) = 2*pi*Hd./Td.^2/g
Total wave steepness (zero-upcrossing wave) = 2*pi*Hu./Tu.^2/g
The definition of g, Ac,At, Tcf, etc. are given in gravity and
wafo.definitions.
Example
-------
>>> import wafo.data as wd
>>> import wafo.objects as wo
>>> x = wd.sea()
>>> ts = wo.mat2timeseries(x)
>>> for i in xrange(-3,4):
... S, H = ts.wave_height_steepness(method=i)
... print(S[:2],H[:2])
(array([ 0.01186982, 0.04852534]), array([ 0.69, 0.86]))
(array([ 0.02918363, 0.06385979]), array([ 0.69, 0.86]))
(array([ 0.27797411, 0.33585743]), array([ 0.69, 0.86]))
(array([ 0.60835634, 0.60930197]), array([ 0.42, 0.78]))
(array([ 0.60835634, 0.60930197]), array([ 0.42, 0.78]))
(array([ 0.10140867, 0.06141156]), array([ 0.42, 0.78]))
(array([ 0.01821413, 0.01236672]), array([ 0.42, 0.78]))
>>> import pylab as plt
>>> h = plt.plot(S,H,'.')
>>> h = plt.xlabel('S')
>>> h = plt.ylabel('Hd [m]')
See also
--------
wafo.definitions
'''
# Ac,At = crest and trough amplitude, respectively
# Tcf,
# Tcb = Crest front and crest (rear) back period, respectively
# z_ind = indices to the zero-crossings (d,u) of the defining
# trough to trough waves (tw). If M>1 then
# z_ind=[N1 z_ind1 N2 z_ind2 ...NM z_indM] where
# Ni = length(z_indi) and z_indi are the indices to
# the zero-crossings of xi, i=1,2...M.
#
# yn = interpolated signal
#
# xn = [ti x1 x2 ... xM], where
# ti = time and x1 x2 ... xM are M column vectors of
# sampled surface elevation.
#[S,H,z_ind2,AC1,AT1,TFRONT1,TREAR1]=deal([]); % Initialize to []
dT = self.sampling_period()
#[N M] = size(xx);
if g is None:
g = gravity() #% acceleration of gravity
if rate>1:
dT = dT/rate
t0, tn = self.args[0], self.args[-1]
n = len(self.args)
ti = linspace(t0, tn, int(rate*n))
xi = interp1d(self.args , self.data.ravel(), kind='cubic')(ti)
else:
ti, xi = self.args, self.data.ravel()
#for ix=2:M
tc_ind, z_ind = findtc(xi,v=0,kind='tw')
tc_a = xi[tc_ind]
tc_t = ti[tc_ind]
AC = tc_a[1::2] # crest amplitude
AT= -tc_a[0::2] # trough amplitude
if (0<= method and method <=2):
# time between zero-upcrossing and crest [s]
tu = ecross(ti, xi, z_ind[1:-1:2],v=0)
TFRONT = tc_t[1::2]-tu
TFRONT[(TFRONT==0)]=dT # avoiding division by zero
if (0 >= method and method>=-2):
# time between crest and zero-downcrossing [s]
td = ecross(ti,xi, z_ind[2::2],v=0)
TREAR = td-tc_t[1::2]
TREAR[(TREAR==0)]=dT; #% avoiding division by zero
if method==0:
# max(Vcf, Vcr) and the corresponding wave height Hd or Hu in H
HU = AC+AT[1:]
HD = AC+AT[:-1]
T = np.where(TFRONT<TREAR, TFRONT, TREAR)
S = AC/T
H = np.where(TFRONT<TREAR, HD, HU)
elif method==1: # extracting crest front velocity [m/s] and
# Zero-downcrossing wave height [m]
H = AC+AT[:-1] # Hd
S = AC/TFRONT
elif method== -1: # extracting crest rear velocity [m/s] and
# Zero-upcrossing wave height [m]
H = AC+AT[1:] #Hu
S = AC/TREAR
elif method== 2: #crest front steepness in S and the wave height Hd in H.
H = AC + AT[:-1] #Hd
T = diff(ecross(ti,xi, z_ind[::2],v=0))
S = 2*pi*AC/T/TFRONT/g
elif method== -2: # crest back steepness in S and the wave height Hu in H.
H = AC+AT[1:]
T = diff(ecross(ti,xi, z_ind[1::2],v=0))
S = 2*pi*AC/T/TREAR/g
elif method==3: # total steepness in S and the wave height Hd in H
# for zero-doewncrossing waves.
H = AC+AT[:-1]
T = diff(ecross(ti,xi , z_ind[::2],v=0))# Period zero-downcrossing waves
S = 2*pi*H/T**2/g
elif method== -3: # total steepness in S and the wave height Hu in H for
# zero-upcrossing waves.
H = AC+AT[1::]
T = diff(ecross(ti, xi, z_ind[1::2],v=0))# Period zero-upcrossing waves
S = 2*pi*H/T**2/g
return S, H
def _trdata_cdf(self, **options): def _trdata_cdf(self, **options):
''' '''
Estimate transformation, g, from observed marginal CDF. Estimate transformation, g, from observed marginal CDF.
@ -1293,8 +1449,8 @@ class TimeSeries(WafoData):
except: except:
t = ind t = ind
mean = self.data.mean() mean = self.data.mean()
stdev = self.data.std() sigma = self.data.std()
return TurningPoints(self.data[ind], t, mean=mean, stdev=stdev) return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma)
def trough_crest(self, v=None, wavetype=None): def trough_crest(self, v=None, wavetype=None):
""" """
@ -1323,8 +1479,8 @@ class TimeSeries(WafoData):
except: except:
t = ind t = ind
mean = self.data.mean() mean = self.data.mean()
stdev = self.data.std() sigma = self.data.std()
return TurningPoints(self.data[ind], t, mean=mean, stdev=stdev) return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma)
def wave_periods(self, vh=None, pdef='d2d', wdef=None, index=None, rate=1): def wave_periods(self, vh=None, pdef='d2d', wdef=None, index=None, rate=1):
""" """
@ -1502,7 +1658,7 @@ class TimeSeries(WafoData):
# TODO: finish reconstruct # TODO: finish reconstruct
pass pass
def plot_wave(self, sym1='k.', ts=None, sym2='k+', nfig=None, nsub=None, def plot_wave(self, sym1='k.', ts=None, sym2='k+', nfig=None, nsub=None,
stdev=None, vfact=3): sigma=None, vfact=3):
''' '''
Plots the surface elevation of timeseries. Plots the surface elevation of timeseries.
@ -1520,7 +1676,7 @@ class TimeSeries(WafoData):
changed to nsub=min(6,ceil(nsub/nfig)) changed to nsub=min(6,ceil(nsub/nfig))
nfig : scalar integer nfig : scalar integer
Number of figures. By default nfig=ceil(Nsub/6). Number of figures. By default nfig=ceil(Nsub/6).
stdev : real scalar sigma : real scalar
standard deviation of data. standard deviation of data.
vfact : real scalar vfact : real scalar
how large in stdev the vertical scale should be (default 3) how large in stdev the vertical scale should be (default 3)
@ -1552,22 +1708,22 @@ class TimeSeries(WafoData):
xn2 = ts.data xn2 = ts.data
tn2 = ts.args tn2 = ts.args
if stdev is None: if sigma is None:
stdev = xn[indg].std() sigma = xn[indg].std()
if nsub is None: if nsub is None:
nsub = int(floor(len(xn2) / (2 * nw))) + 1 # about Nw mdc waves in each plot nsub = int(len(xn2) / (2 * nw)) + 1 # about Nw mdc waves in each plot
if nfig is None: if nfig is None:
nfig = int(ceil(nsub / 6)) nfig = int(ceil(nsub / 6))
nsub = min(6, int(ceil(nsub / nfig))) nsub = min(6, int(ceil(nsub / nfig)))
n = len(xn) n = len(xn)
Ns = int(floor(n / (nfig * nsub))) Ns = int(n / (nfig * nsub))
ind = r_[0:Ns] ind = r_[0:Ns]
if all(xn >= 0): if all(xn >= 0):
vscale = [0, 2 * stdev * vfact] vscale = [0, 2 * sigma * vfact]
else: else:
vscale = array([-1, 1]) * vfact * stdev vscale = array([-1, 1]) * vfact * sigma
XlblTxt = 'Time [sec]' XlblTxt = 'Time [sec]'
@ -1580,7 +1736,7 @@ class TimeSeries(WafoData):
dT = 1 / 60 dT = 1 / 60
XlblTxt = 'Time (minutes)' XlblTxt = 'Time (minutes)'
if np.max(abs(xn[indg])) > 5 * stdev: if np.max(abs(xn[indg])) > 5 * sigma:
XlblTxt = XlblTxt + ' (Spurious data since max > 5 std.)' XlblTxt = XlblTxt + ' (Spurious data since max > 5 std.)'
plot = plotbackend.plot plot = plotbackend.plot
@ -1602,7 +1758,7 @@ class TimeSeries(WafoData):
#plotbackend.axis([h_scale*dT, v_scale]) #plotbackend.axis([h_scale*dT, v_scale])
for iy in [-2, 2]: for iy in [-2, 2]:
plot(h_scale * dT, iy * stdev * ones(2), ':') plot(h_scale * dT, iy * sigma * ones(2), ':')
ind = ind + Ns ind = ind + Ns
#end #end
@ -1630,7 +1786,6 @@ class TimeSeries(WafoData):
>>> ts = wafo.objects.mat2timeseries(x[0:500,...]) >>> ts = wafo.objects.mat2timeseries(x[0:500,...])
>>> h = ts.plot_sp_wave(np.r_[6:9,12:18]) >>> h = ts.plot_sp_wave(np.r_[6:9,12:18])
See also See also
-------- --------
plot_wave, findtc plot_wave, findtc

4
pywafo/src/wafo/spectrum/core.py
Unescape Escape View File

@ -2088,9 +2088,9 @@ class SpecData1D(WafoData):
print(' ... be patient this may take a while') print(' ... be patient this may take a while')
rep = int(np.floor(ns * cases / maxsize) + 1) rep = int(ns * cases / maxsize) + 1
Nstep = np.floor(cases / rep); Nstep = int(cases / rep)
acf = self.tocovdata() acf = self.tocovdata()
#R = spec2cov(S); #R = spec2cov(S);

2
pywafo/src/wafo/wafodata.py
Unescape Escape View File

@ -76,7 +76,7 @@ class WafoData(object):
self.labels = AxisLabels(**kwds) self.labels = AxisLabels(**kwds)
if not self.plotter: if not self.plotter:
self.setplotter() self.setplotter(kwds.get('plotmethod', None))
def plot(self, *args, **kwds): def plot(self, *args, **kwds):
tmp = None tmp = None

Write Preview
Loading…
Cancel
Save