Updated tutorial ipython notebook scripts

master
Per.Andreas.Brodtkorb 11 years ago
parent a4d26ef265
commit d3080556a1

@ -11,34 +11,50 @@
"cell_type": "heading", "cell_type": "heading",
"level": 1, "level": 1,
"metadata": {}, "metadata": {},
"source": "CHAPTER 1 demonstrates some applications of WAFO" "source": [
"CHAPTER 1 demonstrates some applications of WAFO"
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "CHAPTER1 gives an overview through examples some of the capabilities of WAFO. WAFO is a toolbox of Matlab routines for statistical analysis and simulation of random waves and loads. The commands are edited for fast computation.\n" "source": [
"CHAPTER1 gives an overview through examples some of the capabilities of WAFO. WAFO is a toolbox of Matlab routines for statistical analysis and simulation of random waves and loads. The commands are edited for fast computation.\n"
]
}, },
{ {
"cell_type": "heading", "cell_type": "heading",
"level": 2, "level": 2,
"metadata": {}, "metadata": {},
"source": "Section 1.4 Some applications of WAFO" "source": [
"Section 1.4 Some applications of WAFO"
]
}, },
{ {
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Section 1.4.1 Simulation from spectrum, estimation of spectrum " "source": [
"Section 1.4.1 Simulation from spectrum, estimation of spectrum "
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Simulation of the sea surface from spectrum. The following code generates 200 seconds of data sampled with 10Hz from the Torsethaugen spectrum." "source": [
"Simulation of the sea surface from spectrum. The following code generates 200 seconds of data sampled with 10Hz from the Torsethaugen spectrum."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo.spectrum.models as wsm\nS = wsm.Torsethaugen(Hm0=6, Tp=8);\nS1 = S.tospecdata()\nS1.plot()\nshow()", "input": [
"import wafo.spectrum.models as wsm\n",
"S = wsm.Torsethaugen(Hm0=6, Tp=8);\n",
"S1 = S.tospecdata()\n",
"S1.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -52,7 +68,13 @@
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo.objects as wo\nxs = S1.sim(ns=2000, dt=0.1)\nts = wo.mat2timeseries(xs)\nts.plot_wave('-')\nshow()", "input": [
"import wafo.objects as wo\n",
"xs = S1.sim(ns=2000, dt=0.1)\n",
"ts = wo.mat2timeseries(xs)\n",
"ts.plot_wave('-')\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -67,17 +89,31 @@
"cell_type": "heading", "cell_type": "heading",
"level": 4, "level": 4,
"metadata": {}, "metadata": {},
"source": "Estimation of spectrum " "source": [
"Estimation of spectrum "
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "A common situation is that one wants to estimate the spectrum for wave measurements. The following code simulate 20 minutes signal sampled at 4Hz and compare the spectral estimate with the original Torsethaugen spectum.\n" "source": [
"A common situation is that one wants to estimate the spectrum for wave measurements. The following code simulate 20 minutes signal sampled at 4Hz and compare the spectral estimate with the original Torsethaugen spectum.\n"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nFs = 4; \nxs = S1.sim(ns=fix(20 * 60 * Fs), dt=1. / Fs) \nts = wo.mat2timeseries(xs) \nSest = ts.tospecdata(L=400)\nS1.plot()\nSest.plot('--')\naxis([0, 3, 0, 5]) # This may depend on the simulation\nshow()", "input": [
"clf()\n",
"Fs = 4; \n",
"xs = S1.sim(ns=fix(20 * 60 * Fs), dt=1. / Fs) \n",
"ts = wo.mat2timeseries(xs) \n",
"Sest = ts.tospecdata(L=400)\n",
"S1.plot()\n",
"Sest.plot('--')\n",
"axis([0, 3, 0, 5]) # This may depend on the simulation\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -92,17 +128,35 @@
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Section 1.4.2 Probability distributions of wave characteristics." "source": [
"Section 1.4.2 Probability distributions of wave characteristics."
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Probability distribution of wave trough period: WAFO gives the possibility of computing the exact probability distributions for a number of characteristics given a spectral density. In the following example we study the trough period extracted from the time series and compared with the theoretical density computed with exact spectrum, S1, and the estimated spectrum, Sest.\n" "source": [
"Probability distribution of wave trough period: WAFO gives the possibility of computing the exact probability distributions for a number of characteristics given a spectral density. In the following example we study the trough period extracted from the time series and compared with the theoretical density computed with exact spectrum, S1, and the estimated spectrum, Sest.\n"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nimport wafo.misc as wm\ndtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)\ndtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)\n\nT, index = ts.wave_periods(vh=0, pdef='d2u')\nbins = wm.good_bins(T, num_bins=25, odd=True)\nwm.plot_histgrm(T, bins=bins, normed=True)\n\ndtyex.plot()\ndtyest.plot('-.')\naxis([0, 10, 0, 0.35])\nshow()", "input": [
"clf()\n",
"import wafo.misc as wm\n",
"dtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)\n",
"dtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)\n",
"\n",
"T, index = ts.wave_periods(vh=0, pdef='d2u')\n",
"bins = wm.good_bins(T, num_bins=25, odd=True)\n",
"wm.plot_histgrm(T, bins=bins, normed=True)\n",
"\n",
"dtyex.plot()\n",
"dtyest.plot('-.')\n",
"axis([0, 10, 0, 0.35])\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -117,17 +171,36 @@
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Section 1.4.3 Directional spectra" "source": [
"Section 1.4.3 Directional spectra"
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Here are a few lines of code, which produce directional spectra with frequency independent and frequency dependent spreading." "source": [
"Here are a few lines of code, which produce directional spectra with frequency independent and frequency dependent spreading."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nplotflag = 1\nNt = 101; # number of angles\nth0 = pi / 2; # primary direction of waves\nSp = 15; # spreading parameter\n\nD1 = wsm.Spreading(type='cos', theta0=th0, method=None) # frequency independent\nD12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent\n\nSD1 = D1.tospecdata2d(S1)\nSD12 = D12.tospecdata2d(S1)\nSD1.plot()\nSD12.plot()#linestyle='dashdot')\nshow()", "input": [
"clf()\n",
"plotflag = 1\n",
"Nt = 101; # number of angles\n",
"th0 = pi / 2; # primary direction of waves\n",
"Sp = 15; # spreading parameter\n",
"\n",
"D1 = wsm.Spreading(type='cos', theta0=th0, method=None) # frequency independent\n",
"D12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent\n",
"\n",
"SD1 = D1.tospecdata2d(S1)\n",
"SD12 = D12.tospecdata2d(S1)\n",
"SD1.plot()\n",
"SD12.plot()#linestyle='dashdot')\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -142,23 +215,38 @@
"cell_type": "heading", "cell_type": "heading",
"level": 4, "level": 4,
"metadata": {}, "metadata": {},
"source": "3D Simulation of the sea surface " "source": [
"3D Simulation of the sea surface "
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "The simulations show that frequency dependent spreading leads to much more irregular surface so the orientation of waves is less transparent compared to the frequency independent case." "source": [
"The simulations show that frequency dependent spreading leads to much more irregular surface so the orientation of waves is less transparent compared to the frequency independent case."
]
}, },
{ {
"cell_type": "heading", "cell_type": "heading",
"level": 5, "level": 5,
"metadata": {}, "metadata": {},
"source": "Frequency independent spreading" "source": [
"Frequency independent spreading"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "#plotflag = 1; iseed = 1;\n#\n#Nx = 2 ^ 8;Ny = Nx;Nt = 1;dx = 0.5; dy = dx; dt = 0.25; fftdim = 2;\n#randn('state', iseed)\n#Y1 = seasim(SD1, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);\n#wafostamp('', '(ER)')\n#axis('fill')\n#disp('Block = 6'), pause(pstate)", "input": [
"#plotflag = 1; iseed = 1;\n",
"#\n",
"#Nx = 2 ^ 8;Ny = Nx;Nt = 1;dx = 0.5; dy = dx; dt = 0.25; fftdim = 2;\n",
"#randn('state', iseed)\n",
"#Y1 = seasim(SD1, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);\n",
"#wafostamp('', '(ER)')\n",
"#axis('fill')\n",
"#disp('Block = 6'), pause(pstate)"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []
@ -167,12 +255,19 @@
"cell_type": "heading", "cell_type": "heading",
"level": 5, "level": 5,
"metadata": {}, "metadata": {},
"source": "Frequency dependent spreading" "source": [
"Frequency dependent spreading"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "#randn('state', iseed)\n#Y12 = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);\n#wafostamp('', '(ER)')\n#axis('fill')", "input": [
"#randn('state', iseed)\n",
"#Y12 = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);\n",
"#wafostamp('', '(ER)')\n",
"#axis('fill')"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []
@ -181,17 +276,37 @@
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Estimation of directional spectrum" "source": [
"Estimation of directional spectrum"
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "The figure is not shown in the Tutorial" "source": [
"The figure is not shown in the Tutorial"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "# Nx = 3; Ny = 2; Nt = 2 ^ 12; dx = 10; dy = 10;dt = 0.5;\n# F = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, 1, 0); \n# Z = permute(F.Z, [3 1 2]);\n# [X, Y] = meshgrid(F.x, F.y);\n# N = Nx * Ny;\n# types = repmat(sensortypeid('n'), N, 1);\n# bfs = ones(N, 1);\n# pos = [X(:), Y(:), zeros(N, 1)];\n# h = inf;\n# nfft = 128;\n# nt = 101;\n# SDe = dat2dspec([F.t Z(:, :)], [pos types, bfs], h, nfft, nt);\n#plotspec(SDe), hold on\n#plotspec(SD12, '--'), hold off\n#disp('Block = 8'), pause(pstate)\n", "input": [
"# Nx = 3; Ny = 2; Nt = 2 ^ 12; dx = 10; dy = 10;dt = 0.5;\n",
"# F = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, 1, 0); \n",
"# Z = permute(F.Z, [3 1 2]);\n",
"# [X, Y] = meshgrid(F.x, F.y);\n",
"# N = Nx * Ny;\n",
"# types = repmat(sensortypeid('n'), N, 1);\n",
"# bfs = ones(N, 1);\n",
"# pos = [X(:), Y(:), zeros(N, 1)];\n",
"# h = inf;\n",
"# nfft = 128;\n",
"# nt = 101;\n",
"# SDe = dat2dspec([F.t Z(:, :)], [pos types, bfs], h, nfft, nt);\n",
"#plotspec(SDe), hold on\n",
"#plotspec(SD12, '--'), hold off\n",
"#disp('Block = 8'), pause(pstate)\n"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []
@ -200,17 +315,37 @@
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Section 1.4.4 Fatigue, Load cycles and Markov models" "source": [
"Section 1.4.4 Fatigue, Load cycles and Markov models"
]
}, },
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Switching Markow chain of turningpoints.\nIn fatigue applications the exact sample path is not important, but only the tops and bottoms of the load, called the sequence of turning points (TP). From the turning points one can extract load cycles, from which damage calculations and fatigue life predictions can be performed.\n\nThe commands below computes the intensity of rainflowcycles for the Gaussian model with spectrum S1 using the Markov approximation. \nThe rainflow cycles found in the simulated load signal are shown in the figure." "source": [
"Switching Markow chain of turningpoints.\n",
"In fatigue applications the exact sample path is not important, but only the tops and bottoms of the load, called the sequence of turning points (TP). From the turning points one can extract load cycles, from which damage calculations and fatigue life predictions can be performed.\n",
"\n",
"The commands below computes the intensity of rainflowcycles for the Gaussian model with spectrum S1 using the Markov approximation. \n",
"The rainflow cycles found in the simulated load signal are shown in the figure."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "#clf()\n#paramu = [-6 6 61];\n#frfc = spec2cmat(S1, [], 'rfc', [], paramu);\n#pdfplot(frfc);\n#hold on\n#tp = dat2tp(xs);\n#rfc = tp2rfc(tp);\n#plot(rfc(:, 2), rfc(:, 1), '.')\n#wafostamp('', '(ER)')\n#hold off\n#disp('Block = 9'), pause(pstate)", "input": [
"#clf()\n",
"#paramu = [-6 6 61];\n",
"#frfc = spec2cmat(S1, [], 'rfc', [], paramu);\n",
"#pdfplot(frfc);\n",
"#hold on\n",
"#tp = dat2tp(xs);\n",
"#rfc = tp2rfc(tp);\n",
"#plot(rfc(:, 2), rfc(:, 1), '.')\n",
"#wafostamp('', '(ER)')\n",
"#hold off\n",
"#disp('Block = 9'), pause(pstate)"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []
@ -219,19 +354,32 @@
"cell_type": "heading", "cell_type": "heading",
"level": 3, "level": 3,
"metadata": {}, "metadata": {},
"source": "Section 1.4.5 Extreme value statistics" "source": [
"Section 1.4.5 Extreme value statistics"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nimport wafo.data as wd\nxn = wd.yura87()\n#xn = load('yura87.dat'); \nsubplot(211) \nplot(xn[::30, 0] / 3600, xn[::30, 1], '.')\ntitle('Water level')\nylabel('(m)')", "input": [
"clf()\n",
"import wafo.data as wd\n",
"xn = wd.yura87()\n",
"#xn = load('yura87.dat'); \n",
"subplot(211) \n",
"plot(xn[::30, 0] / 3600, xn[::30, 1], '.')\n",
"title('Water level')\n",
"ylabel('(m)')"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "pyout", "output_type": "pyout",
"prompt_number": 10, "prompt_number": 10,
"text": "<matplotlib.text.Text at 0x730d070>" "text": [
"<matplotlib.text.Text at 0x730d070>"
]
}, },
{ {
"output_type": "display_data", "output_type": "display_data",
@ -243,12 +391,24 @@
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Formation of 5 min maxima" "source": [
"Formation of 5 min maxima"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "yura = xn[:85500, 1]\nyura = np.reshape(yura, (285, 300)).T\nmaxyura = yura.max(axis=0)\nsubplot(212)\nplot(xn[299:85500:300, 0] / 3600, maxyura, '.')\nxlabel('Time (h)')\nylabel('(m)')\ntitle('Maximum 5 min water level')\nshow()", "input": [
"yura = xn[:85500, 1]\n",
"yura = np.reshape(yura, (285, 300)).T\n",
"maxyura = yura.max(axis=0)\n",
"subplot(212)\n",
"plot(xn[299:85500:300, 0] / 3600, maxyura, '.')\n",
"xlabel('Time (h)')\n",
"ylabel('(m)')\n",
"title('Maximum 5 min water level')\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -262,19 +422,31 @@
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Estimation of GEV for yuramax" "source": [
"Estimation of GEV for yuramax"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nimport wafo.stats as ws\nphat = ws.genextreme.fit2(maxyura, method='ml')\nphat.plotfitsummary()\nshow()\n#disp('Block = 11, Last block')", "input": [
"clf()\n",
"import wafo.stats as ws\n",
"phat = ws.genextreme.fit2(maxyura, method='ml')\n",
"phat.plotfitsummary()\n",
"show()\n",
"#disp('Block = 11, Last block')"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stderr", "stream": "stderr",
"text": "c:\\pab\\workspace\\pywafo_svn\\pywafo\\src\\wafo\\stats\\estimation.py:1080: UserWarning: P-value is on the conservative side (i.e. too large) due to ties in the data!\n warnings.warn('P-value is on the conservative side (i.e. too large) due to ties in the data!')\n" "text": [
"c:\\pab\\workspace\\pywafo_svn\\pywafo\\src\\wafo\\stats\\estimation.py:1080: UserWarning: P-value is on the conservative side (i.e. too large) due to ties in the data!\n",
" warnings.warn('P-value is on the conservative side (i.e. too large) due to ties in the data!')\n"
]
}, },
{ {
"output_type": "display_data", "output_type": "display_data",
@ -282,14 +454,6 @@
} }
], ],
"prompt_number": 12 "prompt_number": 12
},
{
"cell_type": "code",
"collapsed": false,
"input": "",
"language": "python",
"metadata": {},
"outputs": []
} }
], ],
"metadata": {} "metadata": {}

@ -10,12 +10,44 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "CHAPTER2 Modelling random loads and stochastic waves\n====================================================\n\nChapter2 contains the commands used in Chapter 2 of the tutorial and present some tools for analysis of random functions with respect to their correlation, spectral and distributional properties. The presentation is divided into three examples: \n\nExample1 is devoted to estimation of different parameters in the model.\nExample2 deals with spectral densities and\nExample3 presents the use of WAFO to simulate samples of a Gaussian process.\n\nSome of the commands are edited for fast computation. \n\nSection 2.1 Introduction and preliminary analysis\n=================================================\n\nExample 1: Sea data\n-------------------\nObserved crossings compared to the expected for Gaussian signals\n" "source": [
"CHAPTER2 Modelling random loads and stochastic waves\n",
"====================================================\n",
"\n",
"Chapter2 contains the commands used in Chapter 2 of the tutorial and present some tools for analysis of random functions with respect to their correlation, spectral and distributional properties. The presentation is divided into three examples: \n",
"\n",
"Example1 is devoted to estimation of different parameters in the model.\n",
"Example2 deals with spectral densities and\n",
"Example3 presents the use of WAFO to simulate samples of a Gaussian process.\n",
"\n",
"Some of the commands are edited for fast computation. \n",
"\n",
"Section 2.1 Introduction and preliminary analysis\n",
"=================================================\n",
"\n",
"Example 1: Sea data\n",
"-------------------\n",
"Observed crossings compared to the expected for Gaussian signals\n"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo\nimport wafo.objects as wo\nxx = wafo.data.sea()\nme = xx[:, 1].mean()\nsa = xx[:, 1].std()\nxx[:, 1] -= me\nts = wo.mat2timeseries(xx)\ntp = ts.turning_points()\n\ncc = tp.cycle_pairs()\nlc = cc.level_crossings()\nlc.plot()\nshow()", "input": [
"import wafo\n",
"import wafo.objects as wo\n",
"xx = wafo.data.sea()\n",
"me = xx[:, 1].mean()\n",
"sa = xx[:, 1].std()\n",
"xx[:, 1] -= me\n",
"ts = wo.mat2timeseries(xx)\n",
"tp = ts.turning_points()\n",
"\n",
"cc = tp.cycle_pairs()\n",
"lc = cc.level_crossings()\n",
"lc.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -29,19 +61,29 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Average number of upcrossings per time unit\n----------------------------------------------\nNext 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. \n" "source": [
"Average number of upcrossings per time unit\n",
"----------------------------------------------\n",
"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. \n"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "T = xx[:, 0].max() - xx[:, 0].min()\nf0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency \nprint('f0 = %g' % f0)", "input": [
"T = xx[:, 0].max() - xx[:, 0].min()\n",
"f0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency \n",
"print('f0 = %g' % f0)"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stdout", "stream": "stdout",
"text": "f0 = 0.224071\n" "text": [
"f0 = 0.224071\n"
]
} }
], ],
"prompt_number": 3 "prompt_number": 3
@ -49,19 +91,29 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Turningpoints and irregularity factor\n----------------------------------------" "source": [
"Turningpoints and irregularity factor\n",
"----------------------------------------"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "fm = len(tp.data) / (2 * T) # frequency of maxima\nalfa = f0 / fm # approx Tm24/Tm02\n\nprint('fm = %g, alpha = %g, ' % (fm, alfa))", "input": [
"fm = len(tp.data) / (2 * T) # frequency of maxima\n",
"alfa = f0 / fm # approx Tm24/Tm02\n",
"\n",
"print('fm = %g, alpha = %g, ' % (fm, alfa))"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stdout", "stream": "stdout",
"text": "fm = 0.456159, alpha = 0.491212, \n" "text": [
"fm = 0.456159, alpha = 0.491212, \n"
]
} }
], ],
"prompt_number": 4 "prompt_number": 4
@ -69,12 +121,23 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Visually examine data\n------------------------\nWe finish this section with some remarks about the quality of the measured data. Especially sea surface measurements can be of poor quality. We shall now check the quality of the dataset {\\tt xx}. It is always good practice to visually examine the data before the analysis to get an impression of the quality, \nnon-linearities and narrow-bandedness of the data.First we shall plot the data and zoom in on a specific region. A part of sea data is visualized with the following commands" "source": [
"Visually examine data\n",
"------------------------\n",
"We finish this section with some remarks about the quality of the measured data. Especially sea surface measurements can be of poor quality. We shall now check the quality of the dataset {\\tt xx}. It is always good practice to visually examine the data before the analysis to get an impression of the quality, \n",
"non-linearities and narrow-bandedness of the data.First we shall plot the data and zoom in on a specific region. A part of sea data is visualized with the following commands"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nts.plot_wave('k-', tp, '*', nfig=1, nsub=1)\n\naxis([0, 2, -2, 2])\nshow()", "input": [
"clf()\n",
"ts.plot_wave('k-', tp, '*', nfig=1, nsub=1)\n",
"\n",
"axis([0, 2, -2, 2])\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -88,19 +151,39 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Finding possible spurious points\n------------------------------------\nHowever, if the amount of data is too large for visual examinations one could use the following criteria to find possible spurious points. One must be careful using the criteria for extremevalue analysis, because\nit might remove extreme waves that are OK and not spurious." "source": [
"Finding possible spurious points\n",
"------------------------------------\n",
"However, if the amount of data is too large for visual examinations one could use the following criteria to find possible spurious points. One must be careful using the criteria for extremevalue analysis, because\n",
"it might remove extreme waves that are OK and not spurious."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo.misc as wm\ndt = ts.sampling_period()\n# dt = np.diff(xx[:2,0])\ndcrit = 5 * dt\nddcrit = 9.81 / 2 * dt * dt\nzcrit = 0\ninds, indg = wm.findoutliers(ts.data, zcrit, dcrit, ddcrit, verbose=True)", "input": [
"import wafo.misc as wm\n",
"dt = ts.sampling_period()\n",
"# dt = np.diff(xx[:2,0])\n",
"dcrit = 5 * dt\n",
"ddcrit = 9.81 / 2 * dt * dt\n",
"zcrit = 0\n",
"inds, indg = wm.findoutliers(ts.data, zcrit, dcrit, ddcrit, verbose=True)"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stdout", "stream": "stdout",
"text": "Found 0 spurious positive jumps of Dx\nFound 0 spurious negative jumps of Dx\nFound 37 spurious positive jumps of D^2x\nFound 200 spurious negative jumps of D^2x\nFound 244 consecutive equal values\nFound the total of 1152 spurious points\n" "text": [
"Found 0 spurious positive jumps of Dx\n",
"Found 0 spurious negative jumps of Dx\n",
"Found 37 spurious positive jumps of D^2x\n",
"Found 200 spurious negative jumps of D^2x\n",
"Found 244 consecutive equal values\n",
"Found the total of 1152 spurious points\n"
]
} }
], ],
"prompt_number": 6 "prompt_number": 6
@ -108,12 +191,23 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Section 2.2 Frequency Modeling of Load Histories\n---------------------------------------------------\nPeriodogram: Raw spectrum" "source": [
"Section 2.2 Frequency Modeling of Load Histories\n",
"---------------------------------------------------\n",
"Periodogram: Raw spectrum"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nLmax = 9500\nS = ts.tospecdata(L=Lmax)\nS.plot()\naxis([0, 5, 0, 0.7])\nshow()", "input": [
"clf()\n",
"Lmax = 9500\n",
"S = ts.tospecdata(L=Lmax)\n",
"S.plot()\n",
"axis([0, 5, 0, 0.7])\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -127,19 +221,27 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Calculate moments \n-------------------" "source": [
"Calculate moments \n",
"-------------------"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "mom, text = S.moment(nr=4)\nprint('sigma = %g, m0 = %g' % (sa, sqrt(mom[0])))", "input": [
"mom, text = S.moment(nr=4)\n",
"print('sigma = %g, m0 = %g' % (sa, sqrt(mom[0])))"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stdout", "stream": "stdout",
"text": "sigma = 0.472955, m0 = 0.472955\n" "text": [
"sigma = 0.472955, m0 = 0.472955\n"
]
} }
], ],
"prompt_number": 8 "prompt_number": 8
@ -147,12 +249,26 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Section 2.2.1 Random functions in Spectral Domain - Gaussian processes\n--------------------------------------------------------------------------\nSmoothing of spectral estimate \n----------------------------------\nBy decreasing Lmax the spectrum estimate becomes smoother." "source": [
"Section 2.2.1 Random functions in Spectral Domain - Gaussian processes\n",
"--------------------------------------------------------------------------\n",
"Smoothing of spectral estimate \n",
"----------------------------------\n",
"By decreasing Lmax the spectrum estimate becomes smoother."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nLmax0 = 200; Lmax1 = 50\nS1 = ts.tospecdata(L=Lmax0)\nS2 = ts.tospecdata(L=Lmax1)\nS1.plot('-.')\nS2.plot()\nshow()", "input": [
"clf()\n",
"Lmax0 = 200; Lmax1 = 50\n",
"S1 = ts.tospecdata(L=Lmax0)\n",
"S2 = ts.tospecdata(L=Lmax1)\n",
"S1.plot('-.')\n",
"S2.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -166,12 +282,28 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": " Estimated autocovariance\n----------------------------\nObviously knowing the spectrum one can compute the covariance\nfunction. The following code will compute the covariance for the \nunimodal spectral density S1 and compare it with estimated \ncovariance of the signal xx." "source": [
" Estimated autocovariance\n",
"----------------------------\n",
"Obviously knowing the spectrum one can compute the covariance\n",
"function. The following code will compute the covariance for the \n",
"unimodal spectral density S1 and compare it with estimated \n",
"covariance of the signal xx."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nLmax = 85\nR1 = S1.tocovdata(nr=1) \nRest = ts.tocovdata(lag=Lmax)\nR1.plot('.')\nRest.plot()\naxis([0, 25, -0.1, 0.25])\nshow()", "input": [
"clf()\n",
"Lmax = 85\n",
"R1 = S1.tocovdata(nr=1) \n",
"Rest = ts.tocovdata(lag=Lmax)\n",
"R1.plot('.')\n",
"Rest.plot()\n",
"axis([0, 25, -0.1, 0.25])\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -185,12 +317,20 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "We can see in Figure below that the covariance function corresponding to the spectral density S2 significantly differs from the one estimated directly from data. It can be seen in Figure above that the covariance corresponding to S1 agrees much better with the estimated covariance function." "source": [
"We can see in Figure below that the covariance function corresponding to the spectral density S2 significantly differs from the one estimated directly from data. It can be seen in Figure above that the covariance corresponding to S1 agrees much better with the estimated covariance function."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nR2 = S2.tocovdata(nr=1)\nR2.plot('.')\nRest.plot()\nshow()", "input": [
"clf()\n",
"R2 = S2.tocovdata(nr=1)\n",
"R2.plot('.')\n",
"Rest.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -204,12 +344,22 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Section 2.2.2 Transformed Gaussian models\n-------------------------------------------\nWe begin with computing skewness and kurtosis for the data set xx and compare it with the second order wave approximation proposed by Winterstein:" "source": [
"Section 2.2.2 Transformed Gaussian models\n",
"-------------------------------------------\n",
"We begin with computing skewness and kurtosis for the data set xx and compare it with the second order wave approximation proposed by Winterstein:"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo.stats as ws\nrho3 = ws.skew(xx[:, 1])\nrho4 = ws.kurtosis(xx[:, 1])\n\nsk, ku = S1.stats_nl(moments='sk')", "input": [
"import wafo.stats as ws\n",
"rho3 = ws.skew(xx[:, 1])\n",
"rho4 = ws.kurtosis(xx[:, 1])\n",
"\n",
"sk, ku = S1.stats_nl(moments='sk')"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
@ -218,12 +368,26 @@
{ {
"cell_type": "raw", "cell_type": "raw",
"metadata": {}, "metadata": {},
"source": "Comparisons of 3 transformations" "source": [
"Comparisons of 3 transformations"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nimport wafo.transform.models as wtm\ngh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata()\ng = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation \nglc, gemp = lc.trdata(mean=me, sigma=sa)\n\nglc.plot('b-') #! Transf. estimated from level-crossings\ngh.plot('b-.') #! Hermite Transf. estimated from moments\ng.plot('r')\ngrid('on')\nshow()", "input": [
"clf()\n",
"import wafo.transform.models as wtm\n",
"gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata()\n",
"g = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation \n",
"glc, gemp = lc.trdata(mean=me, sigma=sa)\n",
"\n",
"glc.plot('b-') #! Transf. estimated from level-crossings\n",
"gh.plot('b-.') #! Hermite Transf. estimated from moments\n",
"g.plot('r')\n",
"grid('on')\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -237,19 +401,37 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Test Gaussianity of a stochastic process\n------------------------------------------\nTESTGAUSSIAN 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.\nIf 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.\n\nAs 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. " "source": [
"Test Gaussianity of a stochastic process\n",
"------------------------------------------\n",
"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.\n",
"If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.\n",
"\n",
"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. "
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\ntest0 = glc.dist2gauss()\n# the following test takes time\nN = len(xx)\ntest1 = S1.testgaussian(ns=N, cases=50, test0=test0)\nis_gaussian = sum(test1 > test0) > 5 \nprint(is_gaussian)\nshow()", "input": [
"clf()\n",
"test0 = glc.dist2gauss()\n",
"# the following test takes time\n",
"N = len(xx)\n",
"test1 = S1.testgaussian(ns=N, cases=50, test0=test0)\n",
"is_gaussian = sum(test1 > test0) > 5 \n",
"print(is_gaussian)\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"output_type": "stream", "output_type": "stream",
"stream": "stdout", "stream": "stdout",
"text": "False\n" "text": [
"False\n"
]
}, },
{ {
"output_type": "display_data", "output_type": "display_data",
@ -261,12 +443,21 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Normalplot of data xx\n------------------------\nindicates that the underlying distribution has a \"heavy\" upper tail and a \"light\" lower tail." "source": [
"Normalplot of data xx\n",
"------------------------\n",
"indicates that the underlying distribution has a \"heavy\" upper tail and a \"light\" lower tail."
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nimport pylab\nws.probplot(ts.data.ravel(), dist='norm', plot=pylab)\nshow()", "input": [
"clf()\n",
"import pylab\n",
"ws.probplot(ts.data.ravel(), dist='norm', plot=pylab)\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -280,12 +471,23 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Section 2.2.3 Spectral densities of sea data\n-----------------------------------------------\nExample 2: Different forms of spectra" "source": [
"Section 2.2.3 Spectral densities of sea data\n",
"-----------------------------------------------\n",
"Example 2: Different forms of spectra"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "import wafo.spectrum.models as wsm\nclf()\nHm0 = 7; Tp = 11;\nspec = wsm.Jonswap(Hm0=Hm0, Tp=Tp).tospecdata()\nspec.plot()\nshow()", "input": [
"import wafo.spectrum.models as wsm\n",
"clf()\n",
"Hm0 = 7; Tp = 11;\n",
"spec = wsm.Jonswap(Hm0=Hm0, Tp=Tp).tospecdata()\n",
"spec.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -299,12 +501,23 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Directional spectrum and Encountered directional spectrum\n=========================================================\nDirectional spectrum\n---------------------" "source": [
"Directional spectrum and Encountered directional spectrum\n",
"=========================================================\n",
"Directional spectrum\n",
"---------------------"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nD = wsm.Spreading('cos2s')\nSd = D.tospecdata2d(spec)\nSd.plot()\nshow()", "input": [
"clf()\n",
"D = wsm.Spreading('cos2s')\n",
"Sd = D.tospecdata2d(spec)\n",
"Sd.plot()\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -318,12 +531,95 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Encountered directional spectrum\n--------------------------------- " "source": [
"Encountered directional spectrum\n",
"--------------------------------- "
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "#clf()\n#Se = spec2spec(Sd,'encdir',0,10);\n#plotspec(Se), hold on\n#plotspec(Sd,1,'--'), hold off\n##!wafostamp('','(ER)')\n#disp('Block = 17'),pause(pstate)\n#\n##!#! Frequency spectra\n#clf\n#Sd1 =spec2spec(Sd,'freq');\n#Sd2 = spec2spec(Se,'enc');\n#plotspec(spec), hold on\n#plotspec(Sd1,1,'.'),\n#plotspec(Sd2),\n##!wafostamp('','(ER)')\n#hold off\n#disp('Block = 18'),pause(pstate)\n#\n##!#! Wave number spectrum\n#clf\n#Sk = spec2spec(spec,'k1d')\n#Skd = spec2spec(Sd,'k1d')\n#plotspec(Sk), hold on\n#plotspec(Skd,1,'--'), hold off\n##!wafostamp('','(ER)')\n#disp('Block = 19'),pause(pstate)\n#\n##!#! Effect of waterdepth on spectrum\n#clf\n#plotspec(spec,1,'--'), hold on\n#S20 = spec;\n#S20.S = S20.S.*phi1(S20.w,20);\n#S20.h = 20;\n#plotspec(S20), hold off\n##!wafostamp('','(ER)')\n#disp('Block = 20'),pause(pstate)\n#\n##!#! Section 2.3 Simulation of transformed Gaussian process\n##!#! Example 3: Simulation of random sea \n##! The reconstruct function replaces the spurious points of seasurface by\n##! simulated data on the basis of the remaining data and a transformed Gaussian\n##! process. As noted previously one must be careful using the criteria \n##! for finding spurious points when reconstructing a dataset, because\n##! these criteria might remove the highest and steepest waves as we can see\n##! in this plot where the spurious points is indicated with a '+' sign:\n##!\n#clf\n#[y, grec] = reconstruct(xx,inds);\n#waveplot(y,'-',xx(inds,:),'+',1,1)\n#axis([0 inf -inf inf])\n##!wafostamp('','(ER)')\n#disp('Block = 21'),pause(pstate)\n#\n##! Compare transformation (grec) from reconstructed (y) \n##! with original (glc) from (xx)\n#clf\n#trplot(g), hold on\n#plot(gemp(:,1),gemp(:,2))\n#plot(glc(:,1),glc(:,2),'-.')\n#plot(grec(:,1),grec(:,2)), hold off \n#disp('Block = 22'),pause(pstate)\n#\n##!#!\n#clf\n#L = 200;\n#x = dat2gaus(y,grec);\n#Sx = dat2spec(x,L);\n#disp('Block = 23'),pause(pstate)\n# \n##!#!\n#clf\n#dt = spec2dt(Sx)\n#Ny = fix(2*60/dt) #! = 2 minutes\n#Sx.tr = grec;\n#ysim = spec2sdat(Sx,Ny);\n#waveplot(ysim,'-')\n##!wafostamp('','(CR)')\n#disp('Block = 24'),pause(pstate)\n\n", "input": [
"#clf()\n",
"#Se = spec2spec(Sd,'encdir',0,10);\n",
"#plotspec(Se), hold on\n",
"#plotspec(Sd,1,'--'), hold off\n",
"##!wafostamp('','(ER)')\n",
"#disp('Block = 17'),pause(pstate)\n",
"#\n",
"##!#! Frequency spectra\n",
"#clf\n",
"#Sd1 =spec2spec(Sd,'freq');\n",
"#Sd2 = spec2spec(Se,'enc');\n",
"#plotspec(spec), hold on\n",
"#plotspec(Sd1,1,'.'),\n",
"#plotspec(Sd2),\n",
"##!wafostamp('','(ER)')\n",
"#hold off\n",
"#disp('Block = 18'),pause(pstate)\n",
"#\n",
"##!#! Wave number spectrum\n",
"#clf\n",
"#Sk = spec2spec(spec,'k1d')\n",
"#Skd = spec2spec(Sd,'k1d')\n",
"#plotspec(Sk), hold on\n",
"#plotspec(Skd,1,'--'), hold off\n",
"##!wafostamp('','(ER)')\n",
"#disp('Block = 19'),pause(pstate)\n",
"#\n",
"##!#! Effect of waterdepth on spectrum\n",
"#clf\n",
"#plotspec(spec,1,'--'), hold on\n",
"#S20 = spec;\n",
"#S20.S = S20.S.*phi1(S20.w,20);\n",
"#S20.h = 20;\n",
"#plotspec(S20), hold off\n",
"##!wafostamp('','(ER)')\n",
"#disp('Block = 20'),pause(pstate)\n",
"#\n",
"##!#! Section 2.3 Simulation of transformed Gaussian process\n",
"##!#! Example 3: Simulation of random sea \n",
"##! The reconstruct function replaces the spurious points of seasurface by\n",
"##! simulated data on the basis of the remaining data and a transformed Gaussian\n",
"##! process. As noted previously one must be careful using the criteria \n",
"##! for finding spurious points when reconstructing a dataset, because\n",
"##! these criteria might remove the highest and steepest waves as we can see\n",
"##! in this plot where the spurious points is indicated with a '+' sign:\n",
"##!\n",
"#clf\n",
"#[y, grec] = reconstruct(xx,inds);\n",
"#waveplot(y,'-',xx(inds,:),'+',1,1)\n",
"#axis([0 inf -inf inf])\n",
"##!wafostamp('','(ER)')\n",
"#disp('Block = 21'),pause(pstate)\n",
"#\n",
"##! Compare transformation (grec) from reconstructed (y) \n",
"##! with original (glc) from (xx)\n",
"#clf\n",
"#trplot(g), hold on\n",
"#plot(gemp(:,1),gemp(:,2))\n",
"#plot(glc(:,1),glc(:,2),'-.')\n",
"#plot(grec(:,1),grec(:,2)), hold off \n",
"#disp('Block = 22'),pause(pstate)\n",
"#\n",
"##!#!\n",
"#clf\n",
"#L = 200;\n",
"#x = dat2gaus(y,grec);\n",
"#Sx = dat2spec(x,L);\n",
"#disp('Block = 23'),pause(pstate)\n",
"# \n",
"##!#!\n",
"#clf\n",
"#dt = spec2dt(Sx)\n",
"#Ny = fix(2*60/dt) #! = 2 minutes\n",
"#Sx.tr = grec;\n",
"#ysim = spec2sdat(Sx,Ny);\n",
"#waveplot(ysim,'-')\n",
"##!wafostamp('','(CR)')\n",
"#disp('Block = 24'),pause(pstate)\n",
"\n"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []
@ -331,12 +627,24 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Estimated spectrum compared to Torsethaugen spectrum\n-------------------------------------------------------" "source": [
"Estimated spectrum compared to Torsethaugen spectrum\n",
"-------------------------------------------------------"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "clf()\nfp = 1.1;dw = 0.01\nH0 = S1.characteristic('Hm0')[0]\nSt = wsm.Torsethaugen(Hm0=H0,Tp=2*pi/fp).tospecdata(np.arange(0,5+dw/2,dw)) \nS1.plot()\nSt.plot('-.')\naxis([0, 6, 0, 0.4])\nshow()\n", "input": [
"clf()\n",
"fp = 1.1;dw = 0.01\n",
"H0 = S1.characteristic('Hm0')[0]\n",
"St = wsm.Torsethaugen(Hm0=H0,Tp=2*pi/fp).tospecdata(np.arange(0,5+dw/2,dw)) \n",
"S1.plot()\n",
"St.plot('-.')\n",
"axis([0, 6, 0, 0.4])\n",
"show()\n"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -350,12 +658,29 @@
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": "Transformed Gaussian model compared to Gaussian model\n-------------------------------------------------------\n" "source": [
"Transformed Gaussian model compared to Gaussian model\n",
"-------------------------------------------------------\n"
]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "dt = St.sampling_period()\nva, sk, ku = St.stats_nl(moments='vsk' )\n#sa = sqrt(va)\ngh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku, ysigma=sa)\n \nysim_t = St.sim(ns=240, dt=0.5)\nxsim_t = ysim_t.copy()\nxsim_t[:,1] = gh.gauss2dat(ysim_t[:,1])\n\nts_y = wo.mat2timeseries(ysim_t)\nts_x = wo.mat2timeseries(xsim_t)\nts_y.plot_wave(sym1='r.', ts=ts_x, sym2='b', sigma=sa, nsub=5, nfig=1)\nshow()", "input": [
"dt = St.sampling_period()\n",
"va, sk, ku = St.stats_nl(moments='vsk' )\n",
"#sa = sqrt(va)\n",
"gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku, ysigma=sa)\n",
" \n",
"ysim_t = St.sim(ns=240, dt=0.5)\n",
"xsim_t = ysim_t.copy()\n",
"xsim_t[:,1] = gh.gauss2dat(ysim_t[:,1])\n",
"\n",
"ts_y = wo.mat2timeseries(ysim_t)\n",
"ts_x = wo.mat2timeseries(xsim_t)\n",
"ts_y.plot_wave(sym1='r.', ts=ts_x, sym2='b', sigma=sa, nsub=5, nfig=1)\n",
"show()"
],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@ -369,7 +694,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"collapsed": false, "collapsed": false,
"input": "", "input": [],
"language": "python", "language": "python",
"metadata": {}, "metadata": {},
"outputs": [] "outputs": []

@ -0,0 +1,73 @@
{
"metadata": {
"name": "WAFO Chapter 3"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"CHAPTER3 Demonstrates distributions of wave characteristics\n",
"=============================================================\n",
"\n",
"Chapter3 contains the commands used in Chapter3 in the tutorial.\n",
" \n",
"Some of the commands are edited for fast computation. \n",
"\n",
"Section 3.2 Estimation of wave characteristics from data\n",
"----------------------------------------------------------\n",
"Example 1\n",
"~~~~~~~~~~ "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"speed = 'fast'\n",
"#speed = 'slow'\n",
"\n",
"import wafo.data as wd\n",
"import wafo.misc as wm\n",
"import wafo.objects as wo\n",
"xx = wd.sea() \n",
"xx[:,1] = wm.detrendma(xx[:,1],len(xx))\n",
"ts = wo.mat2timeseries(xx)\n",
"Tcrcr, ix = ts.wave_periods(vh=0, pdef='c2c', wdef='tw', rate=8)\n",
"Tc, ixc = ts.wave_periods(vh=0, pdef='u2d', wdef='tw', rate=8)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"ename": "AssertionError",
"evalue": "",
"output_type": "pyerr",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mAssertionError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-12-5b70e90102e6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mxx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mwm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdetrendma\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 9\u001b[0m \u001b[0mts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mwo\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmat2timeseries\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m \u001b[0mTcrcr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mix\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mts\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwave_periods\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvh\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpdef\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'c2c'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwdef\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'tw'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrate\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 11\u001b[0m \u001b[0mTc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mixc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mts\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwave_periods\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvh\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpdef\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'u2d'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwdef\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'tw'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrate\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\pab\\workspace\\pywafo_svn\\pywafo\\src\\wafo\\objects.pyc\u001b[0m in \u001b[0;36mwave_periods\u001b[1;34m(self, vh, pdef, wdef, index, rate)\u001b[0m\n\u001b[0;32m 1980\u001b[0m \u001b[0mn\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mceil\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mrate\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1981\u001b[0m \u001b[0mti\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlinspace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1982\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstineman_interp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mti\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1983\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1984\u001b[0m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mC:\\Python27\\lib\\site-packages\\matplotlib\\mlab.pyc\u001b[0m in \u001b[0;36mstineman_interp\u001b[1;34m(xi, x, y, yp)\u001b[0m\n\u001b[0;32m 2932\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloat_\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2933\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloat_\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2934\u001b[1;33m \u001b[1;32massert\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 2935\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2936\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0myp\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mAssertionError\u001b[0m: "
]
}
],
"prompt_number": 12
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long
Loading…
Cancel
Save