Subscribe to our weekly newsletter

A Quick Overview on Geospatial Data Visualization using PyGMT

Utpal Kumar   5 minute read      

We will see how to plot a topographic map, overlay earthquake data on topographic maps, plot focal mechanism solutions and plot tomographic results on a geographic map.

Install libraries

Using python env

python -m venv geoviz
source geoviz/bin/activate
pip install pygmt
  • Try: 
    python -c "import pygmt"

    if there’s no ImportError, then you are good to go.

NOTE: If there’s any pygmt import problem, install GMT separately and link the libgmt.dylib file to the file python is looking for!

  • One way to install GMT is conda install gmt -c conda-forge
ln -s ~/miniconda3/envs/boxgmt/lib/libgmt.dylib ~/miniconda3/envs/geoviz/lib/libgmt.dylib

conda create --name geoviz --channel conda-forge pandas pygmt jupyter notebook

Import Libraries

import pygmt
import pandas as pd
import numpy as np
import xarray as xr
from scipy.interpolate import griddata

Plotting a topographic map

#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)

# define plot geographical range
minlon, maxlon = 60, 95
minlat, maxlat = 0, 25

# 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 continents, shorelines, rivers, and borders
fig.coast(
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i',
    shorelines=True,
    frame=True
    )

# plot the topographic contour lines
fig.grdcontour(
    grid=topo_data,
    interval=4000,
    annotation="4000+f6p",
    limit="-8000/0", #to only display it below 
    pen="a0.15p"
    )

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

# save figure
save_fig = 0
if not save_fig:
    fig.show() 
    #fig.show(method='external') #open with the default pdf reader
else:
    fig.savefig("topo-plot.png", crop=True, dpi=300, transparent=True)
#     fig.savefig("topo-plot.pdf", crop=True, dpi=720)
    print('Figure saved!')
Plotting a topographic map
Plotting a topographic map

Plotting data points on a topographic map

We plot 10 randomly generated coordinates on a topographic map with red circles. More symbol options can be found at the GMT site.

## 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 plot geographical range
minlon, maxlon = 60, 95
minlat, maxlat = 0, 25

# 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 continents, shorelines, rivers, and borders
fig.coast(
    region=[minlon, maxlon, minlat, maxlat],
    projection='M4i',
    shorelines=True,
    frame=True
    )

# plot the topographic contour lines
fig.grdcontour(
    grid=topo_data,
    interval=4000,
    annotation="4000+f6p",
    limit="-8000/0", #to only display it only for the regions below the shorelines
    pen="a0.15p"
    )

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

# save figure
fig.show() 
#fig.show(method='external') #open with the default pdf reader
Plotting data points on a topographic map
Plotting data points on a topographic map

Plotting focal mechanism on a map

Focal mechanisms in Harvard CMT convention

Using strike, dip, rake and magnitude

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)


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.show()
Focal mechanisms in Harvard CMT convention
Focal mechanisms in Harvard CMT convention

Focal mechanisms for Seismic moment tensor (Harvard CMT, with zero trace)

This example script needs formatted GCMT soln. Download the file from here

# Define geographical range
minlon, maxlon = -180, -20
minlat, maxlat = 0, 90


## Read GCMT sol
df_gcmt = pd.read_csv('gcmt_sol2.csv')

# ## Subset the data for the given geographical range (will make the program a bit faster to process)
# df_gcmt = df_gcmt[(df_gcmt['evlon'] >= minlon) & (df_gcmt['evlon'] <= maxlon) \
#                   & (df_gcmt['evlat'] >= minlat) & (df_gcmt['evlat'] <= maxlat)]
df_gcmt = df_gcmt[['evmag', 'evlat', 'evlon', 'evdep', 'm_rr','m_tt', 'm_pp', 'm_rt', 'm_rp', 'm_tp']]
df_gcmt.head()

fig = pygmt.Figure()


fig.basemap(region=[minlon, maxlon, minlat, maxlat], projection="Poly/4i", frame=True)

fig.coast(
        land="lightgrey",
        water="white",
        shorelines="0.1p",
        frame="WSNE",
        resolution='h',
        area_thresh=10000
    )

exponent = 16
factor = 10**exponent

#plotting moment tensor sols 
for irow in range(len(df_gcmt)):
#     print(f"{irow}/{len(df_gcmt)-1}")
    m_rr = float(df_gcmt.loc[irow,'m_rr'])/factor
    m_tt = float(df_gcmt.loc[irow,'m_tt'])/factor
    m_pp = float(df_gcmt.loc[irow,'m_pp'])/factor
    m_rt = float(df_gcmt.loc[irow,'m_rt'])/factor
    m_rp = float(df_gcmt.loc[irow,'m_rp'])/factor
    m_tp = float(df_gcmt.loc[irow,'m_tp'])/factor
    evmag = float(df_gcmt.loc[irow,'evmag']) * 0.02
    evdep = float(df_gcmt.loc[irow,'evdep'])
    lat = float(df_gcmt.loc[irow,'evlat'])
    lon = float(df_gcmt.loc[irow,'evlon'])

    # store focal mechanisms parameters in a dict
    focal_mechanism = dict(mrr=m_rr, mtt=m_tt, mff=m_pp, mrt=m_rt, mrf=m_rp, mtf=m_tp, exponent=exponent)
    fig.meca(focal_mechanism, scale=f"{evmag}i", longitude=lon, latitude=lat, depth=evdep,G='blue')

fig.show()
Focal mechanisms for Seismic moment tensor
Focal mechanisms for Seismic moment tensor

Tomographic data on a geographic map

datafile_dcg='dcg_080'

## Read perturbation data
df=pd.read_csv(datafile_dcg,delimiter='\s+', names=['longitude','latitude','pert', 'error'])
lons0=np.array(df['longitude'])
lats0=np.array(df['latitude'])
data=np.array(df['pert'])


coordinates0 = np.column_stack((lons0,lats0))

## Create structured data for plotting
minlon, maxlon = 120., 122.1
minlat, maxlat = 21.8, 25.6
step = 0.01

lons = np.arange(minlon, maxlon, step)
lats = np.arange(minlat, maxlat, step)

## interpolate data on spatial grid
xintrp, yintrp = np.meshgrid(lons, lats)
z1 = griddata(coordinates0, data, (xintrp, yintrp), method='cubic') #cubic interpolation
xintrp = np.array(xintrp, dtype=np.float32)
yintrp = np.array(yintrp, dtype=np.float32)

## xarray dataarray for plotting using pygmt
da = xr.DataArray(z1,dims=("lat", "long"),coords={"long": lons, "lat": lats},)

frame =  ["a1f0.25", "WSen"]

# Visualization
fig = pygmt.Figure()

# make color pallets
lim=abs(max(data.min(),data.max()))
# print(f'{data.min():.2f}/{data.max():.2f}')
pygmt.makecpt(
    cmap='red,white,blue',
    # series=f'{data.min()}/{data.max()}/0.01',
    series=f'-{lim}/{lim}/0.01',
    continuous=True
)

#plot high res topography
fig.grdimage(
    region=[minlon, maxlon, minlat, maxlat],
    grid=da,
    projection='M2i',
    interpolation='l'
    )

# plot coastlines
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    shorelines=True,
    water="#add8e6",
    frame=frame,
    area_thresh=1000
    )


## Plot colorbar
# Default is horizontal colorbar
fig.colorbar(
    frame='+l"Velocity anomaly (%)"'
    )
fig.show()
Tomographic data on a geographic map
Tomographic data on a geographic map

Download notebook and resources

You can download the notebook and the resources from my github repo

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.

Leave a comment