Next: Color patterns on maps Up: Cook-book Previous: Images clipped by coastlines

#
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 field6.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.

    #!/bin/csh
    #        GMT EXAMPLE 18
    #
    #        @(#)job18.csh    1.12  03/11/99
    #
    # Purpose:    Illustrates volumes of grids inside contours and spatial
    #        selection of data
    # GMT progs:    gmtset, gmtselect, grdclip, grdcontour, grdgradient, grdimage,
    #        grdmath, grdvolume, makecpt, pscoast, psscale, pstext, psxy
    # Unix progs:    $AWK, cat, rm
    #
    # Get Sandwell/Smith gravity for the region
    #img2latlongrd world_grav.img.7.2 -R-149/-135/52.5/58 -Gss_grav.grd -T1

    # Use spherical projection since SS data define on sphere
    gmtset ELLIPSOID Sphere

    # Define location of Pratt seamount
    echo "-142.65    56.25" >! pratt.d

    # First generate gravity image w/ shading, label Pratt, and draw a circle
    # of radius = 200 km centered on Pratt.

    makecpt -Crainbow -T-60/60/10 -Z >! grav.cpt
    grdgradient ss_grav.grd -Nt1 -A45 -Gss_grav_i.grd
    grdimage ss_grav.grd -Iss_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

    # Then draw 10 mGal contours and overlay 50 mGal contour in green

    grdcontour ss_grav.grd -JM -C20 -B2f1WSEn -O -K -Y-4.85i -U/-1.25i/-0.75i/"Example 18 in Cookbook" \
       >> example_18.ps
    grdcontour ss_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

    # Now determine centers of each enclosed seamount > 50 mGal but only plot
    # the ones within 200 km of Pratt seamount.

    # First make a simple $AWK script that returns the average position
    # of a file with coordinates x, y (remember to escape the $ sign)

    cat << EOF >! center.awk
    BEGIN {
        x = 0
        y = 0
        n = 0
    }
    {
        x += \$1
        y += \$2
        n++
    }
    END {
        print x/n, y/n
    }
    EOF

    # Now determine mean location of each closed contour and
    # add it to the file centers.d

    \rm -f centers.d
    foreach file (sm_*.xyz)
        $AWK -f center.awk $file >>! centers.d
    end

    # Only plot the ones within 200 km

    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

    # Then report the volume and area of these seamounts only
    # by masking out data outside the 200 km-radius circle
    # and then evaluate area/volume for the 50 mGal contour

    grdmath -R -I2m -F -142.65 56.25 GDIST = mask.grd
    grdclip mask.grd -Sa200/NaN -Sb200/1 -Gmask.grd
    grdmath ss_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

    # Clean up

    \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

Next: Color patterns on maps Up: Cook-book Previous: Images clipped by coastlines