(18) Volumes and Spatial Selections — GMT 5.1.1 documentation

### Navigation

- index
- next
- previous
- Project Home »
- GMT 5.1.1 documentation »
- The GMT historical collection »
- Search:

# (18) Volumes and Spatial Selections

To demonstrate potential usage of the programs gmtspatial, grdvolume and gmtselect we extract a subset of the Sandwell & Smith altimetric gravity field [1] 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 use gmtspatial 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 (gravimetric) volumes of the chosen seamounts. Our illustration is presented in Figure.

    #!/bin/bash
    #		GMT EXAMPLE 18
    #
    # Purpose:	Illustrates volumes of grids inside contours and spatial
    #		selection of data
    # GMT progs:	gmtset, gmtselect, gmtspatial, grdclip, grdcontour, grdgradient, grdimage
    # GMT progs:	grdmath, grdvolume, makecpt, pscoast, psscale, pstext, psxy
    # Unix progs:	$AWK, cat, rm
    #
    ps=example_18.ps

    # Use spherical gmt projection since SS data define on sphere
    gmt gmtset PROJ_ELLIPSOID Sphere FORMAT_FLOAT_OUT %g

    # Define location of Pratt seamount and the 400 km diameter
    echo "-142.65 56.25 400" > pratt.d

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

    gmt makecpt -Crainbow -T-60/60/120 -Z > grav.cpt
    gmt grdgradient AK_gulf_grav.nc -Nt1 -A45 -GAK_gulf_grav_i.nc
    gmt grdimage AK_gulf_grav.nc -IAK_gulf_grav_i.nc -JM5.5i -Cgrav.cpt -B2f1 -P -K -X1.5i \
    	-Y5.85i > $ps
    gmt pscoast -RAK_gulf_grav.nc -J -O -K -Di -Ggray -Wthinnest >> $ps
    gmt psscale -D2.75i/-0.4i/4i/0.15ih -Cgrav.cpt -Bx20f10 -By+l"mGal" -O -K >> $ps
    $AWK '{print $1, $2, "Pratt"}' pratt.d | gmt pstext -R -J -O -K -D0.1i/0.1i \
    	-F+f12p,Helvetica-Bold+jLB >> $ps
    gmt psxy pratt.d -R -J -O -K -SE- -Wthinnest >> $ps

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

    gmt grdcontour AK_gulf_grav.nc -J -C20 -B2f1 -BWSEn -O -K -Y-4.85i >> $ps
    # Save 50 mGal contours to individual files, then plot them
    gmt grdcontour AK_gulf_grav.nc -C10 -L49/51 -Dsm_%d_%c.txt
    gmt psxy -R -J -O -K -Wthin,green sm_*.txt >> $ps
    gmt pscoast -R -J -O -K -Di -Ggray -Wthinnest >> $ps
    gmt psxy pratt.d -R -J -O -K -SE- -Wthinnest >> $ps
    rm -f sm_*_O.txt	# Only consider the closed contours

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

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

    gmt gmtspatial -Q -fg sm_*_C.txt > centers.d

    # Only plot the ones within 200 km

    gmt gmtselect -C200k/pratt.d centers.d -fg | gmt psxy -R -J -O -K -SC0.04i -Gred -Wthinnest >> $ps
    gmt psxy -R -J -O -K -ST0.1i -Gyellow -Wthinnest pratt.d >> $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

    gmt grdmath -R `$AWK '{print $1, $2}' pratt.d` SDIST = mask.nc
    gmt grdclip mask.nc -Sa200/NaN -Sb200/1 -Gmask.nc
    gmt grdmath AK_gulf_grav.nc mask.nc MUL = tmp.nc
    area=`gmt grdvolume tmp.nc -C50 -Sk | cut -f2`
    volume=`gmt grdvolume tmp.nc -C50 -Sk | cut -f3`

    gmt psxy -R -J -A -O -K -L -Wthin -Gwhite >> $ps << END
    -148.5	52.75
    -141	52.75
    -141	53.75
    -148.5	53.75
    END
    gmt pstext -R -J -O -F+f14p,Helvetica-Bold+jLM >> $ps << END
    -148 53.08 Areas: $area km@+2@+
    -148 53.42 Volumes: $volume mGal\264km@+2@+
    END

    # Clean up

    rm -f grav.cpt sm_*.txt *_i.nc tmp.nc mask.nc pratt.d center* gmt.conf

Volumes and geo-spatial selections.

[1] | See http://topex.ucsd.edu/marine_grav/mar_grav.html.

### Navigation

- index
- next
- previous
- Project Home »
- GMT 5.1.1 documentation »
- The GMT historical collection »
- Search: