How to plot earthquakes data on a three-dimensional topographic map

Read the earthquake data from a data file and overlay on a three-dimensional topographic map using PyGMT

Continuing our previous post on overlaying the earthquake location data on topographic map, this time we will overlay the earthquake data on a three-dimensional map. The three-dimensional map will be very similar to the one we built using PyGMT in one of our previous post.

Read the data file

Taiwan is located at a compressive tectonic boundary between the Eurasian Plate and the Philippine Sea Plate. The island straddles the compressional plate boundary between the Eurasian Plate to the west and the north and the Philippine Sea Plate to the east and south. The present convergence rate is about 7 cm/year in a NW-SE direction.

The data for this post was downloaded from Central Weather Bureau of Taiwan, and saved as evtlist.dat. The dataset contains 2317 events with the date range from 1995-07-05 to 2019-12-16. The depth of the events ranges from 11-180 km, and magnitude (Mw) from 3.0-7.1.

We can read the tabular data using pandas’ read_csv method.

import pandas as pd
# dataformat: 1 1995-07-05 17:33:48.600 24.8313 122.1158 27 4.87
data = pd.read_csv('evtlist.dat', names=['SN','year','time','latitude','longitude','depth_km','magnitude'], delimiter='\s+')

Load the topographic data

We load the topographic data directly using the PyGMT load_earth_relief method. The PyGMT will download the specified resolution of the topographic data locally.

minlon, maxlon = 119.5, 124.0
minlat, maxlat = 19.8, 26.0
# Load sample earth relief data
grid = pygmt.datasets.load_earth_relief(resolution="03s", region=[minlon, maxlon, minlat, maxlat])

Define three-dimensional perspective

We first specify the cpt files for the topographic colors, and then we specify the three-dimensional perspective view. The view is set to have the azimuth and elevation angle of the viewpoint to be 150 and 30 degrees. The cover type of the grid is taken as “surface”. The other options are mesh, waterfall, image, etc. See the pygmt.Figure.grdview for details.

We define the region by the bounds of latitudes, longitudes and depths. The projection was taken to be “Mercator” with width of 15cm. The vertical height was taken to be 4cm.

pygmt.makecpt(
        cmap='geo',
        series=f'-{maxdep}/4000/100',
        continuous=True
    )
fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    region=[minlon, maxlon, minlat, maxlat, -maxdep, 4000],
    perspective=[150, 30],
    frame=frame,
    projection="M15c",
    zsize="4c",
    surftype="i",
    plane=f"-{maxdep}+gazure",
    shading=0,
    # Set the contour pen thickness to "1p"
    contourpen="1p",
)

Plot the data on a topographic map

Since the depths of the event ranges up to 180 km, we had to constrain it to show it within the range of our chosen maximum depth for the three-dimensional map. We project the depths in the range of 0 to the maximum depth of the displayed topography.

Another good choice for displaying the actual depth for each earthquake events is displaying on a cross-sectional map. The advantage of displaying the events on the topographic map is understanding its relation with the tectonics (topography or bathymetry).

depthData = data.depth_km.values
depthData = maxdep*(depthData - data.depth_km.min())/(data.depth_km.max()-data.depth_km.min()) #normalize the depth to the maximum of maxdep

# colorbar colormap
pygmt.makecpt(cmap="jet", series=[
              depthData.min(), depthData.max()])

fig.plot3d(
    x=data.longitude,
    y=data.latitude,
    z=-1*depthData,
    size=0.05*data.magnitude, #scale the magnitude
    color=depthData,
    cmap=True,
    style="cc",
    pen="black",
    perspective=True,
)

Add colorbar to the figure

At last, we can add the colorbar to show the values of actual depths of the event.

pygmt.makecpt(cmap="jet", series=[
              data.depth_km.min(), data.depth_km.max()])
fig.colorbar(frame='af+l"Depth (km)"', perspective=True,)

Three dimensional map of Taiwan with overlayed earthquake data
Three dimensional map of Taiwan with overlayed earthquake data

Complete script

import pygmt
import pandas as pd
# 1995-07-05 17:33:48.600 24.8313 122.1158 27 4.87
data = pd.read_csv('evtlist.dat', names=['SN','year','time','latitude','longitude','depth_km','magnitude'], delimiter='\s+')


minlon, maxlon = 119.5, 124.0
minlat, maxlat = 19.8, 26.0
# Load sample earth relief data
grid = pygmt.datasets.load_earth_relief(resolution="03s", region=[minlon, maxlon, minlat, maxlat])

maxdep = 10000
frame =  ["xa1f0.25","ya1f0.25", "z2000+lmeters", "wSEnZ"]

pygmt.makecpt(
        cmap='geo',
        series=f'-{maxdep}/4000/100',
        continuous=True
    )
fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    region=[minlon, maxlon, minlat, maxlat, -maxdep, 4000],
    perspective=[150, 30],
    frame=frame,
    projection="M15c",
    zsize="4c",
    surftype="i",
    plane=f"-{maxdep}+gazure",
    shading=0,
    # Set the contour pen thickness to "1p"
    contourpen="1p",
)


depthData = data.depth_km.values
depthData = maxdep*(depthData - data.depth_km.min())/(data.depth_km.max()-data.depth_km.min())

# colorbar colormap
pygmt.makecpt(cmap="jet", series=[
              depthData.min(), depthData.max()])
# pygmt.makecpt(cmap="jet", series=[
#               data.depth_km.min(), data.depth_km.max()])
fig.plot3d(
    x=data.longitude,
    y=data.latitude,
    z=-1*depthData,
    size=0.05*data.magnitude,
    color=depthData,
    cmap=True,
    style="cc",
    pen="black",
    perspective=True,
)
pygmt.makecpt(cmap="jet", series=[
              data.depth_km.min(), data.depth_km.max()])
fig.colorbar(frame='af+l"Depth (km)"', perspective=True,)

fig.savefig("topo-plot_3d.png", crop=True, dpi=300)

References

  1. Yi-Hsuan Wu, Chien-Chih Chen, Donald L. Turcotte, John B. Rundle, Quantifying the seismicity on Taiwan, Geophysical Journal International, Volume 194, Issue 1, July 2013, Pages 465–469, https://doi.org/10.1093/gji/ggt101
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: 32

Leave a Reply

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