6.3 Spectral estimation and xy-plots

In this example we will show how to use the GMT programs fitcircle, project, sample1d, spectrum1d, psxy, and pstext. Suppose you have (lon, lat, gravity) along a satellite track in a file called sat.xyg, and (lon, lat, gravity) along a ship track in a file called ship.xyg. You want to make a cross-spectral analysis of these data. First, you will have to get the two data sets into equidistantly sampled time-series form. To do this, it will be convenient to project these along the great circle that best fits the sat track. We must use fitcircle to find this great circle and choose the L$_2$ estimates of best pole. We project the data using project to find out what their ranges are in the projected coordinate. The minmax utility will report the minimum and maximum values for multi-column ASCII tables. Use this information to select the range of the projected distance coordinate they have in common. The script prompts you for that information after reporting the values. We decide to make a file of equidistant sampling points spaced 1 km apart from -1167 to +1169, and use the UNIX utility $AWK to accomplish this step. We can then resample the projected data, and carry out the cross-spectral calculations, assuming that the ship is the input and the satellite is the output data. There are several intermediate steps that produce helpful plots showing the effect of the various processing steps (example_3[a-f].ps), while the final plot example_03.ps shows the ship and sat power in one diagram and the coherency on another diagram, both on the same page. Note the extended use of pstext and psxy to put labels and legends directly on the plots. For that purpose we often use -Jx1i and specify positions in inches directly. Thus, the complete automated script reads:

fitcircle sat.xyg -L2 >! report
set cpos = `grep "L2 Average Position" report` 
set ppos = `grep "L2 N Hemisphere" report` 
project  sat.xyg -C$cpos[1]/$cpos[2] -T$ppos[1]/$ppos[2] -S -Fpz -Q >! sat.pg
project ship.xyg -C$cpos[1]/$cpos[2] -T$ppos[1]/$ppos[2] -S -Fpz -Q >! ship.pg
set plotr = `cat sat.pg ship.pg | minmax -I100/25 -C`
gmtset MEASURE_UNIT INCH
psxy -R$plotr[1]/$plotr[2]/$plotr[3]/$plotr[4] -U/-1.75i/-1.25i/"Example 3a in Cookbook" \
   -JX8i/5i -X2i -Y1.5i -K -W1p sat.pg \
   -Ba500f100:"Distance along great circle":/a100f25:"Gravity anomaly (mGal)":WeSn >! example_03a.ps
psxy -R -JX -O -Sp0.03i ship.pg >> example_03a.ps
$AWK '{ if (NR > 1) print $1 - last1; last1 = $1; }' ship.pg | pshistogram  -W0.1 -G0 -JX3i -K \
   -X2i -Y1.5i -B:."Ship": -U/-1.75i/-1.25i/"Example 3b in Cookbook" >! example_03b.ps
$AWK '{ if (NR > 1) print $1 - last1; last1 = $1; }' sat.pg  | pshistogram  -W0.1 -G0 -JX3i -O \
   -X5i -B:."Sat": >> example_03b.ps
head -1 ship.pg >! ship.pg.extr
head -1 sat.pg >! sat.pg.extr
paste ship.pg.extr sat.pg.extr | $AWK '{ if ($1 > $3) print int($1); else print int($3); }' \
   >! sampr1
tail -1 ship.pg >! ship.pg.extr
tail -1 sat.pg >! sat.pg.extr 
paste ship.pg.extr sat.pg.extr | $AWK '{ if ($1 < $3) print int($1); else print int($3); }' \
   >! sampr2
set sampr = `paste sampr1 sampr2`
$AWK 'BEGIN { for (i = '$sampr[1]'; i <= '$sampr[2]'; i++) print i }' /dev/null >! samp.x
sample1d sat.pg -Nsamp.x >! samp_sat.pg
filter1d ship.pg -Fm1 -T$sampr[1]/$sampr[2]/1 -E | sample1d -Nsamp.x >! samp_ship.pg
psxy -R$plotr[1]/$plotr[2]/$plotr[3]/$plotr[4] -JX8i/5i -X2i -Y1.5i -K -W1p samp_sat.pg \
   -Ba500f100:"Distance along great circle":/a100f25:"Gravity anomaly (mGal)":WeSn  \
   -U/-1.75i/-1.25i/"Example 3c in Cookbook" >! example_03c.ps
psxy -R -JX -O -Sp0.03i samp_ship.pg >> example_03c.ps
paste samp_ship.pg samp_sat.pg | cut -f2,4 | spectrum1d -S256 -D1 -W -C >& /dev/null
psxy spectrum.coh -Ba1f3p:"Wavelength (km)":/a0.25f0.05:"Coherency@+2@+":WeSn -JX-4il/3.75i \
   -R1/1000/0/1 -U/-2.25i/-1.25i/"Example 3d in Cookbook" -P -K -X2.5i -Sc0.07i -G0 \
   -Ey/2 -Y1.5i >! example_03.ps
echo "3.85 3.6 18 0.0 1 11 Coherency@+2@+" | pstext -R0/4/0/3.75 -Jx1i -O -K >> example_03.ps
cat << END >! box.d
2.375    3.75
2.375    3.25
4    3.25
END
psxy -R -Jx -O -K -W1.5p box.d >> example_03.ps
psxy -Ba1f3p/a1f3p:"Power (mGal@+2@+km)"::."Ship and Satellite Gravity":WeSn spectrum.xpower \
   -St0.07i -O -R1/1000/0.1/10000 -JX-4il/3.75il -Y4.2i -K -Ey/2 >> example_03.ps
psxy spectrum.ypower -R -JX -O -K -G0 -Sc0.07i -Ey/2 >> example_03.ps
echo "3.9 3.6 18 0.0 1 11 Input Power" | pstext -R0/4/0/3.75 -Jx -O -K >> example_03.ps
psxy -R -Jx -O -K -W1.5p box.d >> example_03.ps
psxy -R -Jx -O -K -G240 -L -W1.5p << END >> example_03.ps
0.25    0.25
1.4    0.25
1.4    0.9
0.25    0.9
END
echo "0.4 0.7" | psxy -R -Jx -O -K -St0.07i -G0 >> example_03.ps
echo "0.5 0.7 14 0.0 1 5 Ship" | pstext -R -Jx -O -K >> example_03.ps
echo "0.4 0.4" | psxy -R -Jx -O -K -Sc0.07i -G0 >> example_03.ps
echo "0.5 0.4 14 0.0 1 5 Satellite" | pstext -R -Jx -O >> example_03.ps
trend1d -Fxw -N2r samp_ship.pg >! samp_ship.xw
psxy -R$plotr[1]/$plotr[2]/$plotr[3]/$plotr[4] -JX8i/4i -X2i -Y1.5i -K -Sp0.03i \
   -Ba500f100:"Distance along great circle":/a100f25:"Gravity anomaly (mGal)":WeSn \
   -U/-1.75i/-1.25i/"Example 3e in Cookbook" samp_ship.pg >! example_03d.ps
psxy -R$plotr[1]/$plotr[2]/0/1.1 -JX8i/1.1i -O -Y4.25i -Bf100/a0.5f0.1:"Weight":Wesn -Sp0.03i \
   samp_ship.xw >> example_03d.ps
trend1d -Fxrw -N2r samp_ship.pg | $AWK '{ if ($3 > 0.6) print $1, $2 }' | sample1d -Nsamp.x >! \
   samp2_ship.pg
trend1d -Fxrw -N2r samp_sat.pg  | $AWK '{ if ($3 > 0.6) print $1, $2 }' | sample1d -Nsamp.x >! \
   samp2_sat.pg
set plotr = `cat samp2_sat.pg samp2_ship.pg | minmax -I100/25 -C`
psxy -R$plotr[1]/$plotr[2]/$plotr[3]/$plotr[4] -JX8i/5i -X2i -Y1.5i -K -W1p \
   -Ba500f100:"Distance along great circle":/a50f25:"Gravity anomaly (mGal)":WeSn \
   -U/-1.75i/-1.25i/"Example 3f in Cookbook" samp2_sat.pg >! example_03e.ps
psxy -R -JX -O -Sp0.03i samp2_ship.pg >> example_03e.ps
paste samp2_ship.pg samp2_sat.pg | cut -f2,4 | spectrum1d -S256 -D1 -W -C >& /dev/null
psxy spectrum.coh -Ba1f3p:"Wavelength (km)":/a0.25f0.05:"Coherency@+2@+":WeSn -JX-4il/3.75i \
   -R1/1000/0/1 -U/-2.25i/-1.25i/"Example 3g in Cookbook" -P -K -X2.5i -Sc0.07i -G0 \
   -Ey/2 -Y1.5i >! example_03f.ps
echo "3.85 3.6 18 0.0 1 11 Coherency@+2@+" | pstext -R0/4/0/3.75 -Jx -O -K >> example_03f.ps
cat << END >! box.d
2.375    3.75
2.375    3.25
4    3.25
END
psxy -R -Jx -O -K -W1.5p box.d >> example_03f.ps
psxy -Ba1f3p/a1f3p:"Power (mGal@+2@+km)"::."Ship and Satellite Gravity":WeSn spectrum.xpower \
   -St0.07i -O -R1/1000/0.1/10000 -JX-4il/3.75il -Y4.2i -K -Ey/2 >> example_03f.ps
psxy spectrum.ypower -R -JX -O -K -G0 -Sc0.07i -Ey/2 >> example_03f.ps
echo "3.9 3.6 18 0.0 1 11 Input Power" | pstext -R0/4/0/3.75 -Jx -O -K >> example_03f.ps
psxy -R -Jx -O -K -W1.5p box.d >> example_03f.ps
psxy -R -Jx -O -K -G240 -L -W1.5p << END >> example_03f.ps
0.25    0.25
1.4    0.25
1.4    0.9
0.25    0.9
END
echo "0.4 0.7" | psxy -R -Jx -O -K -St0.07i -G0 >> example_03f.ps
echo "0.5 0.7 14 0.0 1 5 Ship" | pstext -R -Jx -O -K >> example_03f.ps
echo "0.4 0.4" | psxy -R -Jx -O -K -Sc0.07i -G0 >> example_03f.ps
echo "0.5 0.4 14 0.0 1 5 Satellite" | pstext -R -Jx -O >> example_03f.ps
\rm -f box.d report samp* *.pg *.extr spectrum.* .gmtcommands

The final illustration (Figure 6.3) shows that the ship gravity anomalies have more power than altimetry derived gravity for short wavelengths and that the coherency between the two signals improves dramatically for wavelengths $>$ 20 km.

Figure 6.3: Spectral estimation and $x/y$-plots