Seafloor age contours in pygmt

Q&A PyGMT Q&A
1

Hello,

I am new to pygmt and trying to create an age contour map similar to Example 49 in the GMT gallery (https://docs.generic-mapping-tools.org/latest/gallery/ex49.html). My question is about makecpt and how to get the contours from ‘@earth_age_03m’ I don’t understand how to write makecpt for the color code contours and how to make the contours in pygmt.

Thank you!

I’ve pasted my script below (.txt files are lat longs to stations):

import pygmt

GISBORNE_COORDS = (178, -38.65)

KWARGS = dict(
    grid='@earth_relief_15s',
    region=[175,230,-45,-15],
    projection= 'M4i',
    frame= True,)

KWARGS = dict(
    grid='@earth_age_03m',
    region=[175,230,-45,-15],
    projection= 'M4i',
    frame= True,)

fig = pygmt.Figure() 
fig.grdimage(**KWARGS)
fig.grdimage(shading=True, **KWARGS)  # Add illumination!
with fig.inset(position = "jTL+w2.0c+o0.05c", margin = 0, box = "+pred"):
    fig.coast(region="g", projection = "G180/-35/2.0c",land = "gray", water = "white",dcw = "NZ+gred")

inventory = pd.read_csv("Good23ORCAdb.txt", sep='\s+')

pygmt.makecpt()

POINT_FILL = 'blue'
fig.plot(
    x = inventory.lon,
    y = inventory.lat,
    style='c0.05i',
    color=POINT_FILL,
    pen='black',
    label=f'{POINT_FILL}',)

fig.text(x=175, y=-37,position=None,text='Gisborne', angle=0, pen='black',
        font='7p,Helvetica-Bold,black', justify='LM')
fig.plot(*GISBORNE_COORDS, style='s0.15i', color='red', pen=True, label='Gisborne')
fig.show()

inventory = pd.read_csv("Bad7ORCAdb.txt", sep='\s+')
POINT_FILL = 'red'
fig.plot(
    x = inventory.lon,
    y = inventory.lat,
    style='c0.05i',
    color=POINT_FILL,
    pen='black',
    label=f'{POINT_FILL}',)

fig.show(method='external')
2

Here’s an example for plotting the relief and colored age contours in the bottom panel from example 49:

import pygmt
fig = pygmt.Figure()
fig.grdimage(grid="@earth_relief_02m",region=[-30,5,-30,-5],projection= "M4i",frame=True)
pygmt.makecpt(cmap="hot", series=[0, 100, 10], output="t.cpt")
fig.grdcontour(grid="@earth_age_02m",interval="t.cpt",annotation="+f14p",pen="a0.1p+c",label_placement="L30W/22S/5E/13S")
fig.show()
3

Kia Ora @joey! Just to expand on @maxrjones’s answer, please make sure that pygmt.makecpt() is called after the figure is initialized with fig = pygmt.Figure() or you’ll encounter problems with colourmaps getting mixed around as in https://github.com/GenericMappingTools/pygmt/issues/372.