<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN">
<!--Converted with LaTeX2HTML 2002-2-1 (1.71)
original version by: Nikos Drakos, CBLU, University of Leeds
* revised and updated by: Marcus Hennecke, Ross Moore, Herb Swan
* with significant contributions from:
Jens Lippmann, Marek Rouchal, Martin Wilck and others -->
<HTML>
<HEAD>
<TITLE>7.3 Spectral estimation and xy-plots</TITLE>
<META NAME="description" CONTENT="7.3 Spectral estimation and xy-plots">
<META NAME="keywords" CONTENT="GMT_Docs">
<META NAME="resource-type" CONTENT="document">
<META NAME="distribution" CONTENT="global">
<META NAME="Generator" CONTENT="LaTeX2HTML v2002-2-1">
<META HTTP-EQUIV="Content-Style-Type" CONTENT="text/css">
<LINK REL="STYLESHEET" HREF="GMT_Docs.css">
<LINK REL="next" HREF="node125.html">
<LINK REL="previous" HREF="node123.html">
<LINK REL="up" HREF="node121.html">
<LINK REL="next" HREF="node125.html">
</HEAD>
<BODY bgcolor="#ffffff">
<!--Navigation Panel-->
<A NAME="tex2html4392"
HREF="node125.html">
<IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="next.png"></A>
<A NAME="tex2html4386"
HREF="node121.html">
<IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="up.png"></A>
<A NAME="tex2html4380"
HREF="node123.html">
<IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="prev.png"></A>
<A NAME="tex2html4388"
HREF="node1.html">
<IMG WIDTH="65" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="contents" SRC="contents.png"></A>
<A NAME="tex2html4390"
HREF="node256.html">
<IMG WIDTH="43" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="index" SRC="index.png"></A>
<BR>
<B> Next:</B> <A NAME="tex2html4393"
HREF="node125.html">7.4 A 3-D perspective</A>
<B> Up:</B> <A NAME="tex2html4387"
HREF="node121.html">7. Creating GMT Graphics</A>
<B> Previous:</B> <A NAME="tex2html4381"
HREF="node123.html">7.2 Image presentations</A>
<B> <A NAME="tex2html4389"
HREF="node1.html">Contents</A></B>
<B> <A NAME="tex2html4391"
HREF="node256.html">Index</A></B>
<BR>
<BR>
<!--End of Navigation Panel-->
<H1><A NAME="SECTION001530000000000000000"></A>
<A NAME="24310"></A>
<A NAME="24311"></A>
<BR>
7.3 Spectral estimation and xy-plots
</H1>
<P>
In this example we will show how to use the <A NAME="tex2html1341"
HREF="http://gmt.soest.hawaii.edu"><B>GMT</B></A> programs
<A NAME="tex2html1342"
HREF="../man/fitcircle.html"><I><B>fitcircle</B></I></A><A NAME="24888"></A>, <A NAME="tex2html1343"
HREF="../man/project.html"><I><B>project</B></I></A><A NAME="24893"></A>, <A NAME="tex2html1344"
HREF="../man/sample1d.html"><I><B>sample1d</B></I></A><A NAME="24898"></A>,
<A NAME="tex2html1345"
HREF="../man/spectrum1d.html"><I><B>spectrum1d</B></I></A><A NAME="24903"></A>, <A NAME="tex2html1346"
HREF="../man/psxy.html"><I><B>psxy</B></I></A><A NAME="24908"></A>, and <A NAME="tex2html1347"
HREF="../man/pstext.html"><I><B>pstext</B></I></A><A NAME="24913"></A>.
Suppose you have (lon, lat, gravity) along a satellite track
in a file called <U>sat.xyg</U>, and (lon, lat, gravity)
along a ship track in a file called <U>ship.xyg</U>.
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 <A NAME="tex2html1348"
HREF="../man/fitcircle.html"><I><B>fitcircle</B></I></A><A NAME="24920"></A> to find this great circle and choose
the L<IMG
WIDTH="11" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
SRC="img26.png"
ALT="$_2$"> estimates of best pole. We project the data using
<A NAME="tex2html1349"
HREF="../man/project.html"><I><B>project</B></I></A><A NAME="24925"></A> to find out what their ranges are in the projected
coordinate. The <A NAME="tex2html1350"
HREF="../man/minmax.html"><I><B>minmax</B></I></A><A NAME="24930"></A> 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 <I>UNIX</I> utility <I>$AWK</I><A NAME="24935"></A> 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 (<U>example_03[a-f].ps</U>), while the final plot
<U>example_03.ps</U> 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
<A NAME="tex2html1351"
HREF="../man/pstext.html"><I><B>pstext</B></I></A><A NAME="24941"></A> and <A NAME="tex2html1352"
HREF="../man/psxy.html"><I><B>psxy</B></I></A><A NAME="24946"></A> to put labels and legends directly on
the plots. For that purpose we often use <B>-Jx</B>1i and specify
positions in inches directly. Thus, the complete automated
script reads:
<P>
<BR CLEAR="ALL">
<HR>
<BR>
<PRE>#!/bin/sh
# GMT EXAMPLE 03
#
# Purpose: Resample track data, do spectral analysis, and plot
# GMT progs: filter1d, fitcircle, minmax, project, sample1d
# GMT progs: spectrum1d, trend1d, pshistogram, psxy, pstext
# Unix progs: $AWK, cat, echo, head, paste, rm, tail
#
# This example begins with data files "ship.xyg" and "sat.xyg" which
# are measurements of a quantity "g" (a "gravity anomaly" which is an
# anomalous increase or decrease in the magnitude of the acceleration
# of gravity at sea level). g is measured at a sequence of points "x,y"
# which in this case are "longitude,latitude". The "sat.xyg" data were
# obtained by a satellite and the sequence of points lies almost along
# a great circle. The "ship.xyg" data were obtained by a ship which
# tried to follow the satellite's path but deviated from it in places.
# Thus the two data sets are not measured at the same of points,
# and we use various GMT tools to facilitate their comparison.
# The main illustration (example_03.ps) are accompanied with 5 support
# plots (03a-f) showing data distributions and various intermediate steps.
#
# First, we use "fitcircle" to find the parameters of a great circle
# most closely fitting the x,y points in "sat.xyg":
#
ps=example_03.ps
fitcircle sat.xyg -L2 > report
cposx=`grep "L2 Average Position" report | cut -f1`
cposy=`grep "L2 Average Position" report | cut -f2`
pposx=`grep "L2 N Hemisphere" report | cut -f1`
pposy=`grep "L2 N Hemisphere" report | cut -f2`
#
# Now we use "project" to project the data in both sat.xyg and ship.xyg
# into data.pg, where g is the same and p is the oblique longitude around
# the great circle. We use -Q to get the p distance in kilometers, and -S
# to sort the output into increasing p values.
#
project sat.xyg -C$cposx/$cposy -T$pposx/$pposy -S -Fpz -Q > sat.pg
project ship.xyg -C$cposx/$cposy -T$pposx/$pposy -S -Fpz -Q > ship.pg
#
# The minmax utility will report the minimum and maximum values for all columns.
# We use this information first with a large -I value to find the appropriate -R
# to use to plot the .pg data.
#
plotr=`cat sat.pg ship.pg | minmax -I100/25`
psxy $plotr -U/-1.75i/-1.25i/"Example 3a in Cookbook" \
-Ba500f100:"Distance along great circle":/a100f25:"Gravity anomaly (mGal)":WeSn \
-JX8i/5i -X2i -Y1.5i -K -Wthick sat.pg > example_03a.ps
psxy -R -JX -O -Sp0.03i ship.pg >> example_03a.ps
#
# From this plot we see that the ship data have some "spikes" and also greatly
# differ from the satellite data at a point about p ~= +250 km, where both of
# them show a very large anomaly.
#
# To facilitate comparison of the two with a cross-spectral analysis using "spectrum1d",
# we resample both data sets at intervals of 1 km. First we find out how the data are
# typically spaced using $AWK to get the delta-p between points and view it with
# "pshistogram".
#
$AWK '{ if (NR > 1) print $1 - last1; last1=$1; }' ship.pg | pshistogram -W0.1 -Gblack -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 -Gblack -JX3i -O \
-X5i -B:."Sat": >> example_03b.ps
#
# This experience shows that the satellite values are spaced fairly evenly, with
# delta-p between 3.222 and 3.418. The ship values are spaced quite unevelnly, with
# delta-p between 0.095 and 9.017. This means that when we want 1 km even sampling,
# we can use "sample1d" to interpolate the sat data, but the same procedure applied
# to the ship data could alias information at shorter wavelengths. So we have to use
# "filter1d" to resample the ship data. Also, since we observed spikes in the ship
# data, we use a median filter to clean up the ship values. We will want to use "paste"
# to put the two sampled data sets together, so they must start and end at the same
# point, without NaNs. So we want to get a starting and ending point which works for
# both of them. Thus we need to start at max( min(ship.p), min(sat.p) ) and end
# conversely. "minmax" can't do this easily since it will return min( min(), min() ),
# so we do a little head, paste $AWK to get what we want.
#
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
sampr1=`cat sampr1`
sampr2=`cat sampr2`
#
# Now we can use sampr in $AWK to make a sampling points file for sample1d:
$AWK 'BEGIN { for (i='$sampr1'; i <= '$sampr2'; i++) print i }' /dev/null > samp.x
#
# Now we can resample the projected satellite data:
#
sample1d sat.pg -Nsamp.x > samp_sat.pg
#
# For reasons above, we use filter1d to pre-treat the ship data. We also need to sample it
# because of the gaps > 1 km we found. So we use filter1d | sample1d. We also use the -E
# on filter1d to use the data all the way out to sampr1/sampr2 :
#
filter1d ship.pg -Fm1 -T$sampr1/$sampr2/1 -E | sample1d -Nsamp.x > samp_ship.pg
#
# Now we plot them again to see if we have done the right thing:
#
psxy $plotr -JX8i/5i -X2i -Y1.5i -K -Wthick 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
#
# Now to do the cross-spectra, assuming that the ship is the input and the sat is the output
# data, we do this:
#
paste samp_ship.pg samp_sat.pg | cut -f2,4 | spectrum1d -S256 -D1 -W -C > /dev/null
#
# Now we want to plot the spectra. The following commands will plot 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:
#
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 3 in Cookbook" -P -K -X2.5i -Sc0.07i -Gblack \
-Ey/0.5p -Y1.5i > $ps
echo "3.85 3.6 18 0.0 1 TR Coherency@+2@+" | pstext -R0/4/0/3.75 -Jx1i -O -K >> $ps
cat > box.d << END
2.375 3.75
2.375 3.25
4 3.25
END
psxy -R -Jx -O -K -Wthicker box.d >> $ps
psxy -Ba1f3p/a1f3p:"Power (mGal@+2@+km)"::."Ship and Satellite Gravity":WeSn spectrum.xpower \
-Gblack -ST0.07i -O -R1/1000/0.1/10000 -JX-4il/3.75il -Y4.2i -K -Ey/0.5p >> $ps
psxy spectrum.ypower -R -JX -O -K -Gblack -Sc0.07i -Ey/0.5p >> $ps
echo "3.9 3.6 18 0.0 1 TR Input Power" | pstext -R0/4/0/3.75 -Jx -O -K >> $ps
psxy -R -Jx -O -K -Wthicker box.d >> $ps
psxy -R -Jx -O -K -Glightgray -L -Wthicker >> $ps << END
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 -Gblack >> $ps
echo "0.5 0.7 14 0.0 1 ML Ship" | pstext -R -Jx -O -K >> $ps
echo "0.4 0.4" | psxy -R -Jx -O -K -Sc0.07i -Gblack >> $ps
echo "0.5 0.4 14 0.0 1 ML Satellite" | pstext -R -Jx -O >> $ps
#
# Now we wonder if removing that large feature at 250 km would make any difference.
# We could throw away a section of data with $AWK or sed or head and tail, but we
# demonstrate the use of "trend1d" to identify outliers instead. We will fit a
# straight line to the samp_ship.pg data by an iteratively-reweighted method and
# save the weights on output. Then we will plot the weights and see how things
# look:
#
trend1d -Fxw -N2r samp_ship.pg > samp_ship.xw
psxy $plotr -JX8i/4i -X2i -Y1.5i -K -Sp0.03i \
-Ba500f100:"Distance along great circle":/a100f25:"Gravity anomaly (mGal)":WeSn \
-U/-1.75i/-1.25i/"Example 3d in Cookbook" samp_ship.pg > example_03d.ps
plotr=`minmax samp_ship.xw -I100/1.1`
psxy $plotr -JX8i/1.1i -O -Y4.25i -Bf100/a0.5f0.1:"Weight":Wesn -Sp0.03i samp_ship.xw \
>> example_03d.ps
#
# From this we see that we might want to throw away values where w < 0.6. So we try that,
# and this time we also use trend1d to return the residual from the model fit (the
# de-trended data):
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
#
# We plot these to see how they look:
#
plotr=`cat samp2_sat.pg samp2_ship.pg | minmax -I100/25`
psxy $plotr -JX8i/5i -X2i -Y1.5i -K -Wthick \
-Ba500f100:"Distance along great circle":/a50f25:"Gravity anomaly (mGal)":WeSn \
-U/-1.75i/-1.25i/"Example 3e in Cookbook" samp2_sat.pg > example_03e.ps
psxy -R -JX -O -Sp0.03i samp2_ship.pg >> example_03e.ps
#
# Now we do the cross-spectral analysis again. Comparing this plot (example_03e.ps) with
# the previous one (example_03d.ps) we see that throwing out the large feature has reduced
# the power in both data sets and reduced the coherency at wavelengths between 20--60 km.
#
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 3f in Cookbook" -P -K -X2.5i -Sc0.07i -Gblack \
-Ey/0.5p -Y1.5i > example_03f.ps
echo "3.85 3.6 18 0.0 1 TR Coherency@+2@+" | pstext -R0/4/0/3.75 -Jx -O -K >> example_03f.ps
cat > box.d << END
2.375 3.75
2.375 3.25
4 3.25
END
psxy -R -Jx -O -K -Wthicker 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/0.5p >> example_03f.ps
psxy spectrum.ypower -R -JX -O -K -Gblack -Sc0.07i -Ey/0.5p >> example_03f.ps
echo "3.9 3.6 18 0.0 1 TR Input Power" | pstext -R0/4/0/3.75 -Jx -O -K >> example_03f.ps
psxy -R -Jx -O -K -Wthicker box.d >> example_03f.ps
psxy -R -Jx -O -K -Glightgray -L -Wthicker >> example_03f.ps << END
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 -Gblack >> example_03f.ps
echo "0.5 0.7 14 0.0 1 ML Ship" | pstext -R -Jx -O -K >> example_03f.ps
echo "0.4 0.4" | psxy -R -Jx -O -K -Sc0.07i -Gblack >> example_03f.ps
echo "0.5 0.4 14 0.0 1 ML Satellite" | pstext -R -Jx -O >> example_03f.ps
#
rm -f box.d report samp* *.pg *.extr spectrum.* .gmt*
</PRE>
<BR CLEAR="ALL">
<HR>
<P>
The final illustration (Figure <A HREF="#fig:GMT_example_03">7.3</A>) 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 <IMG
WIDTH="16" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
SRC="img13.png"
ALT="$>$"> 20 km.
<DIV ALIGN="CENTER"><A NAME="fig:GMT_example_03"></A><A NAME="24953"></A>
<TABLE>
<CAPTION ALIGN="BOTTOM"><STRONG>Figure 7.3:</STRONG>
Spectral estimation and <IMG
WIDTH="26" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
SRC="img6.png"
ALT="$x/y$">-plots.</CAPTION>
<TR><TD>
<DIV ALIGN="CENTER"><IMG
WIDTH="333" HEIGHT="554" ALIGN="BOTTOM" BORDER="0"
SRC="img157.png"
ALT="\includegraphics[scale=0.5]{scripts/example_03}"></DIV></TD></TR>
</TABLE>
</DIV>
<P>
<A NAME="24333"></A>
<A NAME="24334"></A>
<P>
<HR>
<!--Navigation Panel-->
<A NAME="tex2html4392"
HREF="node125.html">
<IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="next.png"></A>
<A NAME="tex2html4386"
HREF="node121.html">
<IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="up.png"></A>
<A NAME="tex2html4380"
HREF="node123.html">
<IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="prev.png"></A>
<A NAME="tex2html4388"
HREF="node1.html">
<IMG WIDTH="65" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="contents" SRC="contents.png"></A>
<A NAME="tex2html4390"
HREF="node256.html">
<IMG WIDTH="43" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="index" SRC="index.png"></A>
<BR>
<B> Next:</B> <A NAME="tex2html4393"
HREF="node125.html">7.4 A 3-D perspective</A>
<B> Up:</B> <A NAME="tex2html4387"
HREF="node121.html">7. Creating GMT Graphics</A>
<B> Previous:</B> <A NAME="tex2html4381"
HREF="node123.html">7.2 Image presentations</A>
<B> <A NAME="tex2html4389"
HREF="node1.html">Contents</A></B>
<B> <A NAME="tex2html4391"
HREF="node256.html">Index</A></B>
<!--End of Navigation Panel-->
<ADDRESS>
Paul Wessel
2010-07-14
</ADDRESS>
</BODY>
</HTML>
distrib
>
Fedora
>
14
>
x86_64
>
media
>
os
>
by-pkgid
>
4b194777eaf705e440bb2ce282d32772
>
files
>
547
