PyGMT: High-Resolution Topographic Map in Python (codes included)

A simple tutorial on how to plot high resolution topographic map using GMT tools in Python

Introduction

If you have been working in seismology, then you must have come across Generic Mapping Tools (GMT) software. It is widely used software, not only in seismology but across the Earth, Ocean, and Planetary sciences and beyond. It is a free, open-source software used to generate publication quality maps or illustrations, process data and make animations. Recently, GMT built API (Application Programming Interface) for MATLAB, Julia and Python. In this post, we will explore the Python wrapper library for the GMP API – PyGMT. Using the GMT from Python script allows enormous capabilities.

The API reference for PyGMT can be accessed from here and is strongly recommended. Although PyGMT project is still in completion, there are many functionalities available.

Step-by-step guide for PyGMT using an example

In this post, we will demonstrate the PyGMT implementation by plotting the topographic map of southern India. We will also plot some markers on the map.

Importing Libraries

The first thing I like to do is to import all the necessary libraries for the task. This keeps the code organized.

import pygmt


minlon, maxlon = 60, 95
minlat, maxlat = 0, 25

Define topographic data source

The topographic data can be accessed from various sources. In this snippet below, I listed a couple of sources.

#define etopo data file
# topo_data = 'path_to_local_data_file'
topo_data = '@earth_relief_30s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)
# topo_data = '@earth_relief_15s' #15 arc second global relief (SRTM15+V2.1)
# topo_data = '@earth_relief_03s' #3 arc second global relief (SRTM3S)

Initialize the pyGMT figure

Similar to the matplotlib’s fig = plt.Figure(), PyGMT begins with the creation of Figure instance.

# Visualization
fig = pygmt.Figure()

Define CPT file

# make color pallets
pygmt.makecpt(
    cmap='topo',
    series='-8000/8000/1000',
    continuous=True
)

Plot the high resolution topography from the data source

Now, we provide the topo_data, region and the projection for the figure to plot. The region can also be provided in the form of ISO country code strings, e.g. TW for Taiwan, IN for India, etc. For more ISO codes, check the wikipedia page here. In this example, we used the projection of M4i, which specifies four-inch wide Mercator projection. For more projection options, check here.

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i'
    )

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i',
    frame=True
    )

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i',
    shading=True,
    frame=True
    )

Plot the coastlines/shorelines on the map

Figure.coast can be used to plot continents, shorelines, rivers, and borders on maps. For details, visit pygmt.Figure.coast.

fig.coast(
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i',
    shorelines=True,
    frame=True
    )

Plot the topographic contour lines

We can also plot the topographic contour lines to emphasize the change in topography. Here, I used the contour intervals of 4000 km and only show contours with elevation less than 0km.

fig.grdcontour(
    grid=topo_data,
    interval=4000,
    annotation="4000+f6p",
    limit="-8000/0",
    pen="a0.15p"
    )

Plot data on the topographic map

## Generate fake coordinates in the range for plotting
lons = minlon + np.random.rand(10)*(maxlon-minlon)
lats = minlat + np.random.rand(10)*(maxlat-minlat)


# plot data points
fig.plot(
    x=lons,
    y=lats,
    style='c0.1i',
    color='red',
    pen='black',
    label='something',
    )

We plot the locations by red cicles. You can change the markers to any other marker supported by GMT, for example a0.1i will produce stars of size 0.1 inches.

Plot colorbar for the topography

Default is horizontal colorbar

# Plot colorbar
fig.colorbar(
    frame='+l"Topography"'
    )

For vertical colorbar:

# For vertical colorbar
fig.colorbar(
frame='+l"Topography"',
position="x11.5c/6.6c+w6c+jTC+v"
)

We can define the location of the colorbar using the string x11.5c/6.6c+w6c+jTC+v. +v specifies the vertical colorbar.

Output the figure to a file

Similar to matplotlib, PyGMT shows the figure by

# save figure
fig.show() #fig.show(method='external')

To save figure to png. PyGMT crops the figure by default and has output figure resolution of 300 dpi for png and 720 dpi for pdf. There are several other output formats available as well.

# save figure
fig.savefig("topo-plot.png", crop=True, dpi=300, transparent=True)


fig.savefig("topo-plot.pdf", crop=True, dpi=720)

If you want to save all formats (e.g., pdf, eps, tif) then

allformat = 1

if allformat:
    fig.savefig("topo-plot.pdf", crop=True, dpi=720)
    fig.savefig("topo-plot.eps", crop=True, dpi=300)
    fig.savefig("topo-plot.tif", crop=True, dpi=300, anti_alias=True)

Complete Script

import numpy as np 
import pygmt

np.random.seed(0)

minlon, maxlon = 60, 95
minlat, maxlat = 0, 25

## Generate fake coordinates in the range for plotting
lons = minlon + np.random.rand(10)*(maxlon-minlon)
lats = minlat + np.random.rand(10)*(maxlat-minlat)


#define etopo data file
topo_data = '@earth_relief_30s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)


# Visualization
fig = pygmt.Figure()

# make color pallets
pygmt.makecpt(
    cmap='topo',
    series='-8000/8000/1000',
    continuous=True
)

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M4i',
    shading=True,
    frame=True
    )

# plot coastlines
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M4i', 
    shorelines=True,
    frame=True
    )

# plot topo contour lines
fig.grdcontour(
    grid=topo_data,
    interval=4000,
    annotation="4000+f6p",
    # annotation="1000+f6p",
    limit="-8000/0",
    pen="a0.15p"
    )

# plot data points
fig.plot(
    x=lons,
    y=lats, 
    style='c0.1i', 
    color='red', 
    pen='black', 
    label='something',
    )

## Plot colorbar
# Default is horizontal colorbar
fig.colorbar(
    frame='+l"Topography"'
    )


# save figure as pdf
fig.savefig("topo-plot.pdf", crop=True, dpi=720)

Plot Focal Mechanism on a Map

How would you plot the focal mechanism on a map? This can be simply done in GMT using the command psmeca. But PyGMT do not support psmeca so far. So, we use a workaround. We call the GMT module using the PyGMT’s pygmt.clib.Session class. This we do using the context manager with in Python.

import pygmt
import numpy as np

minlon, maxlon = 70, 100
minlat, maxlat = 0, 35

## Generate fake coordinates in the range for plotting
num_fm = 15
lons = minlon + np.random.rand(num_fm)*(maxlon-minlon)
lats = minlat + np.random.rand(num_fm)*(maxlat-minlat)

strikes = np.random.randint(low = 0, high = 360, size = num_fm)
dips = np.random.randint(low = 0, high = 90, size = num_fm)
rakes = np.random.randint(low = 0, high = 360, size = num_fm)
magnitudes = np.random.randint(low = 5, high = 9, size = num_fm)

#define etopo data file
topo_data = '@earth_relief_30s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)


fig = pygmt.Figure()

# make color pallets
pygmt.makecpt(
    cmap='topo',
    series='-8000/11000/1000',
    continuous=True
)

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region='IN', 
    projection='M4i',
    shading=True,
    frame=True
    )

fig.coast( region='IN',
    projection='M4i',
    frame=True,
    shorelines=True,
    borders=1, #political boundary
)

for lon, lat, st, dp, rk, mg in zip(lons, lats, strikes, dips, rakes, magnitudes):
    with pygmt.helpers.GMTTempFile() as temp_file: 
        with open(temp_file.name, 'w') as f:
            f.write(f'{lon} {lat} 0 {st} {dp} {rk} {mg} 0 0') #moment tensor: lon, lat, depth, strike, dip, rake, magnitude
        with pygmt.clib.Session() as session:
            session.call_module('meca', f'{temp_file.name} -Sa0.2i')

fig.savefig("fm-plot.pdf", crop=True, dpi=720)
Randomly generated focal mechanisms plotted on topographic map of India

Plotting interstation paths between two stations

The following script uses the sel_pair_csv file for the coordinates of the pairs of stations. The sel_pair_csv file is a csv file that is formatted like:

   stn1     stlo1    stla1  stn2     stlo2    stla2
0  A002  121.4669  25.1258  B077  120.7874  23.8272
1  A002  121.4669  25.1258  B117  120.4684  24.1324
2  A002  121.4669  25.1258  B123  120.5516  24.0161
3  A002  121.4669  25.1258  B174  120.8055  24.2269
4  A002  121.4669  25.1258  B183  120.4163  23.8759


df_info = pd.read_csv(sel_pair_csv, names=["stn1", "stlo1", "stla1", "stn2", "stlo2", "stla2"], header=None)
event_fig = os.path.join("." , "interstation_map.png")
plot_map = 1
if plot_map:
    topo_data = "@earth_relief_15s"
    res = "f"
    figtitle="Interstation-Paths"
    frame =  ["a1f0.25", f"WSen+t{figtitle}"]
    colormap = "geo"
    minlon, maxlon, minlat, maxlat = (
        min(df_info["stlo1"].min(),df_info["stlo2"].min()),
        max(df_info["stlo1"].max(),df_info["stlo2"].max()),
        min(df_info["stla1"].min(),df_info["stla2"].min()),
        max(df_info["stla1"].max(),df_info["stla2"].max()),
    )

    dcoord = 0.5
    minlon, maxlon, minlat, maxlat = (
        minlon - dcoord,
        maxlon + dcoord,
        minlat - dcoord,
        maxlat + dcoord,
    )
    fig = pygmt.Figure()
    fig.basemap(region=[minlon, maxlon, minlat, maxlat], projection="M6i", frame=frame)
    fig.grdimage(
        grid=topo_data,
        shading=True,
        cmap=colormap,
    )

    fig.coast(
        frame=frame,
        resolution=res,
        shorelines=["1/0.2p,black", "2/0.05p,gray"],
        borders=1,
    )
    
    
    fig.plot(
        x=df_info["stlo1"].values,
        y=df_info["stla1"].values,
        style="i10p",
        color="blue",
        pen="black",
    )
    fig.plot(
            x=df_info["stlo2"].values,
            y=df_info["stla2"].values,
            style="i10p",
            color="blue",
            pen="black",
        )
    print("Plotting paths now...")
    for stn1, stn2, stlo1, stla1, stlo2, stla2 in zip(df_info["stn1"].values, df_info["stn2"].values, df_info["stlo1"].values,df_info["stla1"].values,df_info["stlo2"].values,df_info["stla2"].values):
        print(f"-->{stn1}-{stn2}")

        fig.plot(
            x=[stlo1,stlo2],
            y=[stla1,stla2],
            pen="0.1p,red,-",
            straight_line=1,
        )

    fig.savefig(event_fig, crop=True, dpi=300)

Event-Station Map

import pygmt
import pandas as pd
def plot_event_moll_map(df_info, event, event_fig, clon=None, colormap='geo', topo_data = "@earth_relief_20m"):
    '''
    Utpal Kumar
    2021, March
    param df_info: pandas dataframe containing the event and station coordinates (type: pandas DataFrame)
    param event: event name (type: str)
    param event_fig: output figure name (type str)
    param clon: central longitude (type float)
    '''
    res = "f"
    if not colormap:
        colormap = "geo"
    
    if not clon:
        clon = df_info["stlo"].mean()
    proj = f"W{clon:.1f}/20c"

    fig = pygmt.Figure()
    fig.basemap(region="g", projection=proj, frame=True)
    fig.grdimage(
        grid=topo_data,
        shading=True,
        cmap=colormap,
    )

    fig.coast(
        resolution=res,
        shorelines=["1/0.2p,black", "2/0.05p,gray"],
        borders=1,
    )
    fig.plot(
        x=df_info["stlo"].values,
        y=df_info["stla"].values,
        style="i2p",
        color="blue",
        pen="black",
        label="Station",
    )
    fig.plot(
        x=df_info["evlo"].values[0],
        y=df_info["evla"].values[0],
        style="a15p",
        color="red",
        pen="black",
        label="Event",
    )
    for stlo, stla in zip(df_info["stlo"].values, df_info["stla"].values):
        fig.plot(
            x=[df_info["evlo"].values[0], stlo],
            y=[df_info["evla"].values[0], stla],
            pen="red",
            straight_line=1,
        )

    fig.savefig(event_fig, crop=True, dpi=300)

if __name__=="__main__":
    event="test_event"
    data_info_file = f"data_info_{event}.txt"
    event_fig = f"event_map_{event}.png"

    df_info = pd.read_csv(data_info_file)
    plot_event_moll_map(df_info, event, event_fig)

Here, I used the Mollweide projection but this code can be applied for any other projections with a little tweak.

The first few lines of the data file are:

network,station,channel,stla,stlo,stel,evla,evlo,evdp,starttime,endtime,samplingRate,dist,baz
TW,GWUB,HHZ,24.5059,121.1131,2159.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000002Z,2018-12-11T03:56:31.990002Z,100.0,15445.41,205.14
TW,LXIB,HHZ,24.0211,121.4133,1327.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,15409.58,204.75
TW,VCHM,HHZ,23.2087,119.4295,60.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,15242.27,205.4
TW,VDOS,HHZ,20.701,116.7306,5.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,14871.63,205.6
TW,VNAS,HHZ,10.3774,114.365,2.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,13726.26,203.23
TW,VWDT,HHZ,23.7537,121.1412,2578.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,15371.08,204.77
TW,VWUC,HHZ,24.9911,119.4492,42.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,15420.85,206.25
TW,YD07,HHZ,25.1756,121.62,442.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.000000Z,2018-12-11T03:56:31.990000Z,100.0,15534.34,205.2
TW,HOPB,HHZ,24.3328,121.6929,160.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.003130Z,2018-12-11T03:56:31.993130Z,100.0,15452.83,204.75
TW,SYNB,HHZ,23.2482,120.9862,2352.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.003131Z,2018-12-11T03:56:31.993131Z,100.0,15313.6,204.62
TW,DYSB,HHZ,24.8208,121.4049,1000.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.008393Z,2018-12-11T03:56:31.998393Z,100.0,15489.54,205.14
TW,FUSB,HHZ,24.7597,121.5875,690.0,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.008394Z,2018-12-11T03:56:31.998394Z,100.0,15491.23,205.01
TW,HGSD,HHZ,23.4921,121.4239,134.8,-58.5981,-26.4656,164.66,2018-12-11T02:25:32.008393Z,2018-12-11T03:56:31.998393Z,100.0,15356.77,204.5
---

References

  1. Uieda, L., Wessel, P., 2019. PyGMT: Accessing the Generic Mapping Tools from Python. AGUFM 2019, NS21B–0813.
  2. Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F., Tian, D., 2019. The Generic Mapping Tools Version 6. Geochemistry, Geophys. Geosystems 20, 5556–5564. https://doi.org/10.1029/2019GC008515
Views: 216
Tags
Utpal Kumar
Utpal Kumar

Geophysicist | Geodesist | Seismologist | Open-source Developer
I am a geophysicist with a background in computational geophysics, currently working as a postdoctoral researcher at UC Berkeley. My research focuses on seismic data analysis, structural health monitoring, and understanding deep Earth structures. I have had the opportunity to work on diverse projects, from investigating building characteristics using smartphone data to developing 3D models of the Earth's mantle beneath the Yellowstone hotspot.

In addition to my research, I have experience in cloud computing, high-performance computing, and single-board computers, which I have applied in various projects. This includes working with platforms like AWS, Docker, and Kubernetes, as well as supercomputing environments such as STAMPEDE2, ANVIL, Savio and PERLMUTTER (and CORI). My work involves developing innovative solutions for structural health monitoring and advancing real-time seismic response analysis. I am committed to applying these skills to further research in computational seismology and structural health monitoring.

Articles: 29

Related Posts

13 Comments

  2. Hey just wanted to give you a quick heads up. The words in your post seem to be running off the screen in Chrome. I’m not sure if this is a format issue or something to do with browser compatibility but I thought I’d post to let you know. The layout look great though! Hope you get the problem solved soon. Thanks

  4. Hello there! I could have sworn I’ve been to this blog before but after looking at many of the posts I realized it’s new to me. Regardless, I’m certainly pleased I discovered it and I’ll be bookmarking it and checking back often!

  9. Hi there! This is kind of off topic but I need some help from an established blog. Is it hard to set up your own blog? I’m not very techincal but I can figure things out pretty fast. I’m thinking about creating my own but I’m not sure where to start. Do you have any tips or suggestions? Cheers

Leave a Reply

Your email address will not be published. Required fields are marked *