Plot moment tensors in google earth

I usually need to plot the moment tensors on the map, but plot them in the google earth is really cool. This week, I will blog how to do that based on this blog from my friend Thomas. The original code in his blog is still working fine in Python 2, here I run it in Python 3, which need a little change of the code. You can find all the code on Qingkai's Github. You need to install pykml package first if you don't have it. 

Note: pykml is written in Python 2, and if you are using Python 3 as I am doing it here, then you need to add parentheses for the print function in the factory.py file in the pykml. 

!pip install pykml

from lxml import etree      #will manage the identation of the kml script
from pykml.factory import KML_ElementMaker as KML #import pykml library 
import datetime as date
import numpy as np

def beachball(data):
    """function to draw beachball using obspy library"""
    import matplotlib.pyplot as plt
    from obspy.imaging.beachball import beachball
 
    mt =data[:, 9:]
    event=data[:,0]                 #index to identify the beachball created
    for j in range(len(event)):
        #create and save beachball in a outfile in the directory where the .py file is 
        beach = beachball(mt[j, :],outfile=str(int(event[j]))) 

The first thing is that we need to generate a file that contains the list of events with the moment tensors, here the format is: index, yyyy, mm, dd, hr, mn, ss, lat, lon, mrr, mtt, mpp, mrt, mrp, mtp

I generated a file called mt_list, let's see the first 3 lines:

!head -n 3 ./mt_list

# index, yyyy, mm, dd, hr, mn, ss, lat, lon, mrr, mtt, mpp, mrt, mrp, mtp
0.00000,2018.00000,5.00000,29.00000,11.00000,56.00000,11.00000,19.41200,-155.28383,-98650000000000000.00000,53510000000000000.00000,45140000000000000.00000,-75790000000000000.00000,9460000000000000.00000,14910000000000000.00000
1.00000,2018.00000,5.00000,30.00000,20.00000,53.00000,50.00000,19.41083,-155.28583,-107650000000000000.00000,45970000000000000.00000,61690000000000000.00000,-78720000000000000.00000,67750000000000000.00000,22930000000000000.00000

#import your data
# the first row isn't used
data=np.loadtxt('./mt_list',skiprows=1, delimiter=',') 
latitude = data[:,7]
longitude = data[:,8]
yyyy,mm,dd=data[:,1],data[:,2],data[:,3]
hr,mn,ss = data[:,4],data[:,5],data[:,6]
index=data[:,0]

#call beachball function 
beachball(data)        

##########################################################################
# create a document element with multiple label style
kmlobj = KML.kml(
    KML.Document(
    )
)   
 
for j in range(len(yyyy)):  #create the ref icons we will use
    kmlobj.Document.append(     
        KML.Style(             
            KML.IconStyle(
                KML.Icon(
                    KML.href('%s.png'%str(int(index[j]))),
                    KML.scale(0.6),   #scale the beachball in googleEarth
                ),
                KML.heading(0.0),
            ),
        id='beach_ball_%i'%j    #gives the icon a ref that will be used later
        ),
    )

# add images to the document element
for i in range(len(yyyy)):
    datum = str(date.date(int(yyyy[i]),int(mm[i]),int(dd[i])))
    ev_time = str(date.time(int(hr[i]),int(mn[i]),int(ss[i])))
 
    kmlobj.Document.append(
        KML.Placemark(
            #~ KML.name('%s'%str(int(index[i]))),   #uncomment this to add a name to the placemark (will always appear in GoogleEarth)
            KML.ExtendedData(                       #I add information about the earthquake, it appears in a table ('info' : value)
                KML.Data(                           
                    KML.value('%s'%datum),          #add value of the specific info
                name ='date'                        #name of'info' you add.
                ),
                KML.Data(
                    KML.value('%s'%ev_time),        #add value of the specific info 
                name ='time'                        #name of 'info' you add.
                ),                                     #more data can be added, following the same structure (line 65-68)
            ),
            KML.styleUrl('#beach_ball_%i'%i),       #get the correct beachball in the directory as marker
            KML.Point(
                KML.coordinates(longitude[i],',',latitude[i]),
            ),
 
        ),
    )

print(etree.tostring(etree.ElementTree(kmlobj),pretty_print=True))

outfile = open('Focal_mechanism_devy.kml', 'wb') #save the kml structure code
outfile.write(etree.tostring(kmlobj, pretty_print=True))