/usr/share/doc/gmt-examples/examples/ex03/job03.sh is in gmt-examples 4.5.11-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 | #!/bin/bash
# GMT EXAMPLE 03
# $Id: job03.sh 9545 2011-07-27 19:31:54Z pwessel $
#
# 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":
#
. ../functions.sh
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. This is a job for gmtmath UPPER/LOWER.
#
head -1 ship.pg > tmp
head -1 sat.pg >> tmp
sampr1=`gmtmath tmp -Ca -Sf -F0 UPPER CEIL =`
tail -1 ship.pg > tmp
tail -1 sat.pg >> tmp
sampr2=`gmtmath tmp -Ca -Sf -F0 LOWER FLOOR =`
#
# Now we can use sampr1|2 in gmtmath to make a sampling points file for sample1d:
gmtmath -T$sampr1/$sampr2/1 -N1/0 T = 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 tmp samp* *.pg *.extr spectrum.* .gmt*
|