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