Using folium - 2: Use customized icon (Plot moment tensor)

In this script, we will try to use a customize icon for marker on the folium map. It is actulaly really simple to use a customized icon, you can either refer to a file on your local drive, or give it a url so that it will link the file. The following scripts are showing both, the first one is ploting iron man who visit Berkeley (link the file online), and the second one is plotting a moment tensor from a file we generated on our drive. You can find the script at Qingkai's Github.

In [1]:

import folium

In [2]:

url_base = 'http://vignette3.wikia.nocookie.net/'
url_folder = 'the-adventures-of-the-gladiators-of-cybertron/'
url_file = 'images/6/67/Iron_Man.png/revision/latest?cb=20140720074813'
logo_url = url_base + url_folder + url_file
UCB_lat = 37.8719
UCB_lon = -122.2585

map_1 = folium.Map(location=[UCB_lat, UCB_lon], zoom_start=12,\
                      control_scale = True, tiles='Stamen Terrain')

icon = folium.features.CustomIcon(logo_url,\
                                  icon_size=(200, 200))

folium.Marker([UCB_lat, UCB_lon],
          popup='Iron man visit Berkeley',
          icon=icon).add_to(map_1)

map_1

Out[2]:

Plot moment tensor on the map¶

The following is how we can plot the moment tensors using obspy, and place it on the folium map using custom icon. Note that you need Obspy as the dependences. Let's plot the moment tensor for the 2014 south Napa earthquake here at Bay area, you can find details of this event at USGS. I won't show the maps here, but save it as a html file you can see it.

In [3]:

from obspy.imaging.beachball import beachball

# latitude and longitude of the earthquake
evla = 38.215
evlo = -122.312

# South Napa EQ moment tensor
mt = [247, 82, 8]
beachball(mt, size=200, linewidth=2, facecolor='b', outfile= './beachball.png')

# Add the USGS style map
url_base = 'http://server.arcgisonline.com/'
tiles = 'ArcGIS/rest/services/NatGeo_World_Map/MapServer/tile/{z}/{y}/{x}'
tileset = url_base + tiles
map_1 = folium.Map(location=[38, -122], zoom_start=8,\
                      control_scale = True, tiles=tileset, attr='USGS style')

# Add the Moment Tensor as an icon
icon_url = './beachball.png'
icon = folium.features.CustomIcon('/beachball.png',icon_size=(45, 45))

folium.Marker([evla, evlo],
          popup='Napa Earthquake',
          icon=icon).add_to(map_1)

map_1.save('moment_tensors.html')

Using folium - 1: Use different map tiles (USGS style map)

Folium is an awesome tool that can make interactive maps built on Leaflet. I found it very easy to use. I will create a series of blog entries to show some of the plots that I usually use for my research.
Let's start with a simple map that change the tile to the style that USGS using for their earthquake page. I think the map tile looks good, here's an example from USGS. Also, there's a good website that you can checkout different map tiles that can be used in folium - Here. For the following code, you can download the notebook from Qingkai's Github

In [1]:

import folium

In [2]:

# Add the USGS style tile
url_base = 'http://server.arcgisonline.com/ArcGIS/rest/services/'
service = 'NatGeo_World_Map/MapServer/tile/{z}/{y}/{x}'
tileset = url_base + service

map_1 = folium.Map(location=[37.8716, -122.2727], zoom_start=10,\
        control_scale = True, tiles=tileset, attr='USGS style')

map_1.add_children(
    folium.Marker([37.8716, -122.2727], popup = 'I am here'))
map_1

Out[2]:

Add progress bar to your python loops

Today I want to add a progress bar on my python loops, after I searched online, I found a nice simple module - tqdm can do the job easily. All you need do is install it via 'pip install tqdm', and then you can use it. You can find the code below here.

(Gif grab from tqdm github)

How to Write a Great Research Paper (great talk)

This is a great talk by Professor Simon Peyton Jones that I came across today! Even though it is aimed for Computer Science Conference papers, it is applicable to other fields as well! I will change my way to write paper!

Learning python basics

I have one friend ask me today for resources to learn python. The followings are the best ones I can think of, hope these will also useful to you.

Book:

My faverate book, and will let my students in the future to have as a reference is (if I become a faculty in the future ^)^):

Effective computation in physics

This book not only covers a lot of things in python, but also a lot popular tools that used in science. 

Another good book is: Automate boring stuff with python

This book is very practical, and you can automate a lot of boring stuff with python! Go python!

Webtutorial:

The notebook gallery has lots of interesting notebooks that you can play with, this is the best way to learn python in my opinion, learn from doing. You can start with the python tutorials. A gallery of interesting Ipython Notebooks

Another great website that has a collection of good tutorials is: learning python

Videos:

Scipy conference has a lot of great tutorials online, you can checkout the following link, most of the tutoirals have a link to their repo on github, so you can do this following their video. Scipy 2015 videos

If you go through all these, I think you will be quite fluent in python!

Plotting historical moment tensors from GCMT directly

I created this script for our earthquake of this week (There are some other scripts for different purposes, I may add them sometime later, but you can find them at my Github). The purpose of this script is that, a lot of times, when we look at an earthquake occured recently, we want to see the historical moment tensors in this region. These historical moment tensors can tell us a lot of things, i.e. tectonic setting, fault information, etc. Then using this script, all you need do is to specify the starttime, endtime, and region (for region, right now, I just use one recent earthquake location as the center, and expand them both in latitude and longitude a certain degree).
For example, the 2016 M7.8 Ecuador, if we want to see the historical moment tensors in this region, we can do as the followings. You can find the script on Qingkai's Github.

In [1]:

import warnings
warnings.filterwarnings('ignore')
from obspy import UTCDateTime
import urllib
import urllib2
import numpy as np
from collections import OrderedDict
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from obspy.imaging.beachball import beach
from matplotlib.lines import Line2D
from BeautifulSoup import BeautifulSoup
%pylab inline

Populating the interactive namespace from numpy and matplotlib

In [3]:

def get_hist_mt(t1, t2, llat = '-90', ulat = '90', llon = '-180', ulon = '180', 
    lmw = 0, umw = 10, evla = None, evlo = None, step = 2.0, list_ = '6'):
    '''
    Function to query the GCMT and save the data for ploting. If evla and 
    evlo are provided, it will use this location as the center and adding 
    step to both latitude and longitude for the box. 
    
    t1 - starttime, example: UTCDateTime("1979-01-01")
    t2 - endtime, example: UTCDateTime("2016-01-01")
    llat - lower latitude
    ulat - upper latitude
    llon - left longitude
    ulon - right longitude
    lmw - the minimum magnitude to search
    umw - the maximum magnitude to search
    evla - latitude of the current earthquake
    evlo - longitude of the current earthquake
    step - from current earthquake location, expand lat/lon for a box
    list_ - format of data you want to return from the GCMT
    '''
    
    yr = t1.year
    mo = t1.month
    day = t1.day
    oyr = t2.year
    omo = t2.month
    oday = t2.day
    mat = {}
    locs = locals()  
    
    base_url = 'http://www.globalcmt.org/cgi-bin/globalcmt-cgi-bin/CMT4/form'
    
    #note here, we must use the ordered Dictionary, so that the order in the 
    #url is exactly the same order
    param = OrderedDict()
    param['itype'] = 'ymd'
    param['yr'] = yr
    param['mo'] = mo
    param['day'] = day
    param['otype'] = 'ymd'
    param['oyr'] = oyr
    param['omo'] = omo
    param['oday'] = oday
    param['jyr'] = '1976'
    param['jday'] = '1'
    param['ojyr'] = '1976'
    param['ojday'] = '1'
    
    param['nday'] = '1'
    param['lmw'] = str(lmw)
    param['umw'] = str(umw)
    param['lms'] = '0'
    param['ums'] = '10'
    param['lmb'] = '0'
    param['umb'] = '10'
    
    # now specify the region box
    if evla and evlo is not None:
        llat = evla - step
        ulat = evla + step
        llon = evlo - step
        ulon = evlo + step
    
    # save parameter for query
    param['llat'] = llat
    param['ulat'] = ulat
    param['llon'] = llon
    param['ulon'] = ulon
    
    param['lhd'] = '0'
    param['uhd'] = '1000'
    param['lts'] = '-9999'
    param['uts'] = '9999'
    param['lpe1'] = '0'
    param['upe1'] = '90'
    param['lpe2'] = '0'
    param['upe2'] = '90'
    param['list'] = list_
    
    # build the URL
    url = "?".join((base_url, urllib.urlencode(param)))
    print url
    
    # grab data and parse it
    page = urllib2.urlopen(url)
    parsed_html = BeautifulSoup(page)
    mecs_str = parsed_html.findAll('pre')[1].text.split('\n')

    # string to array
    def convertString(mecs_str):
        return map(float, mecs_str.split()[:9])
        
    psmeca = np.array(map(convertString, mecs_str))
    
    # save the results for plotting
    mat['psmeca'] = psmeca
    mat['url'] = url
    mat['range'] = (llat, ulat, llon, ulon)
    mat['evloc'] = (evla, evlo)
    return mat
    
def plot_hist_mt(psmeca_dict, figsize = (16,24), mt_size = 10, \
                 pretty = False, resolution='l'):
    '''
    Plot the historical moment tensor from the query of GCMT
    
    psmeca_dict - dictionary that returned from get_hist_mt function
    figsize - tuple of the size of the figure
    mt_size - size of the moment tensor on the plot
    pretty - boolean, whether want to plot nice maps
    resolution - low or high as you want
    '''
    
    if psmeca_dict['psmeca'] != []:
        psmeca = psmeca_dict['psmeca']
        #get the latitudes, longitudes, and the 6 independent component
        lats = psmeca[:,1]
        lons = psmeca[:,0]
        focmecs = psmeca[:,3:9]
        depths =  psmeca[:,2]    
        (llat, ulat, llon, ulon) = psmeca_dict['range'] 
        evla = psmeca_dict['evloc'][0]
        evlo = psmeca_dict['evloc'][1]

        plt.figure(figsize=figsize)
        m = Basemap(projection='cyl', lon_0=142.36929, lat_0=38.3215, 
            llcrnrlon=llon,llcrnrlat=llat,urcrnrlon=ulon,urcrnrlat=ulat,\
                    resolution=resolution)
    
        m.drawcoastlines()
        m.drawmapboundary()
    
        if pretty:    
            m.etopo()
        else:
            m.fillcontinents()
    
        llat = float(llat)
        ulat = float(ulat)
        llon = float(llon)
        ulon = float(ulon)
    
        m.drawparallels(np.arange(llat, ulat, (ulat - llat) / 4.0), \
                        labels=[1,0,0,0])
        m.drawmeridians(np.arange(llon, ulon, (ulon - llon) / 4.0), \
                        labels=[0,0,0,1])   
    
        ax = plt.gca()
        x, y = m(lons, lats)
    
        for i in range(len(focmecs)):
        
            if depths[i] <= 50:
                color = '#FFA500'
                #label_
            elif depths[i] > 50 and depths [i] <= 100:
                color = 'g'
            elif depths[i] > 100 and depths [i] <= 200:
                color = 'b'
            else:
                color = 'r'
        
            index = np.where(focmecs[i] == 0)[0]
        
            #note here, if the mrr is zero, then you will have an error
            #so, change this to a very small number 
            if focmecs[i][0] == 0:
                focmecs[i][0] = 0.001
        
            try:
                b = beach(focmecs[i], xy=(x[i], y[i]),width=mt_size, \
                          linewidth=1, facecolor=color)
            except:
                pass
            
            b.set_zorder(10)
            ax.add_collection(b)
        
        # add the current earthquake
        x_0, y_0 = m(evlo, evla)
        m.plot(x_0, y_0, 'r*', markersize=25, zorder = 10) 
        
        # add the legend
        circ1 = Line2D([0], [0], linestyle="none", \
                marker="o", alpha=0.6, markersize=10, markerfacecolor="#FFA500")
        circ2 = Line2D([0], [0], linestyle="none", \
                marker="o", alpha=0.6, markersize=10, markerfacecolor="g")
        circ3 = Line2D([0], [0], linestyle="none", \
                marker="o", alpha=0.6, markersize=10, markerfacecolor="b")
        circ4 = Line2D([0], [0], linestyle="none", \
                marker="o", alpha=0.6, markersize=10, markerfacecolor="r")
        plt.legend((circ1, circ2, circ3, circ4), \
                   ("depth $\leq$ 50 km", "50 km $<$ depth $\leq$ 100 km", \
                    "100 km $<$ depth $\leq$ 200 km", "200 km $<$ depth"), \
                   numpoints=1, loc=3)
        plt.show()
    else:
        print 'No historical MT found!'

In [6]:

t1 = UTCDateTime("1979-01-01")
t2 = UTCDateTime("2016-01-01")

psmeca = get_hist_mt(t1, t2, evla = 0.371, evlo= -79.940, step = 8, lmw = 0)
plot_hist_mt(psmeca, figsize = (12, 10),  mt_size = 0.4, pretty = True, resolution='l')