Organizing multiple subplots

%%capture
%load_ext autoreload
%autoreload 2
import pygmt

from polartoolkit import fetch, maps, regions, utils

define a region to plot

plot_region = regions.ross_ice_shelf

create a dictionary of grids names

%%capture

ice_thickness = fetch.bedmachine(region=plot_region, layer="thickness")
GHF = fetch.ghf(region=plot_region, version="burton-johnson-2020")
crustal_thickness = fetch.crustal_thickness(region=plot_region, version="an-2015")
FA_grav = fetch.gravity(region=plot_region, version="antgg-update", anomaly_type="FA")
BA_grav = fetch.gravity(region=plot_region, version="antgg-update", anomaly_type="BA")
basement = fetch.basement(region=plot_region)

grids = [
    ice_thickness,
    GHF,
    crustal_thickness,
    FA_grav,
    BA_grav,
    basement,
]

Warning 1: The definition of projected CRS EPSG:3031 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

Plot the grids

By default, maps.subplots will try to create a square of subplots from the supplied list of grids. You can also set parameter dims=(number of rows, number of columns) to customize this. Below is the resulting subplots with the default settings:

fig = maps.subplots(grids)
fig.show()

Here are some of the keyword arguments you can provide to customize the plots:

fig = maps.subplots(
    grids,
    plot_region,
    margins="2c",
    coast=True,
    grd2cpt=True,
    fig_title="Plotting multiple grids",
    subplot_titles=[
        "Ice thickness",
        "GHF",
        "Crustal thickness",
        "Free-air gravity",
        "Bouguer gravity",
        "Basement elevation",
    ],
    cmaps=[
        "cool",
        "solar",
        "thermal",
        "balance+h0",
        "balance+h0",
        "batlow",
    ],
    cbar_labels=[
        "BedMachine v1",
        "Burton-Johnson et al. 2020",
        "An et al. 2015",
        "Forsberg et al. 2020",
        "Forsberg et al. 2020",
        "Tankersley et al. 2022",
    ],
    cbar_units=["m", "mW/m@+2@", "m", "mGal", "mGal", "m"],
    autolabel=True,
)


fig.show()

grd2cpt [ERROR]: Making a continuous cpt from a discrete cpt may give unexpected results!
PSL: Warning: Super-scripting not terminated [mW/m@+2@]
PSL: Warning: Super-scripting not terminated [mW/m@+2@]
grd2cpt [ERROR]: Making a continuous cpt from a discrete cpt may give unexpected results!
grd2cpt [ERROR]: Making a continuous cpt from a discrete cpt may give unexpected results!
grd2cpt [ERROR]: Making a continuous cpt from a discrete cpt may give unexpected results!

If you want to further customize the plots, it might be better to use the built-in PyGMT subplot function, and set origin_shift='no_shift'.

To help organize titles and grids, we can add everything to a dictionary.

Here we add titles, a shapefile, and a legend to each subplot.

layers = {
    "Ice thickness": {
        "grid": ice_thickness,
    },
    "GHF": {"grid": GHF},
    "Crustal thickness": {"grid": crustal_thickness},
    "Free-air gravity": {"grid": FA_grav},
    "Bouguer gravity": {"grid": BA_grav},
    "Basement elevation": {"grid": basement},
}

subplot_dimensions = utils.square_subplots(len(layers.items()))

proj, proj_latlon, fig_width, fig_height = utils.set_proj(plot_region, 15)

fig = pygmt.Figure()

# fetch GHF data
df = fetch.ghf(version="burton-johnson-2020", points=True)

with fig.subplot(
    nrows=subplot_dimensions[0],
    ncols=subplot_dimensions[1],
    subsize=(fig_width, fig_height),
    frame="f",
    margins="1.5c",
):
    for i, (k, v) in enumerate(layers.items()):
        with fig.set_panel(panel=i):
            # plot the grids
            maps.plot_grd(
                v["grid"],
                fig=fig,
                title=k,
                origin_shift="no_shift",
                region=plot_region,
                coast=True,
            )
            # plot GHF point measurements
            # make a colorscale
            pygmt.makecpt(cmap="inferno", series=[30, 90])
            # plot the points
            fig.plot(
                x=df.x,
                y=df.y,
                fill=df.GHF,
                cmap=True,
                pen=".8p,black",
                style="c.4c",
            )
            # add a colorbar
            fig.colorbar(
                position="JMR+o1c/0c+w6c/.6c", frame=["x+lGHF", "y+lmW/m@+2@+"]
            )

fig.show()

Alternatively, if you create your own layout by using the origin_shift parameter of maps.plot_grid.

Use kwargs xshift_amount and yshift_amount to give number of figs to shift origin by. (e.g. xshift_amount=3 would be 3 figure widths)

# get dimensions (2 rows x 3 columns)
subplot_dimensions = utils.square_subplots(len(layers.items()))

# iterate through the grids in the dictionary
for i, (k, v) in enumerate(layers.items()):
    # for the first grid, need orgin_shift to be default
    if i == 0:
        fig = maps.plot_grd(
            v["grid"],
            title=k,
            region=plot_region,
            coast=True,
        )
        # plot a shapefile
        fig.plot(
            data=fetch.sample_shp(name="Roosevelt_Island"),
            pen="2p,magenta",
            label="Roosevelt Island Transect",
        )
        # add legend
        fig.legend(position="JTR+jTR", box="+gwhite+p1p")

    # for the end of each row (mulitples of the number of columns) need to
    # shift origin down and back to beginning
    elif i % subplot_dimensions[1] == 0:
        fig = maps.plot_grd(
            v["grid"],
            title=k,
            fig=fig,
            origin_shift="both_shift",
            xshift_amount=-(subplot_dimensions[0]),  # gives -2, shifts plot back by 2
            yshift_amount=-1,  # shift new row down by 1
            region=plot_region,
            coast=True,
        )
        # plot a shapefile
        fig.plot(
            data=fetch.sample_shp(name="Roosevelt_Island"),
            pen="2p,magenta",
            label="Roosevelt Island Transect",
        )
        # add legend
        fig.legend(position="JTR+jTR", box="+gwhite+p1p")

    # for the rest of the grids, just shift to the right like normal
    else:
        fig = maps.plot_grd(
            v["grid"],
            title=k,
            fig=fig,
            origin_shift="xshift",
            region=plot_region,
            coast=True,
        )
        # plot a shapefile
        fig.plot(
            data=fetch.sample_shp(name="Roosevelt_Island"),
            pen="2p,magenta",
            label="Roosevelt Island Transect",
        )
        # add legend
        fig.legend(position="JTR+jTR", box="+gwhite+p1p")

fig.show()