At the end of the article I posted earlier (I wrote a beach ball with Matplotlib + ObsPy), the cooperation between ObsPy and the map library Cartopy has not been confirmed. I wrote that, but it turned out that ** it is possible to draw on a map **, so I will describe this.

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 install Cartopy](#4-Cartopy のインストール方法) 5. Demonstration 6. Conclusion 7. References

2. Purpose of this page

Draw the source sphere of the lower hemisphere projection on the map.

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

cartopy was installed with pip.

4. How to install Cartopy

According to Official description, assuming the environment of Anaconda, it is quick to enter with the conda command conda install -c conda-forge cartopy. Also, if required dependencies is satisfied, you can install it with pip install cartopy. Probably, if you are using Anaconda, installing with pip should automatically resolve this required dependency.

From here on, it's a shame, but if you installed python "directly" on windows 10 via the official installer or store, installing Cartopy can be quite a daunting task. When I tried before switching to Anaconda, when I tried to pip cartopy, it said that the version of proj was low, and when I tried to pip proj, I didn't have the required version (proj4). No, when I try to install proj4, I end up installing OSGeo4w, so it's okay to install OSGeo4w, but how do I use it? I fell into the situation and gave up. Please let me know if there is an article that even I, a beginner, can understand ...

By the way, in Ubuntu, I entered without any problem even if it was not Anaconda. Ubuntu is a good boy.

5. Demonstration

So let's run a suitable script. In My previous work, I found that ObsPy returned a beach ball as a collection of his Matplotlib, so if I knew that I could draw it on a map, I could write it easily.

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

Put this in the same directory as the script below and run it. Regarding drawing a map with Cartopy, please refer to Official script example and This person's blog (I will post it because it was link free ...). I did. If you want to devise a map drawing, please see there.

`cartopy_obspy.py`



import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
from obspy.imaging.beachball import beach

#Map and grid drawing
fig = plt.Figure(figsize=(5,4))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # platecarree (Equirectangular projection)Display with
ax.set_extent((134, 137, 33.5, 36.5)) #Area limitation
grids = ax.gridlines(draw_labels=True) #Draw parallels and meridians with labels
grids.xlocator = mticker.FixedLocator(np.arange(120,150,0.5)) # 0.5 degrees each
grids.ylocator = mticker.FixedLocator(np.arange(20,50,0.5)) 
ax.coastlines(resolution='10m') #Draw coastline with a resolution of 10 m
ax.add_feature(cfeature.LAND, color="lightgrey") #Make the land gray

mfile = pd.read_csv("temp.csv") #Data file
mech = mfile[["strike1", "dip1", "rake1"]] #Focal mechanism solution
loc = mfile[["longitude", "latitude"]] #position
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

for i in range(len(loc)) :
    beach1 = beach(mech.iloc[i,:], xy=loc.iloc[i,:], width=0.2,linewidth=0.4,facecolor=pallet[int(typ.iloc[i,0])])
    ax.add_collection(beach1)

fig.savefig("sample.png ")
plt.close(fig)

The result is the same as Opening.

6. Conclusion

It was easy when I tried it. However, at present, this only allows plotting on a map, and we are still looking for a way to draw a source sphere on a vertical section (probably not). In terms of being able to do that, I feel that GMT has an overwhelming degree of freedom. This script may not be useful except in situations where you have to do it in Python.

7. References

obspy.imaging - Plotting routines for ObsPy — ObsPy Documentation (1.2.0)
https://docs.obspy.org/packages/obspy.imaging.html
1. Beachball Plot — ObsPy Documentation (1.2.0)
  https://docs.obspy.org/tutorial/code_snippets/beachball_plot.html
obspy.imaging.beachball.beach — ObsPy Documentation (1.2.0)
https://docs.obspy.org/packages/autoge/obspy.imaging.beachball.beach.html#obspy.imaging.beachball.beach --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 --Python- (basemap) Introduction of Cartopy | Fushizen
https://www.mitsumatado.com/zen/cartopy/

[Python, ObsPy] I drew a beach ball on the map with Cartopy + ObsPy.