Tutorial 1 - First figure 🚀 and Subplots / layout

Tutorial 1 - First figure 🚀 and Subplots / layout#

This tutorial will cover the fundamental concepts behind making figures with PyGMT:

  • Drawing coastlines

  • Drawing a map frame

  • Choosing a projection

  • Downloading remote dataset

  • Imaging Earth’s relief

  • Creating colorbar

  • Subplot in a figure

Note

This tutorial is part of the AGU24 annual meeting GMT/PyGMT pre-conference workshop (PREWS9) Mastering Geospatial Visualizations with GMT/PyGMT

History

References

Fee free to play around with these code examples 🚀. In case you found any kind of error, just report it by opening an issue or provide a fix via a pull request. Please use the GMT forum to ask questions.

Tip

In Jupyter, some shortcut key will increase your learning curve!

  • Run cell: Shift+Enter

  • Code indentation: TAB

  • Code dedenting: Shift+TAB

  • Writing single-line comments: beginning with # or Command+/

  • Auto-completion: TAB while you writing some code

0️⃣ Importing#

First thing to do is load PyGMT (import) so that we can access its functionality.

PyGMT has a flat package layout, meaning that you can access everything in it with a single import.

import pygmt

1️⃣ Starting your first figure – pygmt.Figure#

Every figure must start with the creation of a pygmt.Figure()

fig = pygmt.Figure()

We have blank canvas in the fig variable now, we will overlay elements and show it on same fig.

1.1 Drawing coastline – pygmt.Figure.coast#

Before plotting, you need to define

  1. region: Controlling the geographic or Cartesian extent of the figure. Regions are specified as lists of [xmin, xmax, ymin, ymax]. Another shortcut, region="g", which specifies a global domain.

  2. projection: Controlling the map projection. GMT (and therefore PyGMT) supports 31 different map projections, from basic Cartesian axes to arbitrary vertical perspectives.

Determining the size of your images 10c, the size will affect how large of your following font size and marker size

length unit c: centimeter (default)

region = [-180, -150, 50, 65]

fig.coast(region=region, projection="M10c", shorelines=True)

# To see the figure looks like
fig.show()
_images/13f0b994a9ed46393bf80406f01993b274dd37a147cee7ff87dc3557564e2586.png

See also

On Jupyter, show will embed a PNG of the figure directly into the notebook. But it can also open a PDF in an external viewer, which is probably what you want if you’re using a plain Python script. See the documentation for pygmt.Figure.show for more information.

See also

If the image size is too small, the default annotation fonts may appear disproportionately large, impacting the overall visual balance and readability.

1.2 Coloring and add map element (frame and ticks)#

  1. axis labels WSne: If an upper-case letter (W, S, N, E) is passed, the axis is plotted with tick marks and annotations. The lower-case version (w, s, n, e) plots the axis only with tick marks. To only plot the axis pass (l, b, t, r).

  2. annotations a: tick labels

  3. frame f: ticks

  4. grid g: grid line in the figure.

frame 
fig.coast(water="lightblue", land="grey", frame=["WSne", "a10f5g5"])
fig.show()
_images/863cc90e60d69e0edf6f3996bbb6dcd12786ffbaa18c6713ac3b2268855236c3.png

1.3 Saving a figure to a file – pygmt.Figure.savefig#

To save the figure to a file you can use the pygmt.Figure.savefig method, and pass the desired file name and extension to the fname parameter.

fig.savefig(fname="my_first_figure.png")

1.4 Stacking approach of GMT / PyGMT#

  1. In GMT/PyGMT, plotting is achieved by layering new elements, meaning that each new element is stacked on top of the previous layers. Therefore, if you draw a black line in an earlier layer and then add a new layer (such as color filling), these new layers might cover the original black line, making it invisible.

  2. In a same figure, once you define region/projection before, you don’t need to define again.

2️⃣ Downloading global dataset – pygmt.datasets#

Before you access remote dataset from PyGMT, you need to define

  1. region: Region of interest, format is [xmin, xmax, ymin, ymax]

  2. resolution: Grid resolution. The suffix d (arc-degrees), m (arc-minutes), and s (arc-seconds).

grid = pygmt.datasets.load_earth_relief(resolution="03m", region=region)

# You also can access then cut the grid
# and this way can save grid (outgrid) as a file in your folder.
# https://www.generic-mapping-tools.org/remote-datasets/

# pygmt.grdcut(
#    grid="@earth_relief_03m", region=region, outgrid="Alaska.grd"
# )

grdblend [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

grdblend [NOTICE]: SRTM15 Earth Relief v2.6 at 03x03 arc minutes reduced by Gaussian Cartesian filtering (15.7 km fullwidth) [Tozer et al., 2019].

grdblend [NOTICE]:   -> Download 90x90 degree grid tile (earth_relief_03m_g): N00W180

A way to confirm your grid information pygmt.grdinfo

print(pygmt.grdinfo(grid=grid))

: Title: 
: Command: 
: Remark: 
: Gridline node registration used [Geographic grid]
: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
: x_min: -180 x_max: -150 x_inc: 0.05 (3 min) name: x n_columns: 601
: y_min: 50 y_max: 65 y_inc: 0.05 (3 min) name: y n_rows: 301
: v_min: -7521.5 v_max: 4634 name: z
: scale_factor: 1 add_offset: 0
: format: classic
: Default CPT:

3️⃣ Visualizing grid – pygmt.Figure.grdimage#

Project and plot grids or images

Now, we start with a new canvas. So you need to set up region/projection again!

Before imaging, you need to define

  1. grid: Specifies the grid file (e.g., a topographic grid) from which data will be sampled.

fig = pygmt.Figure()

fig.grdimage(grid=grid, region=region, projection="M10c")

fig.coast(shorelines="1p,black", frame=["WSne", "a5f2.5"])

fig.show()
_images/ea6e804893275306c0027467e15e56648e4ecf8c68c88befb8596f2177f690d2.png

3.1 Customizing colorbar – pygmt.makecpt#

Making GMT color palette tables (CPT).

You need to define:

  1. cmap: Choosing a colormap to image your grid/dataset.

  1. series: Defining the range of the CPT by giving the lowest and highest z-value. e.g., [0, 500] (without interval) or [0, 500, 10] (defining the interval)

pygmt.makecpt(
    cmap="geo",
    series=[-8000, 7000, 500],
    continuous=True,
    # The file name with extension .cpt to store the generated CPT file
    output="eleva.cpt",
)

Now, re-plot the figure with new colormap and …

3.2 Adding pygmt.Figure.colorbar#

Plot a gray or color scale-bar on maps.

You need to define,

  1. cmap: File name of a CPT file

  2. frame: Setting color bar boundary frame, labels, and axes attributes

    • +l given colorbar title

  3. position: Defining the reference point on the map for the color scale

    Four coordinate systems:

    • g for map (user) coordinates

    • j or J for setting reference point via a 2-character justification code

    • x for plot coordinates: x/y (specific position)

    • n for normalized (0-1) coordinates

    Others

    • +w for length and width of the color bar. e.g., +w10c/0.5c

    • +h horizontal colorbar, if you want to plot vertical one, use +v

    • +m for colorbar title position

position
fig = pygmt.Figure()

fig.grdimage(grid=grid, region=region, projection="M10c", cmap="eleva.cpt")

fig.coast(shorelines="1p,black", frame=["WSne", "a5f2.5"])

fig.colorbar(
    cmap="eleva.cpt",
    frame="a2000f500+lElevation (m)",
    position="JBC+w10c/0.4c+mc+h",
)

fig.show()
_images/d9e3a5b11d774ff087afeb95b8bfd3e7974afe8b1e8f51d12a4736419cdc6974.png

4️⃣ Creating multi-panel figures in a canvas – pygmt.Figure.subplot and pygmt.Figure.set_panel#

You need to define the figure layout with figure.subplot:

  1. nrows: Number of vertical rows

  2. ncols: Number of horizontal columns

  3. figsize or subsize: Entire figure dimensions, e.g., [width, height]

Second, setting each plot with figure.set_panel:

  1. panel: Setting the current subplot panel to plot on. e.g., index or [row, col]

subplot
fig = pygmt.Figure()

grid_crust_age = pygmt.datasets.load_earth_age(
    resolution="30m", region=[-210, -150, 35, 65]
)

with fig.subplot(nrows=1, ncols=2, figsize=["18c", "12c"]):

    with fig.set_panel(panel=0):

        fig.grdimage(
            grid=grid,
            region=region,
            # "?" means that map width automatically determined from the subplot width.
            projection="M?",
            cmap="eleva.cpt",
        )
        fig.coast(shorelines="1p,black", frame=["WSne", "a5f2.5"])
        fig.colorbar(
            cmap="eleva.cpt",
            position="JBC+w8c/0.3c+mc+h",
            frame="a2000f500+lElevation (m)",
        )

    with fig.set_panel(panel=1):

        pygmt.makecpt(cmap="@earth_age.cpt", series=[0, 340])
        fig.grdimage(
            grid=grid_crust_age,
            region=[-210, -150, 35, 65],
            # changing projection method
            projection="L-180/55/35/65/?",
            frame=["a5f2.5", "SEnw"],
            cmap=True,
        )

        fig.colorbar(
            position="JBC+w6c/0.3c+mc+h", frame="a40f20+lSeafloor crustal age (Ma)"
        )

fig.show()

gmtread [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

gmtread [NOTICE]: EarthByte Earth Seafloor Age at 30x30 arc minutes reduced by Gaussian Cartesian filtering (157.3 km fullwidth) [Seton et al., 2020].

gmtread [NOTICE]:   -> Download grid file [205K]: earth_age_30m_g.grd
_images/595c695c25a21e4711d8c268a368ef6fee13a6d9d7918e1ebdb92d8b0bbd4475.png

5️⃣ Basic projection types#

position By default, PyGMT will use an equidistant cylindrical projection if the region seems to be geographic longitude and latitude.

  1. Mercator Projection – mapping the Earth onto a cylinder, preserving angles, making it useful for navigational purposes.

  • Use Case: Ideal for world maps where direction needs to be preserved, such as marine navigation. This projection distorts size, especially near the poles.

  • GMT Command: M[lon]/[lat], specifying the central longitude and latitude.

  1. Conic Projection – projecting the Earth onto a cone.

  • Use Case: for regional maps, especially in mid-latitude countries, where distortion of shapes and areas is minimized within specific latitude ranges.

  • GMT Command: L[lon0]/[lat0]/[lat1]/[lat2], where lon0 and lat0 define the projection center and lat1/lat2 define the standard parallels.

  1. Azimuthal Projection – projecting the Earth onto a plane, preserving directions from a central point.

  • Use Case: Suitable for mapping polar regions or any area where direction from a central point is essential.

  • GMT Command: E[lon]/[lat], with lon and lat defining the central point.

Contents