Skip to content

Commit

Permalink
Add Hawkes earthquake aftershocks study example
Browse files Browse the repository at this point in the history
We need to install basemap to make it work
  • Loading branch information
Mbompr committed Jun 13, 2018
1 parent cfe357d commit c31d333
Show file tree
Hide file tree
Showing 6 changed files with 6,382 additions and 2 deletions.
6,253 changes: 6,253 additions & 0 deletions doc/_static/example_data/earthquakes.csv

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions examples/earthquakes.csv
111 changes: 111 additions & 0 deletions examples/plot_hawkes_earthquake_example.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
"""
==============================================================
Study earthquake aftershocks propagation with Hawkes processes
==============================================================
This experiments aims to detect how aftershocks propagate near Japan. A Hawkes
process has been fit on data from
Ogata, Y., 1988.
Statistical models for earthquake occurrences and residual analysis
for point processes.
Journal of the American Statistical association, 83(401), pp.9-27.
On the left we can see where earthquakes have occurred and on the right
the propagation matrix, i.e. how likely a earthquake in a given zone will
trigger an aftershock in another zone.
We can observe than zone 1, 2 and 3 are tightly linked while zone 4 and
5 are more self-excited.
Note that this Hawkes process does not take the location of an earthquake
to recover this result.
Dataset `'earthquakes.csv'` is available
`here <../_static/example_data/earthquakes.csv>`_.
"""

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap

from tick.hawkes import HawkesEM
from tick.plot import plot_hawkes_kernel_norms
from tick.plot.plot_utilities import get_plot_color

df = pd.read_csv('earthquakes.csv')

lats = df.Latitude.values
lons = df.Longitude.values

# put timestamps in minutes and remove minimum
timestamps = pd.to_datetime(df.DateTime).values.astype(np.float64)
timestamps *= 1e-9 / 60
timestamps -= min(timestamps)

fig, ax_list = plt.subplots(1, 2, figsize=(8, 3.4))
ax_list[0].set_title('Earthquakes near Japan')

# Draw map
lon0, lat0 = (139, 34)
lon1, lat1 = (147, 44)
m = Basemap(projection='merc', llcrnrlon=lon0, llcrnrlat=lat0, urcrnrlon=lon1,
urcrnrlat=lat1, resolution='l', ax=ax_list[0])
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966', lake_color='#99ffff')

# Number of splits in the map
n_groups_lats = 5

group_lats = []
group_lons = []
group_timestamps = []
group_names = []

delta_lats = (lats.max() - lats.min()) / n_groups_lats
delta_lons = lons.max() - lons.min()

for group in range(n_groups_lats):
# Split events into groups
min_lat_group = lats.min() + group * delta_lats
max_lat_group = lats.min() + (group + 1) * delta_lats

mask = (min_lat_group <= lats) & (lats < max_lat_group)

group_lats += [lats[mask]]
group_lons += [lons[mask]]
group_timestamps += [timestamps[mask]]

# Draw earthquakes on map
x, y = m(group_lons[group], group_lats[group])
m.scatter(x, y, 0.3, marker='o', color=get_plot_color(group))

# Draw zone labels on map
group_name = 'Z{}'.format(group)
group_names += [group_name]

center_lat = 0.5 * (min_lat_group + max_lat_group)
x_zone, y_zone = m(lons.max() + 0.1 * delta_lons, center_lat)

zone_bbox_style = dict(boxstyle="round", ec=(0, 0, 0, 0.9), fc=(1., 1, 1,
0.6))
ax_list[0].text(x_zone, y_zone, group_name, fontsize=12, fontweight='bold',
ha='left', va='center', color='k', withdash=True,
bbox=zone_bbox_style)

# Fit Hawkes process on events
events = []
for g in range(n_groups_lats):
events += [group_timestamps[g]]
events[g].sort()

kernel_discretization = np.linspace(0, 10000 / 60, 6)
em = HawkesEM(kernel_discretization=kernel_discretization, tol=1e-9,
max_iter=1000, print_every=200)
em.fit(events)

plot_hawkes_kernel_norms(em, ax=ax_list[1], node_names=group_names)
ax_list[1].set_xticklabels(ax_list[1].get_xticklabels(), fontsize=11)
ax_list[1].set_yticklabels(ax_list[1].get_yticklabels(), fontsize=11)

fig.tight_layout()
plt.show()
13 changes: 11 additions & 2 deletions tick/plot/plot_hawkes.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@


def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
node_names=None):
node_names=None, ax=None):
"""Generic function to plot Hawkes kernel norms.
Parameters
Expand All @@ -33,6 +33,10 @@ def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
node names that will be displayed on axis.
If `None`, node index will be used.
ax : `matplotlib.axes`, default=None
If not None, the figure will be plot on this axis and show will be
set to False. Used only with matplotlib
Notes
-----
Kernels are displayed such that it shows norm of column influence's
Expand All @@ -50,7 +54,12 @@ def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
column_labels = ['$\\rightarrow {}$'.format(i) for i in node_names]

norms = kernel_object.get_kernel_norms()
fig, ax = plt.subplots()

if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
else:
fig = ax.figure
show = False

if pcolor_kwargs is None:
pcolor_kwargs = {}
Expand Down
4 changes: 4 additions & 0 deletions tools/travis/docker_run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,15 @@ set -e -x
cp -R /io src
cd src

apt-get update;
apt-get install -y libgeos*;

eval "$(pyenv init -)"

pyenv global ${PYVER}
pyenv local ${PYVER}

python -m pip install https://github.com/matplotlib/basemap/archive/v1.1.0.tar.gz
python -m pip install yapf --upgrade
python -m yapf --style tools/code_style/yapf.conf -i tick --recursive

Expand Down
2 changes: 2 additions & 0 deletions tools/travis/osx_install.sh
Original file line number Diff line number Diff line change
Expand Up @@ -31,5 +31,7 @@ python -m pip install -r requirements.txt
python -m pip install sphinx pillow
python -m pip install cpplint
python -m pip install tensorflow
brew install geos
python -m pip install https://github.com/matplotlib/basemap/archive/v1.1.0.tar.gz
pyenv rehash

0 comments on commit c31d333

Please sign in to comment.