From e2161e86962accde8387dd7ece4abddd67f72a65 Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Sat, 29 Jan 2011 09:25:42 +0000 Subject: [PATCH] Added wave_height_steepness method to TimeSeries class --- .../src/wafo/doc/tutorial_scripts/chapter1.py | 41 +-- .../src/wafo/doc/tutorial_scripts/chapter2.py | 86 +++---- .../src/wafo/doc/tutorial_scripts/chapter3.py | 199 ++++++++------- pywafo/src/wafo/misc.py | 2 +- pywafo/src/wafo/objects.py | 233 +++++++++++++++--- pywafo/src/wafo/spectrum/core.py | 4 +- pywafo/src/wafo/wafodata.py | 2 +- 7 files changed, 343 insertions(+), 224 deletions(-) diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py index f4d8b19..bbfa013 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py @@ -1,28 +1,6 @@ from scipy 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 #! CHAPTER1 demonstrates some applications of WAFO @@ -42,8 +20,8 @@ from pylab import * #! Simulation of the sea surface from spectrum #! The following code generates 200 seconds of data sampled with 10Hz from #! the Torsethaugen spectrum -import wafo.spectrum.models as sm -S = sm.Torsethaugen(Hm0=6, Tp=8); +import wafo.spectrum.models as wsm +S = wsm.Torsethaugen(Hm0=6, Tp=8); S1 = S.tospecdata() S1.plot() show() @@ -103,8 +81,8 @@ Nt = 101; # number of angles th0 = pi / 2; # primary direction of waves Sp = 15; # spreading parameter -D1 = sm.Spreading(type='cos', theta0=th0, method=None) # frequency independent -D12 = sm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent +D1 = wsm.Spreading(type='cos', theta0=th0, method=None) # frequency independent +D12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent SD1 = D1.tospecdata2d(S1) SD12 = D12.tospecdata2d(S1) @@ -165,11 +143,12 @@ show() #! which damage calculations and fatigue life predictions can be #! performed. #! -#! The matlab commands below computes the intensity of rainflowcycles for -#! the Gaussian model with spectrum S1 using the Markow approximation. +#! The commands below computes the intensity of rainflowcycles for +#! the Gaussian model with spectrum S1 using the Markov approximation. #! The rainflow cycles found in the simulated load signal are shown in the -# figure. -#clf +#! figure. + +#clf() #paramu = [-6 6 61]; #frfc = spec2cmat(S1, [], 'rfc', [], paramu); #pdfplot(frfc); @@ -195,7 +174,7 @@ ylabel('(m)') #! Formation of 5 min maxima yura = xn[:85500, 1] -yura = np.reshape(yura, (285,300)).T +yura = np.reshape(yura, (285, 300)).T maxyura = yura.max(axis=0) subplot(212) plot(xn[299:85500:300, 0] / 3600, maxyura, '.') diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py index e9d4892..28e1456 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py @@ -2,7 +2,7 @@ import numpy as np from scipy import * from pylab import * -# pyreport -o chapter1.html chapter1.py +# pyreport -o chapter2.html chapter2.py #! CHAPTER2 Modelling random loads and stochastic waves #!======================================================= @@ -18,21 +18,7 @@ from pylab import * #! process. #! #! 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 #!==================================================== #! Example 1: Sea data @@ -42,13 +28,12 @@ from pylab import * import wafo import wafo.objects as wo xx = wafo.data.sea() -me = xx[:,1].mean() -sa = xx[:,1].std() -xx[:,1] -= me +me = xx[:, 1].mean() +sa = xx[:, 1].std() +xx[:, 1] -= me ts = wo.mat2timeseries(xx) tp = ts.turning_points() - cc = tp.cycle_pairs() lc = cc.level_crossings() lc.plot() @@ -59,15 +44,18 @@ show() #! Next we compute the mean frequency as the average number of upcrossings #! per time unit of the mean level (= 0); this may require interpolation in the #! crossing intensity curve, as follows. -T = xx[:,0].max()-xx[:,0].min() -f0 = np.interp(0,lc.args,lc.data,0)/T #! zero up-crossing frequency +T = xx[:, 0].max() - xx[:, 0].min() +f0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency +print('f0 = %g' % f0) #! Turningpoints and irregularity factor #!---------------------------------------- -fm = len(tp.data)/(2*T) #! frequency of maxima -alfa = f0/fm #! approx Tm24/Tm02 - +fm = len(tp.data) / (2 * T) # frequency of maxima +alfa = f0 / fm # approx Tm24/Tm02 + +print('fm = %g, alpha = %g, ' % (fm, alfa)) + #! Visually examine data #!------------------------ #! We finish this section with some remarks about the quality @@ -79,7 +67,7 @@ alfa = f0/fm #! approx Tm24/Tm02 #! First we shall plot the data and zoom in on a specific region. #! A part of sea data is visualized with the following commands clf() -ts.plot_wave('k-',tp,'*',nfig=1, nsub=1) +ts.plot_wave('k-', tp, '*', nfig=1, nsub=1) axis([0, 2, -2, 2]) show() @@ -94,10 +82,10 @@ show() import wafo.misc as wm dt = ts.sampling_period() # dt = np.diff(xx[:2,0]) -dcrit = 5*dt -ddcrit = 9.81/2*dt*dt +dcrit = 5 * dt +ddcrit = 9.81 / 2 * dt * dt 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 #!---------------------------------------------------- @@ -112,13 +100,13 @@ show() #! Calculate moments #!------------------- -mom, text= S.moment(nr=4) -[sa, sqrt(mom[0])] +mom, text = S.moment(nr=4) +print('sigma = %g, m0 = %g' % (sa, sqrt(mom[0]))) #! Section 2.2.1 Random functions in Spectral Domain - Gaussian processes #!-------------------------------------------------------------------------- #! Smoothing of spectral estimate -#1---------------------------------- +#!---------------------------------- #! By decreasing Lmax the spectrum estimate becomes smoother. clf() @@ -132,7 +120,7 @@ show() #! Estimated autocovariance #!---------------------------- #! 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 #! covariance of the signal xx. clf() @@ -141,7 +129,7 @@ R1 = S1.tocovdata(nr=1) Rest = ts.tocovdata(lag=Lmax) R1.plot('.') Rest.plot() -axis([0,25,-0.1,0.25]) +axis([0, 25, -0.1, 0.25]) show() #! We can see in Figure below that the covariance function corresponding to @@ -162,33 +150,30 @@ show() #! for the data set xx and compare it with the second order wave approximation #! proposed by Winterstein: import wafo.stats as ws -rho3 = ws.skew(xx[:,1]) -rho4 = ws.kurtosis(xx[:,1]) +rho3 = ws.skew(xx[:, 1]) +rho4 = ws.kurtosis(xx[:, 1]) -sk, ku=S1.stats_nl(moments='sk' ) +sk, ku = S1.stats_nl(moments='sk') #! Comparisons of 3 transformations clf() import wafo.transform.models as wtm -gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku ).trdata() -g = wtm.TrLinear(mean=me, sigma=sa ).trdata() - - -#! Linear transformation - +gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata() +g = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation glc, gemp = lc.trdata(mean=me, sigma=sa) -g.plot('r') + 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') 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 #! 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. -#! 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 #! above 1.00. Thus the data significantly departs from a Gaussian distribution. @@ -196,8 +181,9 @@ clf() test0 = glc.dist2gauss() #! the following test takes time N = len(xx) -test1 = S1.testgaussian(ns=N,cases=50,t0=test0) -sum(test1>test0)<5 +test1 = S1.testgaussian(ns=N, cases=50, test0=test0) +is_gaussian = sum(test1 > test0) < 5 +print(is_gaussian) show() #! Normalplot of data xx @@ -223,7 +209,7 @@ show() #! Directional spectrum clf() D = wsm.Spreading('cos2s') -Sd = D.tospecdata2d(S) +Sd = D.tospecdata2d(spec) Sd.plot() show() diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py index 34dd4ff..541586b 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py @@ -1,40 +1,39 @@ -%% CHAPTER3 Demonstrates distributions of wave characteristics -% -% Chapter3 contains the commands used in Chapter3 in the tutorial. -% -% Some of the commands are edited for fast computation. -% Each set of commands is followed by a 'pause' command. -% - -% Tested on Matlab 5.3 -% History -% Revised pab sept2005 -% Added sections -> easier to evaluate using cellmode evaluation. -% Revised by pab Feb 2005 -% -updated calls to kdetools+spec2XXpdf programs -% 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'; +import numpy as np +from scipy import * +from pylab import * + +#! CHAPTER3 Demonstrates distributions of wave characteristics +#!============================================================= +#! +#! Chapter3 contains the commands used in Chapter3 in the tutorial. +#! +#! Some of the commands are edited for fast computation. +#! +#! Section 3.2 Estimation of wave characteristics from data +#!---------------------------------------------------------- +#! Example 1 +#!~~~~~~~~~~ + 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 -clf -mean(Tc) -max(Tc) +#speed = 'slow' + +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); + kopt = kdeoptset('L2',0); tic ftc1 = kde(Tc,kopt,t); @@ -44,7 +43,7 @@ histgrm(Tc,[],[],1) axis([0 8 0 0.5]) wafostamp([],'(ER)') disp('Block = 2'), pause(pstate) -%% +#!#! tic kopt.inc = 128; ftc2 = kdebin(Tc,kopt); @@ -54,7 +53,7 @@ title('Kernel Density Estimates') hold off disp('Block = 3'), pause(pstate) -%% Extreme waves - model check: the highest and steepest wave +#!#! Extreme waves - model check: the highest and steepest wave clf method = 0; rate = 8; @@ -69,9 +68,9 @@ spwaveplot(yn,[indA indS],'k.') wafostamp([],'(ER)') disp('Block = 5'), pause(pstate) -%% Does the highest wave contradict a transformed Gaussian model? +#!#! Does the highest wave contradict a transformed Gaussian model? clf -inds1 = (5965:5974)'; % points to remove +inds1 = (5965:5974)'; #! points to remove Nsim = 10; [y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ... reconstruct(xx,inds1,Nsim); @@ -87,28 +86,28 @@ wafostamp([],'(ER)') disp('Block = 6'), pause(pstate) -%% Expected value (solid) compared to data removed +#!#! Expected value (solid) compared to data removed clf plot(xx(inds1,1),xx(inds1,2),'+'), hold on mu = tranproc(mu1o,fliplr(grec1)); plot(y1(inds1,1), mu), hold off disp('Block = 7'), pause(pstate) -%% Crest height PDF -% Transform data so that kde works better +#!#! Crest height PDF +#! Transform data so that kde works better clf L2 = 0.6; plotnorm(Ac.^L2) -%% -% +#!#! +#! fac = kde(Ac,{'L2',L2},linspace(0.01,3,200)); pdfplot(fac) wafostamp([],'(ER)') simpson(fac.x{1},fac.f) disp('Block = 8'), pause(pstate) -%% Empirical crest height CDF +#!#! Empirical crest height CDF clf Fac = flipud(cumtrapz(fac.x{1},flipud(fac.f))); Fac = [fac.x{1} 1-Fac]; @@ -117,7 +116,7 @@ axis([0 2 0 1]) wafostamp([],'(ER)') 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 = cumtrapz(facr.x{1},facr.f); hold on @@ -126,7 +125,7 @@ axis([1.25 2.25 0.95 1]) wafostamp([],'(ER)') disp('Block = 10'), pause(pstate) -%% Joint pdf of crest period and crest amplitude +#!#! Joint pdf of crest period and crest amplitude clf kopt2 = kdeoptset('L2',0.5,'inc',256); Tc = Tcf+Tcb; @@ -139,7 +138,7 @@ hold off wafostamp([],'(ER)') disp('Block = 11'), pause(pstate) -%% Example 4: Simple wave characteristics obtained from Jonswap spectrum +#!#! Example 4: Simple wave characteristics obtained from Jonswap spectrum clf S = jonswap([],[5 10]); [m, mt]= spec2mom(S,4,[],0); @@ -150,15 +149,15 @@ spec2bw(S) [ch Sa2] = spec2char(S,[1 3]) disp('Block = 13'), pause(pstate) -%% Section 3.3.2 Explicit form approximations of wave characteristic densities -%% Longuett-Higgins model for Tc and Ac +#!#! Section 3.3.2 Explicit form approximations of wave characteristic densities +#!#! Longuett-Higgins model for Tc and Ac clf t = linspace(0,15,100); h = linspace(0,6,100); flh = lh83pdf(t,h,[m(1),m(2),m(3)]); disp('Block = 14'), pause(pstate) -%% Transformed Longuett-Higgins model for Tc and Ac +#!#! Transformed Longuett-Higgins model for Tc and Ac clf [sk, ku ]=spec2skew(S); 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); disp('Block = 15'), pause(pstate) -%% Cavanie model for Tc and Ac +#!#! Cavanie model for Tc and Ac clf t = linspace(0,10,100); h = linspace(0,7,100); fcav = cav76pdf(t,h,[m(1) m(2) m(3) m(5)],[]); 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 xx = load('sea.dat'); x = xx; @@ -194,7 +193,7 @@ hold off wafostamp([],'(ER)') disp('Block = 17'), pause(pstate) -%% +#!#! clf TC = dat2tc(xx, me); tc = tp2mm(TC); @@ -203,7 +202,7 @@ At = -tc(:,1); AcAt = Ac+At; disp('Block = 18'), pause(pstate) -%% +#!#! clf Fac_h = [fac_h.x{1} cumtrapz(fac_h.x{1},fac_h.f)]; subplot(3,1,1) @@ -232,8 +231,8 @@ title('At+Ac CDF') wafostamp([],'(ER)') disp('Block = 19'), pause(pstate) -%% 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 Exact wave distributions in transformed Gaussian Sea +#!#! Section 3.4.1 Density of crest period, crest length or encountered crest period clf S1 = torsethaugen([],[6 8],1); D1 = spreading(101,'cos',pi/2,[15],[],0); @@ -242,7 +241,7 @@ SD1 = mkdspec(S1,D1); SD12 = mkdspec(S1,D12); disp('Block = 20'), pause(pstate) -%% Crest period +#!#! Crest period clf tic f_tc = spec2tpdf(S1,[],'Tc',[0 11 56],[],4); @@ -252,13 +251,13 @@ wafostamp([],'(ER)') simpson(f_tc.x{1},f_tc.f) disp('Block = 21'), pause(pstate) -%% Crest length +#!#! Crest length if strncmpi(speed,'slow',1) opt1 = rindoptset('speed',5,'method',3); opt2 = rindoptset('speed',5,'nit',2,'method',0); else - % fast + #! fast opt1 = rindoptset('speed',7,'method',3); opt2 = rindoptset('speed',7,'nit',2,'method',0); end @@ -271,7 +270,7 @@ else disp('NIT=5 may take time, running with NIT=3 in the following') NITa = 3; 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); pdfplot(f_Lc,'-.') 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 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1); +#!f_Lc_1 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1); hold on pdfplot(f_Lc_1) wafostamp([],'(ER)') disp('Block = 23'), pause(pstate) -%% +#!#! clf simpson(f_Lc.x{1},f_Lc.f) simpson(f_Lc_1.x{1},f_Lc_1.f) disp('Block = 24'), pause(pstate) -%% +#!#! clf tic 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_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1); -% f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1); +#! f_Lc_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1); +#! f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1); toc pdfplot(f_Lc_d1,'-.'), hold on pdfplot(f_Lc_d12), hold off @@ -307,7 +306,7 @@ wafostamp([],'(ER)') disp('Block = 25'), pause(pstate) -%% +#!#! clf @@ -317,26 +316,26 @@ if strncmpi(speed,'slow',1) f_Lc_d1_5 = spec2tpdf(SD1r,[], 'Lc',[0 300 121],[],5); pdfplot(f_Lc_d1_5), hold on else - % fast + #! fast 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 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_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_2) 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) -%% Section 3.4.2 Density of wave period, wave length or encountered wave period -%% Example 7: Crest period and high crest waves +#!#! Section 3.4.2 Density of wave period, wave length or encountered wave period +#!#! Example 7: Crest period and high crest waves clf tic xx = load('sea.dat'); @@ -354,7 +353,7 @@ t = linspace(0.01,8,200); ftc1 = kde(Tc,{'L2',0},t); pdfplot(ftc1) 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); simpson(f_t.x{1},f_t.f) pdfplot(f_t,'-.') @@ -363,7 +362,7 @@ wafostamp([],'(ER)') toc disp('Block = 27'), pause(pstate) -%% +#!#! clf tic @@ -372,7 +371,7 @@ if strncmpi(speed,'slow',1) else NIT = 2; end -% f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4); +#! f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4); tic f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],Hs/2,NIT); toc @@ -390,14 +389,14 @@ wafostamp([],'(ER)') toc disp('Block = 28'), pause(pstate) -%% Example 8: Wave period for high crest waves -% clf +#!#! Example 8: Wave period for high crest waves +#! clf tic f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1); toc 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],1,5); +#! f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5); simpson(f_tcc3.x{1},f_tcc3.f) pdfplot(f_tcc2,'-.') hold on @@ -406,7 +405,7 @@ toc toc disp('Block = 29'), pause(pstate) -%% +#!#! clf [TC tc_ind v_ind] = dat2tc(yn,[],'dw'); N = length(tc_ind); @@ -418,7 +417,7 @@ spwaveplot(yn,ind(2:4)) wafostamp([],'(ER)') disp('Block = 30'), pause(pstate) -%% +#!#! clf Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1); 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('Some calculations may take some time.') disp('You could experiment with other NIT.') -%return +#!return -%% Section 3.4.3 Joint density of crest period and crest height -%% Example 9. Some preliminary analysis of the data +#!#! Section 3.4.3 Joint density of crest period and crest height +#!#! Example 9. Some preliminary analysis of the data clf tic yy = load('gfaksr89.dat'); @@ -457,7 +456,7 @@ v = v(2) toc disp('Block = 33'), pause(pstate) -%% +#!#! clf tic [TC, tc_ind, v_ind] = dat2tc(yy,v,'dw'); @@ -477,7 +476,7 @@ At = v-yy(t_ind,2); toc disp('Block = 34'), pause(pstate) -%% +#!#! clf tic t = linspace(0.01,15,200); @@ -497,8 +496,8 @@ wafostamp([],'(ER)') toc disp('Block = 35'), pause(pstate) -%% Example 10: Joint characteristics of a half wave: -%% position and height of a crest for a wave with given period +#!#! Example 10: Joint characteristics of a half wave: +#!#! position and height of a crest for a wave with given period clf tic ind = find(4.4v & Mm(:,2)>> import wafo.data >>> import wafo.objects as wo @@ -80,29 +83,30 @@ class LevelCrossings(WafoData): def __init__(self, *args, **kwds): options = dict(title='Level crossing spectrum', xlab='Levels', ylab='Count', + plotmethod='semilogy', plot_args=['b'], plot_args_children=['r--'],) options.update(**kwds) 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.setplotter(plotmethod='step') + #self.setplotter(plotmethod='step') icmax = self.data.argmax() 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)) logcmin = logcros[icmax] logcros = sqrt(2 * abs(logcros - logcmin)) logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1] p = polyfit(self.args[10:-9], logcros[10:-9], 1) #least square fit - if self.stdev is None: - self.stdev = 1.0 / p[0] #estimated standard deviation of x + if self.sigma is None: + self.sigma = 1.0 / p[0] #estimated standard deviation of x if self.mean is None: self.mean = -p[1] / p[0] #self.args[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) self.children = [WafoData(y, self.args)] @@ -212,7 +216,7 @@ class LevelCrossings(WafoData): scr = trapz(lc1 ** 2 * lc, lc1) scr = sqrt(scr - mcr ** 2) - lc2 = LevelCrossings(lc, lc1, mean=mcr, stdev=scr) + lc2 = LevelCrossings(lc, lc1, mean=mcr, sigma=scr) g = lc2.trdata() @@ -350,7 +354,7 @@ class LevelCrossings(WafoData): if mean is None: mean = self.mean if sigma is None: - sigma = self.stdev + sigma = self.sigma 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) @@ -401,7 +405,7 @@ class LevelCrossings(WafoData): 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) - 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: for ix in range(5): @@ -433,7 +437,7 @@ class CyclePairs(WafoData): data : array_like args : vector for 1D - Examples + Examples -------- >>> import wafo.data >>> import wafo.objects as wo @@ -446,7 +450,7 @@ class CyclePairs(WafoData): ''' def __init__(self, *args, **kwds): self.kind = kwds.pop('kind', 'max2min') - self.stdev = kwds.pop('stdev', None) + self.sigma = kwds.pop('sigma', None) self.mean = kwds.pop('mean', None) options = dict(title=self.kind + ' cycle pairs', @@ -590,7 +594,7 @@ class CyclePairs(WafoData): dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx] elif defnr == 3: ## This are upcrossings + minima + maxima 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): ''' @@ -613,7 +617,7 @@ class TurningPoints(WafoData): ''' def __init__(self, *args, **kwds): 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) options = dict(title='Turning points') @@ -672,7 +676,7 @@ class TurningPoints(WafoData): kind = 'max2min' 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): """ @@ -772,7 +776,7 @@ class TimeSeries(WafoData): with attributes: data : ACF vector 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: if R(k)=0 for all lags k>q then an approximation of the variance for large samples due to Bartlett @@ -828,14 +832,15 @@ class TimeSeries(WafoData): t = linspace(0, lag * dt, lag + 1) #cumsum = np.cumsum acf = _wafocov.CovData1D(R[lags], t) - acf.stdev = 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.sigma = sqrt(r_[ 0, r0**2 , r0**2 + 2 * cumsum(R[1:] ** 2)] / Ncens) + acf.children = [WafoData(-2. * acf.sigma[lags], t), WafoData(2. * acf.sigma[lags], t)] acf.plot_args_children = ['r:'] acf.norm = norm 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. Parameters @@ -963,7 +968,7 @@ class TimeSeries(WafoData): R = tsy.tocovdata() if estimate_L: #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 L = min(int(4 * L / 3), n - 2) print('The default L is set to %d' % L) @@ -1025,6 +1030,157 @@ class TimeSeries(WafoData): # S.S = zeros(nf+1,m-1); 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= 0): - vscale = [0, 2 * stdev * vfact] + vscale = [0, 2 * sigma * vfact] else: - vscale = array([-1, 1]) * vfact * stdev + vscale = array([-1, 1]) * vfact * sigma XlblTxt = 'Time [sec]' @@ -1580,7 +1736,7 @@ class TimeSeries(WafoData): dT = 1 / 60 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.)' plot = plotbackend.plot @@ -1602,7 +1758,7 @@ class TimeSeries(WafoData): #plotbackend.axis([h_scale*dT, v_scale]) for iy in [-2, 2]: - plot(h_scale * dT, iy * stdev * ones(2), ':') + plot(h_scale * dT, iy * sigma * ones(2), ':') ind = ind + Ns #end @@ -1630,7 +1786,6 @@ class TimeSeries(WafoData): >>> ts = wafo.objects.mat2timeseries(x[0:500,...]) >>> h = ts.plot_sp_wave(np.r_[6:9,12:18]) - See also -------- plot_wave, findtc diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index 7b11348..a4ce0d9 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -2088,9 +2088,9 @@ class SpecData1D(WafoData): 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() #R = spec2cov(S); diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py index b456d41..3c77b17 100644 --- a/pywafo/src/wafo/wafodata.py +++ b/pywafo/src/wafo/wafodata.py @@ -76,7 +76,7 @@ class WafoData(object): self.labels = AxisLabels(**kwds) if not self.plotter: - self.setplotter() + self.setplotter(kwds.get('plotmethod', None)) def plot(self, *args, **kwds): tmp = None