Easily create .kml file to plot focal mechanisms in Google Earth

Plotting focal mechanisms in Google Earth is an annoying task to do manually. I wrote a script to convert a text file to beach balls using the obspy package and then to create the corresponding KML file. Each ball is a symbol, used by one Placemark. Below, you’ll see a view of the beachball plot result. The behaviour is the same as for normal placemarks. If beachballs are overlapping, you can spread them by simply clicking on the cluster. When the icon scale factor is small, by zooming in, the number of overlapping icons will decrease.  It is easy to add information in the table that appears when clicking on the beachball.

The full code is after the break:

#!/usr/bin/python

#~ warning !: the script will draw beachball and creates the kml file to plot them in googleearth. All created element 
#~ will be placed in the same directory as the one containing this script. The kml file and the beachballs have to 
#~ stay in the same directory. 

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

def beachball(data):
    """function to draw beachball using obspy library"""
    import matplotlib.pyplot as plt
    from obspy.imaging.beachball import Beachball

    strike,dip,rake=data[:,9],data[:,10],data[:,11]
    event=data[:,0]                 #index to identify the beachball created

    for j in range(len(strike)):
        beach = Beachball([strike[j],dip[j],rake[j]],outfile=str(int(event[j])),\
            facecolor='black',edgecolor='black') #create and save beachball in a outfile in the directory where the .py file is 

data=np.loadtxt('path\table_w_your_data.txt',\
    skiprows=1) #import your data. (table format: index yyyy mm dd hr mn ss lat lon strike dip rake) the first row isn't used
latitude = data[:,7]
longitude = data[:,8]
yyyy,mm,dd=data[:,3],data[:,2],data[:,1]
hr,mn,ss = data[:,4],data[:,5],data[:,6]
index=data[:,0]

beachball(data)         #call beachball function   --> put in comment this line if you don't want to draw again all beachballs

###############################################################################################################################
# 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= file('Focal_mechanism_devy.kml','w') #save the kml structure code
outfile.write(etree.tostring(kmlobj, pretty_print=True))

 

17 thoughts on “Easily create .kml file to plot focal mechanisms in Google Earth”

  1. Dimitri,
    nice tool!
    One comment: you should change line 28 into:
    yyyy,mm,dd=data[:,1],data[:,2],data[:,3]

    Is it possible to add a line to keep the beachball fixed to the north? Now if you rotate in Google Earth, also the beachball turns.

    Cheers, Koen

  2. Dimitri,

    Very nice tool! I’m interested to use it, but I think my data is not formatting in a way that numpy likes, I’m getting the following error:

    Traceback (most recent call last):
    File “C:\Python27\DrawFMinGE.py”, line 26, in
    latitude=data[:,7]
    IndexError: invalid index

    Can you recommend exactly how the .txt file should be formatted?

    Thanks,

    Erik

  3. Hello Erik,

    The file format is explained in comment line 25. The 8th column (data[:,7] in python counting system) is the latitude.

    table format: index yyyy mm dd hr mn ss lat lon strike dip rake
    corresponding column value in python 0 1 2 3 4 5 6 7 8 9 10 11

    let me know if it you have any other problems,

    Dimitri

  4. Dimitri,

    Thanks for responding. I don’t have much experience with Python, so I spent significant time yesterday educating myself on some basics. I still have no idea why I’m getting that “invalid index” error. I’m sure I’m missing something very simple and basic, but I don’t know what it is.

    The .txt file I’m trying to use only contains information on one focal mechanism (a test case). The entire contents of my .txt file are:

    index yyyy mm dd hr mn ss lat lon strike dip rake
    1 1978 11 29 14 37 50 30.57 -114.73 189 35 122

    In case formatting didn’t come through, there is a single space between each column a single return separating row 1 (header) from row 2 (data). What have I got wrong here?

    Many thanks,

    Erik

  5. Did anyone ever figure a way to keep the focal mechanisms fixed north?

    Also, I get the following error when I uncomment line 68 (to add a placemark):
    ~ KML.name(‘%s’%str(int(index[i]))), #uncomment this to add a name to the placemark (will always appear in GoogleEarth)

    TypeError: bad operand type for unary ~: ‘lxml.objectify.ObjectifiedElement’

    Any thoughts?

    Cheers,

    JH

  6. Hello Mustafa !

    the format is:

    Line 25: #table format: index yyyy mm dd hr mn ss lat lon strike dip rake) the first row isn’t used

    let me know how it goes

  7. Actually I have made it. then the next message appeared as
    python beachball.py
    Traceback (most recent call last):
    File “beachball.py”, line 26, in
    latitude = data[:,7]
    IndexError: too many indices

  8. Mustafa,

    it should work if you respected the format I mentionned. it might be a delimiter issue, or you don’t have 8 columns in the text file.

    You can test it works by adding a “print data” just above the line 26.

  9. sorry for that.
    I am trying to respect the format.. print data give me this output.
    python beachball.py
    [ 1.00000000e+00 1.97800000e+03 1.10000000e+01 2.90000000e+01
    1.40000000e+01 3.70000000e+01 5.00000000e+01 3.00000000e+01
    -1.14000000e+02 1.89000000e+02 3.50000000e+01 1.22000000e+02]
    Traceback (most recent call last):
    File “beachball.py”, line 27, in
    latitude = data[:,7]
    IndexError: too many indices

  10. After I delete the first character for index. the out put is

    ustafa@mcomoglu:~$ python beachball.py
    [ 2014. 11. 29. 14. 37. 50. 36. 33. 189. 35.
    122.]
    Traceback (most recent call last):
    File “beachball.py”, line 27, in
    latitude = data[:,7]
    IndexError: too many indices

  11. Thanks for the clarifying the things. It works fine now. One thing already we discussed is if it s possible It would be fine to give a name with the new version to the png files a combination of YYYYMMDDHHMMSS? I am adding the name as I described semi automatically.
    Thanks again.

  12. Thank you for this really useful program. I learned enough about Python, which I had not used before, in a few hours to run the program and adapt it to add a few features (e.g. coloured beachballs for different fault types, display depth in info box). It is also easy to convert the kml output file to shapefile format for ArcGIS, using the Arc Toolbox, join the new shapefile to the original data file, and display the beachballs in ArcMap with their attributes.

Leave a Reply

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

*