#!/usr/bin/env python
# coding: utf-8

# # Plot interactive earthquake data using holoviews
# 
# This notebook downloads earthquake data from the USGS database and plots them in different ways using holoviews. Holoviews is wrapper around the bokeh framework that does a lot of the heavy lifting for you. Similar to Bokeh, it also allows for stand alone HTML code that can be embeded into a webpage, but you can also write apps that run on their server.
# 
# A demo for DataLink is provided, where highlighting data points in one plot also highlights them in another. Here, we use it to explore where earthquakes of a given magnitude and depth cluster on the map and vice versa.
# 
# Finally, a Gutenberg-Richter exercise has been added.
# 
# Major credit to http://examples.holoviews.org/Earthquake_Visualization.html, which this notebook builds from.
# 
# #### Warning: May need to clear kernel each time you want to rerun!

# In[1]:


from bokeh.io import output_file, show
import holoviews as hv
from holoviews import dim, opts
import datetime as dt
import numpy as np
import pandas as pd
from matplotlib.image import imread
from mpl_toolkits.basemap import Basemap
import pandas as pd

from io import BytesIO
import requests


# In[2]:


# hv.notebook_extension('bokeh', width=90)
hv.extension('bokeh')
renderer = hv.renderer('bokeh')


# ## Download data from USGS (limit of 20,000 earthquakes)

# In[3]:


query = dict(starttime="2015-05-13", endtime="2020-05-13", minmagnitude="5.0")
query_string = '&'.join('{0}={1}'.format(k, v) for k, v in query.items())
query_url = "http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&" + query_string
response = requests.get(query_url)
df = pd.read_csv(BytesIO(response.content),
                 parse_dates=['time'], index_col='time',
                 infer_datetime_format=True,
                 )
df['Date'] = [str(t)[:19] for t in df.index]
# Make longitude between 0 and 360 instead of -180 to 180
df['longitude'][df['longitude'] < 0] = df['longitude'][df['longitude'] < 0] + 360
df.head()


# ## Set up map structure
# 
# Export NASA's blue marble (https://visibleearth.nasa.gov/collection/1484/blue-marble) to holoviews RGB structure to use as plot background.

# In[4]:


kdims = ['Longitude', 'Latitude']
map = Basemap(llcrnrlon=0,llcrnrlat=-90.,urcrnrlon=359.9,urcrnrlat=90)
img = map.bluemarble()
imgss = np.flipud(img.get_array()[::5, ::5])
lon = np.linspace(img.get_extent()[0],img.get_extent()[1],imgss.shape[1])
lat = np.linspace(img.get_extent()[2],img.get_extent()[3],imgss.shape[0])
rgb = img.get_array()[::5, ::5]
blue_marble = hv.RGB((lon,lat,rgb[:,:,0],rgb[:,:,1],rgb[:,:,2]),
                     kdims=kdims)


# ## Pass the earthquake dataframe into the HoloViews Element
# kdims: independent variables  
# vdims: dependent variables

# In[5]:


# Pass the earthquake dataframe into the HoloViews Element. Call it "earthquakes"
earthquakes = hv.Points(df, kdims=['longitude', 'latitude'],
                        vdims=['place', 'Date', 'depth', 'mag'],
                        group='Earthquakes')


# ## Plot global seismicity as circles colored by depth with circle size corresponding to magnitude
# 

# In[6]:


# Style Plots
get_ipython().run_line_magic('opts', "Points.Earthquakes [responsive=True data_aspect=1 scaling_factor=1.5 colorbar=True toolbar='above' clabel='Depth (km)' tools=['hover', 'crosshair']] (color='depth' cmap='cool' size=1+1.2**(dim('mag')**1.3))")

# Plot blue marble first, then earthquakes
# The " * " opperator overlays components and the " + " creates two separate figures
plot1 = blue_marble * earthquakes.opts(title="Earthquakes > M5.0 (2015-2020)")

# Convert to bokeh figure then save using bokeh (Hugo likes the bokeh HTML better!)
bplot1 = renderer.get_plot(plot1).state
output_file('a1_global_seismicity.html')
show(bplot1)


# ## Explore magnitude-depth relationships with brushing and linking
# 
# Data in the two figures are linked. Selecting earthquakes in the map will highlight them in magnitude-depth space. It also works the other way around – select earthquakes of a desired magnitude and depth to see where they occur on the map.
# 

# In[7]:


# Need to import DataLink to connect the data in the two plots
from holoviews.plotting.links import DataLink
# Redefine earthquakes to avoid styling conflicts with the earlier one...
earthquakes2 = hv.Points(df, kdims=['longitude', 'latitude'],
                        vdims=['place', 'Date', 'depth', 'mag'],
                        group='Earthquakes2')
get_ipython().run_line_magic('opts', "Points.Earthquakes2 [scaling_factor=1.5] (size=1+1.2**(dim('mag')**1.3))")
get_ipython().run_line_magic('opts', "Points [height=int(250*1.3) width=int(400*1.3) tools=['hover', 'lasso_select', 'box_select', 'crosshair'] active_tools=['lasso_select']] (color='indianred', alpha=0.4, nonselection_color='indianred', selection_alpha=0.75, selection_color='deepskyblue')")

# Make magnitude-depth object
mag_depth = hv.Points(earthquakes2, kdims=['mag','depth'], group='Earthquake magnitude vs. depth').opts(size=10, xlabel='Magnitude', ylabel='Depth (km)')

# Link the data in the two objects
DataLink(earthquakes2,mag_depth)

# Make plot
plot2 = (blue_marble * earthquakes2.opts(title="Earthquakes > M5.0 (2015-2020)") + mag_depth.hist(dimension=['depth'])).cols(1)

# Convert to bokeh figure then save using bokeh (Hugo likes the bokeh HTML better!)
bplot2 = renderer.get_plot(plot2).state
output_file('a1_globalEQ_MagDepth.html')
show(bplot2)


# ## Gutenberg-Richter Relationship
# 
# The Gutenber-Richter relationship describes the frequency of occurence of an earthquake of a given magnitude. It is an empirical relationship given by:
# 
# $$ \log_{10}N = a - bM $$
# 
# $M$: Earthquake magnitude  
# $N$: Number of earthquakes greater than or equal to $M$  
# $b$: Shows distribution from small to large (slope of the line, usually ~1)  
# $a$: Measure of seismic activity (y-intercept)

# ## Calculate cumulative earthquake frequency: N

# In[8]:


# Define magnitude bins (includes left edge of first bin and right edge of last bin)
dmag = 0.1  # magnitude increments
mag_bins = np.arange(df.mag.min(),df.mag.max()+dmag,dmag)

# Count up number of earthquakes in each bin (in ascending order)
mag_counts = df['mag'].value_counts(bins=mag_bins, ascending=True)

# Cumulative sum of counts to get N
N_ascend = mag_counts.cumsum()

# Flip to descending magnitude to match mag_bins
N = N_ascend.values[::-1]

# Make sure all values of N are finite (since can't plot zero on log scale)
ind_finite = N > 0
N = N[ind_finite]
M = mag_bins[:-1][ind_finite]


# ## Estimate $a$ and Calculate G-R relationship estimate for $N_{GR}$
# 
# Solving the G-R relationship above for $a$ gives:
# 
# $$ a = \log_{10}N + bM$$
# 
# Then, we can solve for $N_{GR}$ as a function of $M$:
# 
# $$ N = 10^{(a - bM)} $$

# In[9]:


# Rough estimate for "a" assuming b=1
b = 1
a = np.median(np.log10(N) + b*M)

N_GR = 10**(a - b*M)


# ## Plot Guttenberg-Ricther relation using Bokeh

# In[10]:


import bokeh
from bokeh.models import ColumnDataSource, HoverTool
from bokeh.plotting import figure, output_notebook
from bokeh.models.formatters import FuncTickFormatter


# Define figure attributes
p = figure(plot_width=800, plot_height=500, sizing_mode = 'scale_width',
           x_range=([np.min(M)*0.98,np.max(M)*1.02]), y_range=([np.min(N)*0.5, np.max(N)*2]),
           y_axis_type="log",
           x_axis_label='M (magnitude)', y_axis_label='N (number of earthquakes greater than M)', 
           title='Gutenberg-Richter relationship for earthquakes > M5.0 (2015-2020)',
           tools=['save','box_zoom','xwheel_zoom','ywheel_zoom','reset','crosshair','pan'])

# Setup data source
source = ColumnDataSource(data={'N':N, 'N_GR':N_GR, 'M':M})

# Plot Gutenberg-Richter relationship
l1 = p.line('M', 'N_GR', source=source, line_color="firebrick",
             line_width=3, legend_label="N = 10^(a - bM)")
# p.add_tools(HoverTool(renderers=[l1], tooltips=[("Date", "@t_centstr"),("Noise Level","@disp_avg")]))

# Plot N
l2 = p.circle('M', 'N', source=source, line_color="black", fill_color="steelblue",
              alpha=0.3, line_width=1, size=15)  #legend_label="Hourly")
p.add_tools(HoverTool(renderers=[l2], tooltips=[("Magnitude", "@M"),("Earthquakes","@N")]))

# Modify log tick labels
p.yaxis[0].formatter = FuncTickFormatter(code="""
return tick
""")

# Save local file
output_file('a1_GutenbergRichter.html')
# Show plot in browser
show(p)


# In[ ]: