One of the features of Matplotlib is that it is very easy to draw a 3D graph. Taking advantage of this feature, I tried to draw a 3D scatter plot of the epicenter with the map displayed in the map library Cartopy.

The displayed earthquake data is randomly selected and downloaded from the earthquake mechanism (National Research Institute for Earth Science and Disaster Prevention: Search for Mechanism Solution) published by National Research Institute for Earth Science and Disaster Prevention. The images on this page were drawn by individuals as demonstrations and do not require academic significance in the results.
Please use the script on this page at your own risk. Drawing example

1. Table of contents

1. Table of Contents 2. Purpose of this page 3. Usage environment 4. How to put a map into a 3D graph 5. Demonstration 6. Conclusion 7. References

2. Purpose of this page

After plotting the three-dimensional position of the epicenter, draw a map at a depth of 0 km.

3. Usage environment

Windows 10
conda 4.9.2
python 3.8.5
matplotlib 3.3.2
obspy 1.2.2
cartopy 0.18.0

4. How to put a map into a 3D graph

Now, Cartopy doesn't really support 3D plotting in Matplotlib. Therefore, in order to display it on a 3D plot, it is necessary to convert it to some format supported by Matplotlib. Here, we will focus on the format called Collection.

There is precedent in stackoverflow for the basic method. User pelson answers the question "contourf in 3D Cartopy" of user mtb-za, and the script on this page also refers to this. However, as pointed out on the same site, pelson's answer gives an error. This time around, we're outsmarting an ad hoc solution.

The script that pelson answered basically has the following flow.

--Specify the feature to draw from Cartopy --Get only geometry from feature --Convert geometry to the projection you want to draw (here equirectangular projection) --Get the drawing range from the axes of the matplotlib you are drawing and apply the method of intersection between this and the intersecting geometry. --Convert geometry to path --Convert path to polygon --Convert polygon to collection

However, when I make an intersection, I get an error that some geometry is invalid. Imagine that some of the line segments that make up geometry intersect each other. This doesn't work for the script, so I've inserted code before the intersection to eliminate the invalid geometry. It looks like the result is not a problem (I haven't tracked what happens in other areas and projections).

5. Demonstration

The beginning of the data used is as follows.

`temp.csv`


longitude	latitude	depth	strike1	dip1	rake1	strike2	dip2	rake2	mantissa	exponent	type
136.3025	36.2495	11	63	59	94	236	31	84	2.66	21	1
135.5475	35.2028	20	334	85	11	243	79	175	9.92	21	2
136.72	    35.1692	14	3	70	105	145	25	55	3.43	21	1
134.8153	34.4015	11	293	87	28	201	62	177	4.23	21	2

Of these, the type of fault determined from the longitude, latitude, depth, and direction of the principal stress axis is used for drawing.

`cartopy_3d.py`



import itertools
import pandas as pd
import cartopy.feature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from cartopy.mpl.patch import geos_to_path
from matplotlib.collections import PolyCollection
 
mfile = pd.read_csv("temp.csv") #data
loc = mfile[["longitude", "latitude", "depth"]]
typ = mfile[["type"]] #Is it a normal fault type? 0=Normal fault, 1 =Reverse fault, 2 =Strike-slip fault
pallet = ["blue", "red", "lime"] #Color coded by fault type

pallet_list = [] #List colors
for i in range(len(typ)) :
    pallet_list.append(pallet[int(typ.iloc[i,0])])

#Main plot
fig = plt.figure()
ax3d = fig.add_subplot(111, projection='3d')
#Specifying the area
ax3d.set_xlim(134.0,137.0)
ax3d.set_ylim(33.5,36.5)
ax3d.set_zlim(25,0)
#label
ax3d.set_xlabel("longitude")
ax3d.set_ylabel("latitude")
ax3d.set_zlabel("depth")
#Scatter plot
ax3d.scatter(loc["longitude"],loc["latitude"],loc["depth"],"o", c=pallet_list,s=10,depthshade=0)

#Set axis for map
proj_ax = plt.figure().add_subplot(111, projection=ccrs.PlateCarree())
proj_ax.set_xlim(ax3d.get_xlim())
proj_ax.set_ylim(ax3d.get_ylim())
# a usefull function
concat = lambda iterable: list(itertools.chain.from_iterable(iterable))
#Get geometry
target_projection = proj_ax.projection
feature = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
geoms = feature.geometries()
#Area acquisition
boundary = proj_ax._get_extent_geom()
#Projection method conversion
geoms = [target_projection.project_geometry(geom, feature.crs) for geom in geoms]

#Elimination of invalid geometry
geoms2 = []
for i in range(len(geoms)) :
    if geoms[i].is_valid :
        geoms2.append(geoms[i])

geoms = geoms2

# intersection
geoms = [boundary.intersection(geom) for geom in geoms]
# geometry to path to polygon to collection
paths = concat(geos_to_path(geom) for geom in geoms)
polys = concat(path.to_polygons() for path in paths)
lc = PolyCollection(polys, edgecolor='black', facecolor='green', closed=True, alpha=0.2)

#Insert
ax3d.add_collection3d(lc, zs=0)
ax3d.view_init(elev=30, azim=-120)
plt.close(proj_ax.figure) #If not closed plt.Two come out in show
fig.savefig("sample2.png ")
plt.show() 
#plt.close(fig)

Display result

By the way, if you do not intersect, the map will protrude from the drawing area.

6. Conclusion

What did you think. It's fun to move the 3D plot with the mouse (laughs). Other projection methods are also described on the official website of Cartopy, so please use them.

7. References

--National Research Institute for Earth Science and Disaster Prevention Broadband Seismic Observation Network (F-net) Search for Mechanism Solutions
https://www.fnet.bosai.go.jp/event/search.php?LANG=ja

Introduction — cartopy 0.18.0 documentation
https://scitools.org.uk/cartopy/docs/latest/index.html
Using cartopy with matplotlib — cartopy 0.18.0 documentation
https://scitools.org.uk/cartopy/docs/latest/matplotlib/intro.html
matplotlib - contourf in 3D Cartopy - Stack Overflow
https://stackoverflow.com/questions/48269014/contourf-in-3d-cartopy

[Python] I want to make a 3D scatter plot of the epicenter with Cartopy + Matplotlib!