GMT Example #18 (with PyGMT): Volumes and Spatial Selections

Here is the code I used to re-create GMT gallery example 18 using mostly PyGMT. Note that for modules not wrapped in PyGMT (i.e. grdmath) I just used ‘!gmt’ in Jupyter, which let me use GMT as you would normally. Hopefully this example can be updated as more modules are wrapped!

From GMT ex. 18:
“To demonstrate potential usage of the programs gmtspatial, grdvolume and gmtselect we extract a subset of the Sandwell & Smith altimetric gravity field (see footnote #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 18 which also shows the automatic labeling offered by subplot.”

Footnote 1: See http://topex.ucsd.edu/marine_grav/mar_grav.html.

# Example 18. This more complicated script illustrates various search operations in GMT based on geospatial relationships. 
# We seek to find all large seamounts within 200 km of Pratt seamount. 
# We also want to evaluate the combined volume and area of these seamounts. 

import pygmt
from glob import glob
import pandas as pd
import os

# get local version of GMT remote file 'AK_gulf_grav.nc'
# NOTE: "!" allows you to use GMT modules directly without pygmt
!gmt grdcut @AK_gulf_grav.nc -R-149/-135/52.5/58 -GAK_gulf_grav.nc

fig=pygmt.Figure()

# locally set GMT defaults
with pygmt.config(PROJ_ELLIPSOID = 'Sphere', FORMAT_FLOAT_OUT = '%g'): # Use spherical gmt projection since SS data define on sphere 
    with fig.subplot(nrows=2, ncols=1, subsize='15c/9c', autolabel='+JTL+o0.5c',  
                      region='-149/-135/52.5/58', margins='0.5c/1c',
                      frame="WSne", projecton='M14c'): # Note projecton is missing the i in pygmt.?
        # make 1st subplot
        with fig.set_panel(panel=[0,0]): 
            # specify color ramp, plot gravity data and colorbar
            pygmt.makecpt(cmap='rainbow', series='-60/60')
            fig.grdimage(grid='AK_gulf_grav.nc', shading='+d', projection='M14c')
            fig.colorbar(position='JCB+o0/1c', frame=['xaf','y+l"mGal"'])   
            # plot coastline and land
            fig.coast(resolution='i', land='gray', shorelines='thinnest', projection='M14c')
            # make txt file with Pratt location, plot Pratt seamount label and 400 km diameter circle 
            with open('pratt.txt', 'w') as file:
                file.write('-142.65 56.25 400')
            fig.text(x=-142.65, y=56.25, text='Pratt', font='12p,Helvetica-Bold+jLB', offset='8p', projection='M14c')
            fig.plot(data='pratt.txt', style='E-', pen='thinnest', projection='M14c')
        # make 2nd subplot
        with fig.set_panel(panel=[1,0]):
            # plot gravity contours at 20 mGal intervals, 50 mGal contours as greed
            fig.grdcontour(grid='AK_gulf_grav.nc', interval='20', projection='M14c')
            fig.grdcontour(grid='AK_gulf_grav.nc', interval='10', limit='49/51', pen='thin,green', projection='M14c')
            
            # Save 50 mGal contours to individual files and delete the open contours files (with 'O's)
            !gmt grdcontour AK_gulf_grav.nc -C10 -L49/51 -Dsm_%d_%c.txt
            for f in glob('*O.txt'):
                os.remove(f)
            # plot coastline and land    
            fig.coast(resolution='i', land='gray', shorelines='thinnest', projection='M14c')
            # plot 400km circle
            fig.plot(data='pratt.txt', style='E-', pen='thinnest', projection='M14c')
           
            # use 50 mGal closed contour files to calculate area within each contour and centroid location, output to centers.txt
            # I couldn't figure out how to get sm_*_.txt to work. 
            !gmt spatial -Q -fg sm_0_C.txt sm_1_C.txt sm_2_C.txt sm_3_C.txt sm_4_C.txt sm_5_C.txt sm_6_C.txt sm_7_C.txt sm_8_C.txt sm_9_C.txt sm_10_C.txt sm_11_C.txt sm_12_C.txt sm_13_C.txt sm_14_C.txt sm_15_C.txt sm_16_C.txt sm_17_C.txt sm_18_C.txt sm_19_C.txt sm_20_C.txt sm_21_C.txt sm_22_C.txt sm_23_C.txt sm_24_C.txt sm_25_C.txt sm_26_C.txt  > centers.txt
            # select 50 mGal closed contours within 200 km of Pratt and save to centers_close.txt
            !gmt select -Cpratt.txt+d200k centers.txt -fg > centers_close.txt
            # plot seamounts withing 200km of Pratt and plot Pratt itself
            fig.plot(data='centers_close.txt', style='C0.1c', color='red', pen='thinnest', projection='M14c')
            fig.plot(data='pratt.txt', style='T0.25c', color='yellow', pen='thinnest', projection='M14c')
        
            # make a netCDF file with the same regiona and spacing as AK_gulf_grav.nc with values equal to distance from Pratt seamount
            !gmt grdmath -R-149/-135/52.5/58 -I2m pratt.txt POINT SDIST = mask.nc -fg 
            # change to pixle node registration
            !gmt grdsample mask.nc -T -Gmask.nc
            # set every point >200km away to NaN and every point <200km away to 1
            !gmt grdclip mask.nc -Sa200/NaN -Sb200/1 -Gmask.nc 
            # mask gravity data further than 200k from Pratt
            !gmt grdmath AK_gulf_grav.nc mask.nc MUL = tmp.nc 
            
            # calculate area and volume of seamounts within 200km of Pratt and plot the values
            a = !gmt grdvolume tmp.nc -C50 -Sk -o1
            area = int(float(a[0]))
            v = !gmt grdvolume tmp.nc -C50 -Sk -o2
            volume = int(float(v[0]))
            fig.text(x=-149, y=53.25, text='Gravimetric Volumes: {} mGal*km@+2@+'.format(volume), 
                     fill='white', pen='thin', offset='j0.7c', font='14p,Helvetica-Bold+jLB', clearance='8p', projection='M14c')
            fig.text(x=-149, y=52.5, text='Areas: {} km@+2@+'.format(area), 
                     fill='white', pen='thin', offset='j0.7c', font='14p,Helvetica-Bold+jLB', clearance='8p', projection='M14c')

fig.savefig('Ex_18_figure.pdf')
fig.show()[Ex_18_figure.pdf|attachment](upload://eo2Fh6rNZcGwGxNzhmJBjSffI4b.pdf) (205.0 KB) 

Thanks @weiji14 for the help!

2 Likes