|
| 1 | +# (3) Spectral estimation and xy-plots |
| 2 | + |
| 3 | +In this example we will show how to use the GMT programs \myreflink{fitcircle}, \myreflink{project}, |
| 4 | +\myreflink{sample1d}, \myreflink{spectrum1d}, and \myreflink{plot}. Suppose you have (lon, lat, gravity) |
| 5 | +along a satellite track in a file called ``sat.xyg``, and (lon, lat, gravity) along a ship track in a |
| 6 | +file called ``ship.xyg``. You want to make a cross-spectral analysis of these data. First, you will |
| 7 | +have to get the two data sets into equidistantly sampled time-series form. To do this, it will be |
| 8 | +convenient to project these along the great circle that best fits the sat track. We must use |
| 9 | +\myreflink{fitcircle} to find this great circle and choose the L2 estimates of best pole. We project |
| 10 | +the data using \myreflink{project} to find out what their ranges are in the projected coordinate. |
| 11 | +The script computes the common range. We can then resample the projected data, and carry out the |
| 12 | +cross-spectral calculations, assuming that the ship is the input and the satellite is the output data. |
| 13 | +The final plot example_03.ps shows the ship and sat power in one diagram and the coherency on another |
| 14 | +diagram, both on the same page, with an auto-generated legend. Thus, the complete automated script reads: |
| 15 | + |
| 16 | +The final illustration shows that the ship gravity anomalies have more power than altimetry derived gravity |
| 17 | +for short wavelengths and that the coherency between the two signals improves dramatically for wavelengths > 20 km. |
| 18 | + |
| 19 | +\begin{examplefig}{} |
| 20 | +```julia |
| 21 | +using GMT |
| 22 | +GMT.resetGMT() # hide |
| 23 | + |
| 24 | +# First, we use "fitcircle" to find the parameters of a great circle |
| 25 | +# most closely fitting the x,y points in "sat_03.txt": |
| 26 | +cpos = fitcircle("@sat_03.txt", norm=2, coordinates=:mean, par=(IO_COL_SEPARATOR="/",)) |
| 27 | +ppos = fitcircle("@sat_03.txt", norm=2, coordinates=:north, par=(IO_COL_SEPARATOR="/",)) |
| 28 | + |
| 29 | +# Now we use "project" to project the data in both sat_03.txt and ship_03.txt |
| 30 | +# into data.pg, where g is the same and p is the oblique longitude around |
| 31 | +# the great circle. We use "km" to get the p distance in kilometers, and "sort" |
| 32 | +# to sort the output into increasing p values. |
| 33 | +sat_pg = project("@sat_03.txt", origin=cpos, pole=ppos, sort=true, outvars=:pz, km=true) |
| 34 | +ship_pg = project("@ship_03.txt", origin=cpos, pole=ppos, sort=true, outvars=:pz, km=true) |
| 35 | + |
| 36 | +# Get the common extrema values rounded to 1 km |
| 37 | +bounds = round.([max(sat_pg.bbox[1], ship_pg.bbox[1]) min(sat_pg.bbox[2], ship_pg.bbox[2])], digits=0) |
| 38 | + |
| 39 | +samp_x = collect(bounds[1]:bounds[2]) |
| 40 | +# Now we can resample the gmt projected satellite data: |
| 41 | +samp_sat_pg = sample1d(sat_pg, range=samp_x) |
| 42 | + |
| 43 | +# For reasons above, we use filter1d to pre-treat the ship data. We also need to sample |
| 44 | +# it because of the gaps > 1 km we found. So we use filter1d | sample1d. We also |
| 45 | +# use the "ends" on filter1d to use the data all the way out to bounds : |
| 46 | +D = filter1d(ship_pg, filter=:m1, range=@sprintf("%g/%g/1", bounds...), ends=true) |
| 47 | +samp_ship_pg = sample1d(D, range=samp_x) |
| 48 | + |
| 49 | +# Now to do the cross-spectra, assuming that the ship is the input and the sat is the output |
| 50 | +# data, we do this: |
| 51 | +D = spectrum1d([samp_ship_pg[:,2] samp_sat_pg[:,2]], size=256, sample_dist=1, wavelength=true, outputs=(:xpower, :ypower, :coherence)) |
| 52 | + |
| 53 | +# Time to plot spectra |
| 54 | +subplot(grid=(2,1), margins="0.3c", col_axes=(bott=true, label="Wavelength (km)"), |
| 55 | + autolabel=(anchor=:TR, offset="8p"), title="Ship and Satellite Gravity", |
| 56 | + panel_size=10, frame=(axes="WeSn", bg="240/255/240"), Vd=1) |
| 57 | + gmtset(FONT_TAG="18p,Helvetica-Bold",) |
| 58 | + |
| 59 | + subplot(:set, panel=(1,1), autolabel="Input Power") |
| 60 | + plot(D[:Wavelength, :Xpower, :σ_Xpow], proj=:loglog, flip_axes=:x, xaxis=(annot=1, ticks=3, scale=:pow), |
| 61 | + yaxis=(annot=1, ticks=3, scale=:pow, label="Power (mGal@+2@+km)"), fill="red", marker=:Triangle, |
| 62 | + ms="5p", region=(1,1000,0.1,10000), error_bars=(y=true, pen=0.5), legend="Ship") |
| 63 | + |
| 64 | + plot(D[:Wavelength, :Ypower, :σ_Ypow], fill="blue", marker=:circ, ms="5p", |
| 65 | + error_bars=(y=true, pen=0.5), legend="Satellite") |
| 66 | + |
| 67 | + legend(pos=(anchor="BL", offset=0.5), box="+gwhite+pthicker", par=(FONT_ANNOT_PRIMARY="14p,Helvetica-Bold",)) |
| 68 | + subplot(:set, panel=(2,1), autolabel="Coherency@+2@+") |
| 69 | + |
| 70 | + plot(D[:Wavelength, :Coher, :σ_Coher], region=(1,1000,0,1), proj=:logx, flip_axes=:x, |
| 71 | + xaxis=(annot=1, ticks=3, scale=:pow, label="Wavelength (km)"), |
| 72 | + yaxis=(annot=0.25, ticks=0.05, label="Coherency@+2@+"), marker=:circ, ms="5p", |
| 73 | + fill=:purple, error_bars=(y=true, pen=0.5)) |
| 74 | +subplot(:show) |
| 75 | +``` |
| 76 | +\end{examplefig} |
0 commit comments