6.18 Volumes and Spatial Selections

To demonstrate potential usage of the new programs grdvolume and gmtselect we extract a subset of the Sandwell & Smith altimetric gravity field^6.3† for the northern Pacific and decide to isolate all seamounts that (1) exceed 50 mGal in amplitude and (2) are within 200 km of the Pratt seamount. We do this by dumping the 50 mGal contours to disk, then making a simple $AWK script center.awk that returns the mean location of the points making up each closed polygon, and then pass these locations to gmtselect which retains only the points within 200 km of Pratt. We then mask out all the data outside this radius and use grdvolume to determine the combined area and volumes of the chosen seamounts.

gmtset ELLIPSOID Sphere
echo "-142.65    56.25" >! pratt.d
makecpt -Crainbow -T-60/60/10 -Z >! grav.cpt
grdgradient AK_gulf_grav.grd -Nt1 -A45 -GAK_gulf_grav_i.grd
grdimage AK_gulf_grav.grd -IAK_gulf_grav_i.grd -JM5.5i -Cgrav.cpt -B2f1 -P -K -X1.5i -Y5.85i >! example_18.ps
pscoast -R-149/-135/52.5/58 -JM -O -K -Di -G160 -W0.25p >> example_18.ps
psscale -D2.75i/-0.4i/4i/0.15ih -Cgrav.cpt -B20f10/:mGal: -O -K >> example_18.ps
$AWK '{print $1, $2, 12, 0, 1, "LB", "Pratt"}' pratt.d | pstext -R -JM -O -K -D0.1i/0.1i \
   >> example_18.ps
$AWK '{print $1, $2, 0, 200, 200}' pratt.d | psxy -R -JM -O -K -SE -W0.25p >> example_18.ps
grdcontour AK_gulf_grav.grd -JM -C20 -B2f1WSEn -O -K -Y-4.85i -U/-1.25i/-0.75i/"Example 18 in Cookbook" \
   >> example_18.ps
grdcontour AK_gulf_grav.grd -JM -C10 -L49/51 -O -K -Dsm -Wc0.75p/0/255/0 >> example_18.ps
pscoast -R -JM -O -K -Di -G160 -W0.25p >> example_18.ps
$AWK '{print $1, $2, 0, 200, 200}' pratt.d | psxy -R -JM -O -K -SE -W0.25p >> example_18.ps
\rm -f sm_*[0-9].xyz    # Only consider closed contours
cat << EOF >! center.awk
BEGIN {
    x = 0
    y = 0
    n = 0
}
{
    x += \$1
    y += \$2
    n++
}
END {
    print x/n, y/n
}
EOF
\rm -f centers.d
foreach file (sm_*.xyz)
    $AWK -f center.awk $file >>! centers.d
end
gmtselect -R -JM -C200/pratt.d centers.d >! $$
psxy $$ -R -JM -O -K -SC0.04i -G255/0/0 -W0.25p >> example_18.ps
psxy -R -JM -O -K -ST0.1i -G255/255/0 -W0.25p pratt.d >> example_18.ps
grdmath -R -I2m -F -142.65 56.25 GDIST = mask.grd
grdclip mask.grd -Sa200/NaN -Sb200/1 -Gmask.grd
grdmath AK_gulf_grav.grd mask.grd MUL = tmp.grd
set info = `grdvolume tmp.grd -C50 -Sk`
psxy -R -JM -A -O -K -L -W0.75p -G255 << EOF >> example_18.ps
-148.5 52.75
-140.5    52.75
-140.5    53.75
-148.5    53.75
EOF
pstext -R -JM -O << EOF >> example_18.ps
-148 53.08 14 0 1 LM Areas: $info[2] km@+2@+
-148 53.42 14 0 1 LM Volumes: $info[3] mGal\264km@+2@+
EOF
\rm -f $$ grav.cpt sm_*.xyz *_i.grd tmp.grd mask.grd pratt.d center* .gmt*

Our illustration is presented in Figure 6.18.

Figure 6.18: Volumes and geo-spatial selections