Image presentations

	gmt("gmtset GMT_FFT kiss")
	
	cpos  = fitcircle("@sat_03.xyg", L=2, F=:m)
	cposx = cpos[1].data[1,1];		cposy = cpos[1].data[1,2]
	ppos  = fitcircle("@sat_03.xyg", L=2, F=:n)
	pposx = ppos[1].data[1,1];		pposy = ppos[1].data[1,2]

	# Now we use "project" to gmt 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.

	sat_pg  = gmt(@sprintf("project @sat_03.xyg -C%f/%f -T%f/%f -S -Fpz -Q", cposx, cposy, pposx, pposy));
	ship_pg = gmt(@sprintf("project @ship_03.xyg -C%f/%f -T%f/%f -S -Fpz -Q", cposx, cposy, pposx, pposy));
	sat_pg  = sat_pg[1].data;
	ship_pg = ship_pg[1].data;

	# The gmtinfo 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. 
	R = gmtinfo([sat_pg; ship_pg], I=(100,25))
	R = R[1].text[1]
	plot(sat_pg, region=R, U="L/-1.75i/-1.25i/\"Example 3a in Cookbook\"", frame=(axes=:WeSn,
		xaxis=(annot=500, ticks=100, label="Distance along great circle"), yaxis=(annot=100, ticks=25,
        label="Gravity anomaly (mGal)")), x_off=5, y_off=3.75, W=:thick, figsize=(20,12.7))
	plot!(ship_pg, S="p0.03i", savefig="example_03a", show=1)

	# 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 "histogram".

	histogram(diff(ship_pg, dims=1), W=0.1, G=:black, x_off=5, y_off=3.75, frame=0, title="Ship", U="L/-1.75i/-1.25i/\"Example 3b in Cookbook\"", figsize=7.5)
	histogram!(diff(sat_pg, dims=1), W=0.1, G=:black, x_off=12.5, frame=0, title="Ship", figsize=7.5)

	# 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 unevenly, with
	# delta-p between 0.095 and 9.017.  This means that when we want 1 km even sampling,
	# we can use "gmt 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 gmt gmtmath UPPER/LOWER.

	sampr1 = gmt("gmtmath ? -Ca -Sf -o0 UPPER CEIL =",  [ship_pg[1:1,:]; sat_pg[1:1,:]])
	sampr2 = gmt("gmtmath ? -Ca -Sf -o0 LOWER FLOOR =", [ship_pg[end:end,:]; sat_pg[end:end,:]])

	# Now we can use sampr1|2 in gmt gmtmath to make a sampling points file for gmt sample1d:
	samp_x = gmt(@sprintf("gmtmath -T%d/%d/1 -N1/0 T =", sampr1[1].data[1,1], sampr2[1].data[1,1]))

	# Now we can resample the gmt projected satellite data:
	samp_sat_pg = sample1d(sat_pg, samp_x, N=true)

	# For reasons above, we use gmt filter1d to pre-treat the ship data.  We also need to sample
	# it because of the gaps > 1 km we found. So we use gmt filter1d | gmt sample1d.  We also
	# use the -E on gmt filter1d to use the data all the way out to sampr1/sampr2 :
	t = gmt(@sprintf("filter1d -Fm1 -T%d/%d/1 -E", sampr1[1].data[1], sampr2[1].data[1]), ship_pg)
	samp_ship_pg = sample1d(t, samp_x, N=true, savefig="example_03b", show=1)

	# Now we plot them again to see if we have done the right thing:
	plot(sat_pg, region=R, U="L/-1.75i/-1.25i/\"Example 3c in Cookbook\"", frame=(axes=:WeSn,
		xaxis=(annot=500, ticks=100, label="Distance along great circle"), yaxis=(annot=100, ticks=25,
        label="Gravity anomaly (mGal)")), x_off=5, y_off=3.75, W=:thick, figsize=(20,12.7))
	plot!(samp_ship_pg, S="p0.03i", savefig="example_03c", show=1)

	# Now to do the cross-spectra, assuming that the ship is the input and the sat is the output 
	# data, we do this:
	t = [samp_ship_pg[1].data[:,2] samp_sat_pg[1].data[:,2]]
	spects = spectrum1d(t, S=256, D=1, W=true, C=true, N=true)
 
	# 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 gmt pstext and gmt psxy to put labels and legends directly on the
	# plots. For that purpose we often use -Jx1i and specify positions in inches directly:

	plot(spects[1].data[:,[1,16,17]], region=(1,1000,0,1), frame=(axes=WeSn, bg=(240,255,240), xaxis=(annot=1, ticks=3, p=[], label="Wavelength (km)"), yaxis=(annot=0.25, ticks=0.05, label="Coherency@+2@+")), 
		x_off="2.5i", y_off="1.5i", S="c0.07i", G=:purple, E="y/0.5p", figsize=("-4il", "3.75i"))
	
	text!(text_record("Coherency@+2@+"), F="+cTR+f18p,Helvetica-Bold", D="j0.1i")

	plot!(spects[1].data[:,1:3], region=(1,1000,0.1,10000), frame=(axes=WeSn, title="Ship and Satellite Gravity", bg=(240,255,240), xaxis=(annot=1, ticks=3, p=[]), yaxis=(annot=0.25, ticks=0.05, label="Power (mGal@+2@+km)")), G=:red, S="T0.07i", y_off="4.2i", E="y/0.5p", figsize=("-4il","3.75il"))

	plot!(spects[1].data[:,[1,4,5]], G=:blue, S="c0.07i", E="y/0.5p")
	text!(text_record("Input Power"), region=(0,4,0,3.75), F="+cTR+f18p,Helvetica-Bold", D="j0.1i", scale="1i")

	leg = text_record(["S 0.1i T 0.07i red - 0.3i Ship", "S 0.1i c 0.07i blue - 0.3i Satellite"])
	legend!(leg, D="jBL+w1.2i+o0.25i", F="+gwhite+pthicker", par=(FONT_ANNOT_PRIMARY="14,Helvetica-Bold",), savefig="example_03", show=1)

	# 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 "gmt 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:

	samp_ship_xw = trend1d(samp_ship_pg, F=:xw, N="p2+r")
	plot(samp_ship_pg, region=R, U="L/-1.75i/-1.25i/\"Example 3c in Cookbook\"", frame=(axes=:WeSn,
		xaxis=(annot=500, ticks=100, label="Distance along great circle"), yaxis=(annot=100, ticks=25,
        label="Gravity anomaly (mGal)")), x_off=5, y_off=3.75, S="p0.03i", figsize=(20.3,10.2))

	R = gmtinfo(samp_ship_xw, I=(100,1.1))
	plot!(samp_ship_xw, region=R[1].text[1], y_off="4.25i", frame=(axes=:Wesn, xaxis=(ticks=100,), yaxis=(annot=0.5, ticks=0.1, label=:Weight)), S="p0.03i", figsize=("8i","4.25i"), savefig="example_03d", show=1)