GMT Basics

test.py

#!/usr/bin/env python

import typer
from pathlib import Path
from typing import NamedTuple
import numpy as np

# NOTE: You need both! of these packages
# pip install typer pygmt "git+https://github.com/ucgmsim/pygmt_helper.git" (or conda, mamba, whatever)
import pygmt
from pygmt_helper import plotting


app = typer.Typer(context_settings={"ignore_unknown_options": True})


class PlotRegion(NamedTuple):  # Renamed to avoid conflict with other 'Region' types
    min_lon: float
    max_lon: float
    min_lat: float
    max_lat: float


def default_pygmt_config() -> None:
    pygmt.config(
        # thin black line for frame
        MAP_FRAME_TYPE="plain",
        # decimal degrees
        FORMAT_GEO_MAP="ddd.xx",
        # even thinner frame
        MAP_FRAME_PEN="thinner, black",
        # set the font
        FONT_ANNOT_PRIMARY="10p,Helvetica,black",
    )


def plot_basemap(
    fig: pygmt.Figure, region: str | PlotRegion, title: str, width: int = 17
):
    fig.basemap(
        region=region,
        projection=f"M{width}c",
        frame=["af", "xaf+lLongitude", "yaf+lLatitude", f"+t{title}"],
    )
    # Loads the earth topography with a 3 arcsecond resolution. smaller = longer
    # and bigger downloads
    grid = pygmt.datasets.load_earth_relief(resolution="03s", region=region)
    # Clip the land only.
    grid = pygmt.grdclip(grid, below=[0.1, np.nan])
    # Equivalent to colour map in matplotib
    pygmt.makecpt(
        cmap="gray",
        # series = scale for cmap
        # appropriate for NZ topography data
        # -900m for contrast. If you want lighter land, set this lower.
        series=[-900, 3100, 1],
        # Continuous. Turn off if you want discrete colours instead.
        continuous=True,
        # Reverse the colourmap.
        reverse=True,
    )
    # Actually plot the topography
    # grdimage: Draw an "image"
    fig.grdimage(
        grid=grid,
        projection=f"M{width}c",
        region=region,
        # cmap = True uses the *last* defined cpt, which is the one we made earlier.
        # Can sometimes use a preset e.g.
        # cmap = 'viridis'
        cmap=True,
        # apply colour shading. Mostly useful for topography. Provides a 3D effect.
        # Most of the time you can leave this line out.
        shading=True,
        # Sets all nan (i.e. the water we clipped earlier) to transparent.
        # If you want to plot a different colour for the ocean for example.
        nan_transparent=True,
    )
    # fig.coast could be used instead, but the lakes are bad.
    # We maintain a custom high-quality lake dataset.
    nz_map_data = plotting.NZMapData.load()
    # Plot the lakes and rivers in white
    fig.plot(data=nz_map_data.water_df, fill="white")
    # Could add more details like roads and highways (see the documentation)
    # https://quakecoresoft.canterbury.ac.nz/docs/pygmt_helper.plotting.NZMapData.html
    # fig.plot(data=nz_map_data.road_df, fill="yellow")


@app.command()
def main(lat: float, lon: float, output_path: Path, width: int = 17):
    # any argument to main becomes a command line argument!
    default_pygmt_config()
    fig = pygmt.Figure()
    region = "NZ"
    title = "My First GMT Plot"
    # region = PlotRegion(min_lon=-123, max_lon=-121, min_lat=37, max_lat=38.5)
    # Sets up the figure
    plot_basemap(fig, region, title, width)
    # "Plot" can plot polygons, points, lines, most things except images.

    # fig.plot(
    #     # x is EW
    #     x=lon,
    #     # y is NS
    #     y=lat,
    #     # the above could be xy if using Cartesian projection (the 'X17c' projection).
    #     #
    #     # Style specifies what the point looks like
    #     # Format = <shape><size>c (c = centimetres)
    #     # https://www.pygmt.org/dev/gallery/symbols/basic_symbols.html#sphx-glr-gallery-symbols-basic-symbols-py
    #     # a = star (usually, hypocentre)
    #     # i = downwards triangle (stations)
    #     #
    #     # so a0.7c is a star with a width of 0.7cm.
    #     style="a0.7c",
    #     # The border. Format is <size>p,<colour> where p = pixels
    #     # so 0.9p,black is 0.9 pixels black border.
    #     pen="0.9p,black",
    #     # colour on the inside.
    #     fill="yellow",
    # )
    #
    # Plotting instead focal mechanisms
    # Plots a beachball instead of a point.
    focal_mechanism = {"strike": 318, "dip": 89, "rake": -179, "magnitude": 7.75}
    # scale is the size, this is 0.5cm.
    fig.meca(spec=focal_mechanism, scale="0.5c", longitude=lon, latitude=lat, depth=0)
    fig.savefig(output_path)


# tells python to run code inside the if statement
# when invoked like python test.py (or using the play button in you IDE).
# BUT NOT WHEN I IMPORT THE CODE (so I can reuse the plotting stuff in other files).
if __name__ == "__main__":
    app()