diff --git a/Group4_project.ipynb b/Group4_project.ipynb new file mode 100644 index 0000000..d939a26 --- /dev/null +++ b/Group4_project.ipynb @@ -0,0 +1,1106 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Earthquake Data Analysis\n", + "by Rosset Lorenzo, Simionato Giuseppe, Sinigaglia Nicholas and Zanola Andrea" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The aim of this project is to analize a dataset containing records of earthquakes in order to to verify some statistical distributions often used to describe their behaviour, as for example the Gutenberg–Richter Law and the Omori Law of aftershocks. \\\n", + "The dataset used for this analysis is composed of 110271 events detecteced in Southern California in a period of time that goes from Jan 1st 1982 to Jun 29th 2011." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import math\n", + "import pandas as pd\n", + "import pandas_bokeh\n", + "import matplotlib.pyplot as plt\n", + "import datetime\n", + "import imageio\n", + "import networkx as nx\n", + "import os\n", + "import seaborn as sns\n", + "import IPython as ip\n", + "from matplotlib.ticker import MultipleLocator\n", + "from scipy import optimize,stats\n", + "from bokeh.plotting import figure, show\n", + "from bokeh.io import output_notebook\n", + "from bokeh.tile_providers import CARTODBPOSITRON, get_provider\n", + "from bokeh.transform import factor_cmap, factor_mark, linear_cmap, log_cmap\n", + "from bokeh.palettes import Spectral6\n", + "from bokeh.models import markers, Segment, Plot\n", + "from bokeh.models.expressions import CumSum\n", + "\n", + "output_notebook()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "More in detail, for every earthquake in the catalog 7 features are given in 7 columns:\n", + "- 1st column contains the index of the event, from 0 to 112070;\n", + "- 2nd column contains the index of the previous event that triggered the considered one, index is -1 if no ancestor is found. These index are associated to earthquakes by an algorithm based on pandemic models;\n", + "- 3rd column contains the time of the event, in seconds, from 0:00 of Jan 1st, 1982;\n", + "- 4th column contains the magnitude of the event;\n", + "- 5th, 6th and 7th columns contain the x,y,z euclidean coordinates of the event. These coordinates are used in analysis to obtain latitude, longitude and depth of every event." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = pd.read_csv('SouthCalifornia-1982-2011_Physics-of-Data.dat', sep=' ', names=['n','pointer', 't', 'mag', 'x', 'y', 'z'])\n", + "df['date'] = pd.to_datetime(df['t'], unit='s', origin='01/01/1982')\n", + "df=df.sort_values(by='t')\n", + "\n", + "df[1000 : 1005]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following functions allow to convert the cartesian coordinates in latitude, longitude and depth, and also in the Mercator coordiante system using the Elliptical Mercator Projection." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def merc_x(lon):\n", + " r_major=6378137.000\n", + " return r_major*math.radians(lon)\n", + "\n", + "def merc_y(lat):\n", + " if lat>89.5:lat=89.5\n", + " if lat<-89.5:lat=-89.5\n", + " r_major=6378137.000\n", + " r_minor=6356752.3142\n", + " temp=r_minor/r_major\n", + " eccent=math.sqrt(1-temp**2)\n", + " phi=math.radians(lat)\n", + " sinphi=math.sin(phi)\n", + " con=eccent*sinphi\n", + " com=eccent/2\n", + " con=((1.0-con)/(1.0+con))**com\n", + " ts=math.tan((math.pi/2-phi)/2)/con\n", + " y=0-r_major*math.log(ts)\n", + " return y" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 6371000. #Earth radius\n", + "df['lat'] = np.arcsin(df['z']/R)*180./np.pi\n", + "df['long'] = np.arctan2(df['y'], df['x'])*180./np.pi\n", + "df['dep'] = R-np.sqrt(df['x']**2+df['y']**2+df['z']**2)\n", + "df['x_merc'] = df['long'].apply(merc_x)\n", + "df['y_merc'] = df['lat'].apply(merc_y)\n", + "#df.describe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Earthquakes visualization on the map" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tile_provider = get_provider(CARTODBPOSITRON)\n", + "x_range, y_range = (-13600000, -12500000), (3480000, 4500000)\n", + "p1 = figure(x_range=x_range, y_range=y_range,\n", + " x_axis_type='mercator', y_axis_type='mercator', title='Map of earthquake distribution')\n", + "p1.add_tile(tile_provider)\n", + "p1.title.text_font_size = \"30px\"\n", + "\n", + "#Creating a copy of the data for plotting\n", + "df1 = df.sort_values(by='mag') #sorting w.r.t. magnitude for better visualization\n", + "df1['mag'] = df1['mag'].round(0) #round mag values\n", + "df1['mark_size'] = df1['mag']**1.4\n", + "df1['alpha'] = df1['mag']/10.\n", + "\n", + "mapper = log_cmap(field_name='mag', palette=Spectral6, low=2, high=7)\n", + "p1.scatter('x_merc', 'y_merc', source=df1, legend_field='mag', fill_alpha='alpha',\n", + " size='mark_size', fill_color=mapper, line_color=None)\n", + "show(p1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to visualize events over time, we made some 2D gifs. They represent only a small interval of time compared to the almost 30 years covered by dataset, but are useful to visualize concatenation of earthquakes and how unexpected a main event can be (as the 28 June 1992 earthquake).\n", + "\n", + "Every earthquake is shown as a red dot with size proportional to magnitude and that fades while time moves forward. \\\n", + "The x and y coordinates used are respectevely longitude and latitude." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# GIF CREATION\n", + "def gif_maker(df,data,n_frames,dt,index_min,index_max,gif_name):\n", + "# df = dataframe\n", + "# dt = time between two consecutive frames\n", + "# data = gif initial data \n", + "# n_frames = number of frames\n", + "# index_min , index_max = parameters used for read only the part of dataset that cointains the event of interests. \n", + "# They are choosed afterwards and improve performances\n", + "\n", + " # offset used to convert seconds in mm/dd/yyyy \n", + " T0_dataset = datetime.datetime(1982, 1, 2, 2, 0)\n", + " delta_T = (T0_dataset - datetime.datetime(1970,1,1)).total_seconds()\n", + " # gif initial time\n", + " t0 = (data - T0_dataset).total_seconds()\n", + " \n", + " # data needed\n", + " T = np.array(df[\"t\"].tolist())\n", + " M = np.array(df[\"mag\"].tolist())\n", + " Lat = np.array(df[\"lat\"].tolist())\n", + " Long = np.array(df[\"long\"].tolist())\n", + "\n", + " frames = []\n", + " for step in range(n_frames):\n", + " x_points = []\n", + " y_points = []\n", + " size_points = []\n", + " alpha_points = []\n", + " n_events = 0\n", + "\n", + " # looking over all dataset is not necessary, for this reason index_min and index_max are used\n", + " for i in range(index_min, index_max,1):\n", + " \n", + " # also events \"in the past\" are registred in order to make them fade using alpha parameter \n", + " if(T[i]>(step-6)*dt+t0 and T[i]<(step+1)*dt+t0):\n", + " x_points.append(Long[i])\n", + " y_points.append(Lat[i])\n", + " # stronger earthquakes have bigger dots\n", + " size_points.append((M[i]-0.75)**3.33)\n", + " n_events = n_events + 1 \n", + " # alpha parameter is associated, new events have aplha = 1 and are the most visible\n", + " if(T[i]>(step-6)*dt+t0 and T[i]<(step-5)*dt+t0): alpha_points.append(0.1)\n", + " elif(T[i]>(step-5)*dt+t0 and T[i]<(step-4)*dt+t0): alpha_points.append(0.2)\n", + " elif(T[i]>(step-4)*dt+t0 and T[i]<(step-3)*dt+t0): alpha_points.append(0.3)\n", + " elif(T[i]>(step-3)*dt+t0 and T[i]<(step-2)*dt+t0): alpha_points.append(0.4)\n", + " elif(T[i]>(step-2)*dt+t0 and T[i]<(step-1)*dt+t0): alpha_points.append(0.6)\n", + " elif(T[i]>(step-1)*dt+t0 and T[i]<(step)*dt+t0): alpha_points.append(0.8)\n", + " elif(T[i]>(step)*dt+t0 and T[i]<(step+1)*dt+t0): alpha_points.append(1)\n", + " \n", + " # implementin colors (red)\n", + " rgba_colors = np.zeros((n_events,4))\n", + " rgba_colors[:,0] = 1\n", + " rgba_colors[:, 3] = alpha_points \n", + "\n", + " # adding background, which has same x_range and y_range in term of Lat and Long\n", + " img = plt.imread(\"gif_background.png\")\n", + " fig, ax = plt.subplots(figsize=(8, 8))\n", + " ax.imshow(img, extent=[-122.8, -113.46, 29.8, 38])\n", + "\n", + " # creating frame, which is a graph\n", + " ax.scatter(x=x_points,y=y_points,s=size_points,color=rgba_colors)\n", + " name = str(step) + \".png\"\n", + " title = datetime.datetime.fromtimestamp(delta_T+t0+(step+1)*dt).strftime(\"%B %d, %Y %H:%M:%S\")\n", + " ax.set_title(title)\n", + " ax.set_xlim(-122.8, -113.46)\n", + " ax.set_ylim(29.8, 38)\n", + " fig.savefig(name)\n", + " plt.clf()\n", + " plt.close()\n", + " \n", + " # adding file as gif name\n", + " frames.append(name)\n", + " frames.append(name)\n", + " \n", + " # creating gif \n", + " with imageio.get_writer(gif_name, mode='I') as writer:\n", + " for filename in frames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image) \n", + " # deleting file \n", + " for file in set(frames):\n", + " os.remove(file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# with these parameters we can see 28 June 1992 earthquake in detail \n", + "gif_maker(df,datetime.datetime(1992, 6, 27, 22, 0),60,1800,38000,43000,'June_1992_detailed_earthquake.gif')\n", + "ip.display.Image('June_1992_detailed_earthquake.gif')\n", + "\n", + "# with these parameters we can see 22 April 1992 foreshocks \n", + "#gif_maker(df,datetime.datetime(1992, 4, 13, 20, 0),60,40000,35000,40000,'April_1992_foreshocks.gif')\n", + "#ip.display.Image('April_1992_foreshocks.gif')\n", + "\n", + "# with these parameters we can see 28 June 1992 earthquake on a larger time scale\n", + "#gif_maker(df,datetime.datetime(1992, 6, 20, 0, 0),60,42000,35000,51000,'June_1992_earthquake.gif')\n", + "#ip.display.Image('June_1992_earthquake.gif')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Frequency plot\n", + "An interesting feature we can visualize is the distribution of the number of events over time. In the plot below is shown also the magnitude of the single earthquakes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mag_mask = df['mag'] > 3\n", + "c_index = df['mag'][mag_mask].round(0).astype('int32').values - 2\n", + "spectral6 = np.asarray(Spectral6)\n", + "\n", + "gridsize = (3, 1)\n", + "fig = plt.figure(figsize=(12,8))\n", + "\n", + "ax1 = plt.subplot2grid(gridsize, (0, 0), colspan=1, rowspan=2)\n", + "ax1_color = 'xkcd:coral'\n", + "ax1.plot(df['date'], df['n'], lw=2, c= ax1_color)\n", + "#ax1.grid(ls='--')\n", + "ax1.set_ylabel('$Cumulative~sum$', size=15, color=ax1_color)\n", + "ax1.tick_params(axis='y', labelcolor=ax1_color)\n", + "\n", + "ax2 = ax1.twinx()\n", + "ax2_color = 'xkcd:peach'\n", + "ax2.set_ylabel('$Number~of~events$', color=ax2_color, size=15)\n", + "ax2.hist(df['date'], histtype='step', bins=29*365, color=ax2_color);\n", + "ax2.tick_params(axis='y', labelcolor=ax2_color)\n", + "\n", + "ax3 = plt.subplot2grid(gridsize, (2, 0))\n", + "ax3.scatter(df['date'][mag_mask], df['mag'][mag_mask], c=spectral6[c_index])\n", + "ax3.vlines(df['date'][mag_mask], 0, df['mag'][mag_mask], colors=spectral6[c_index])\n", + "ax3.set_yscale('log')\n", + "ax3.set_yticks([3,4,5,6,7])\n", + "ax3.set_yticklabels([3,4,5,6,7])\n", + "ax3.set_xlabel('$Year$', size=15)\n", + "ax3.set_ylabel('$m$', size=15);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It can be easily seen how soon after a very powerful earthquake the density of events over time suddenly increases. This is a clear sign of the fact that earthquackes are correlated over space and time, and tend to organize in swarms." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Networks\n", + "In what follows we are going to visualize the hierarchic networks that arise from two of the main events of the dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_swarm(df, date_range, magth, central_event):\n", + " mask1 = (df['date'] > date_range[0]) & (df['date'] < date_range[1]) & (df['mag'] > magth)\n", + " dn = df[mask1]\n", + " net = nx.from_pandas_edgelist(dn, source='n', target='pointer', create_using=nx.MultiGraph())\n", + " main_nodes = nx.node_connected_component(net, central_event) #nodes of the connected component of the central event. Needs to be converted\n", + " main_net_temp = nx.subgraph(net, main_nodes) #network of the principal connected component\n", + " df_mask = np.asarray(nx.nodes(main_net_temp)) #list of nodes of the connected components\n", + " df_mask = np.delete(df_mask, df_mask == -1)\n", + " main_df = df.loc[df_mask] # df with just the interesting nodes\n", + " return main_df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def network(df, date_range, magth, central_event, directed=False):\n", + " main_df = get_swarm(df, date_range, magth, central_event)\n", + " kind = 0\n", + " if (directed):\n", + " kind = nx.MultiDiGraph()\n", + " else: kind = nx.MultiGraph()\n", + " main_net = nx.from_pandas_edgelist(main_df, source='pointer', target='n', create_using=kind)\n", + " colors = main_df['mag'].round(0)\n", + " colors = colors.reindex(main_net.nodes()) # this is needed for assigning each node to the corresponding color\n", + " return (main_net, colors)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "net1, colors1 = network(df, ('01/05/1992', '01/05/1993'), 3.7, 39805)\n", + "net2, colors2 = network(df, ('01/09/1999', '01/09/2000'), 3.5, 75422)\n", + "fig_net = plt.figure(figsize=(16,8))\n", + "fig_net.add_subplot(121)\n", + "nx.draw(net1, node_size=30, node_color=colors1, cmap=plt.cm.viridis)\n", + "fig_net.add_subplot(122)\n", + "nx.draw(net2, node_size=50, node_color=colors2, cmap=plt.cm.inferno)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Originshocks, Foreshocks, Aftershocks and Gutenberg-Richter Law" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The main goal of this section is to study the behaviour of different types of earthquakes depending on their magnitude. \\\n", + "It's known in literature that exist different phases in the dynamic of earthquakes, that create a kind of earthquakes classification:\n", + "1. At the beginning, due to the stress between different plates an origin event is produced and we call it Originshock. They are marked with -1 in the dataset.\n", + "2. After that an earthquake swarm starts, and maybe it will produce a strong event, a Mainshock;\n", + "3. The earthquakes happened betwenn a Mainshock and its origin are called Foreshocks.\n", + "4. Then, after the main, another earthquake swarm starts and these shocks are called Aftershocks.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First of all, we tried to visualize the distribution followed by the temporal and spatial distances between a Mainshocks and its Origin.\n", + "We considered a Mainshock every earthquake with a magnitude above an arbitrary $m_{th}$ threshold, setted equals to 5, and with the Back Tree Algorithm we found its origin." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Back Tree Algorithm\n", + "def bta(mth):\n", + " mask=(df['mag']>=mth) & (df['pointer']!=-1)\n", + " i1 = df['pointer'][mask].to_numpy()\n", + " origin=[]\n", + " originv=[]\n", + " for i in range(len(i1)):\n", + " mask1=[]\n", + " mask1=[i1[i]]\n", + " while mask1[-1]!=-1:\n", + " i2=df['pointer'][mask1[-1]] \n", + " mask1.append(i2)\n", + " origin=np.insert(origin,0,mask1[-2])\n", + " originv=np.insert(originv,0,mask1[:-1])\n", + " origin=np.flip(origin)\n", + " #all indices that go backwards to generator event \n", + " originv=np.unique(np.flip(originv))\n", + " return origin, originv\n", + "\n", + "#other useful function for fitting\n", + "def q(x,a,b,c):\n", + " q=a+b*x+c*x**2\n", + " return q\n", + "def f(x,a,b):\n", + " f=a+b*x\n", + " return f\n", + "#function that helps to plot the last figures\n", + "def s(i,j):\n", + " if i==1 and j==0:\n", + " p=2\n", + " elif i==1 and j==1:\n", + " p=3\n", + " else:\n", + " p=i+j\n", + " return p" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "#--------------------------------------ORIGINSHOCK---------------\n", + "fig1, ax1 = plt.subplots(figsize=(12,12))\n", + "binss=np.linspace(2,7,21)\n", + "counts1, bins1,_ = ax1.hist(df['mag'][df['pointer']==-1], bins=binss,visible=0,density=1,cumulative=-1)\n", + "mask=counts1>0\n", + "logmB=np.log10(counts1[mask])\n", + "xB=((bins1[:-1]+bins1[1:])/2)[mask]\n", + "B=optimize.curve_fit(f,xB,logmB,p0=[6,-1]) \n", + "resB = stats.linregress(xB, logmB)\n", + "\n", + "#---------------------------------ORIGIN-MAIN PDF---------------\n", + "mth=5 #---CHOOSE MTH HERE-----\n", + "origin,_ =bta(mth)#we call the back tree-alg\n", + "mask=(df['mag']>=mth) & (df['pointer']!=-1) \n", + "dtmc=df['t'][mask].to_numpy()-df['t'][origin].to_numpy()\n", + "x1=(df['x'][mask].to_numpy()-df['x'][origin].to_numpy())**2\n", + "y1=(df['y'][mask].to_numpy()-df['y'][origin].to_numpy())**2\n", + "z1=(df['z'][mask].to_numpy()-df['z'][origin].to_numpy())**2\n", + "drmc=(x1+y1+z1)**0.5\n", + "dfr=pd.DataFrame({'dtmc': dtmc/(60*60*24), 'drmc': drmc/1000,'mag':df['mag'][mask]})\n", + "lev=[0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1] #levels of probability\n", + "\n", + "sns.kdeplot(x=np.log10(dfr['dtmc']), y=np.log10(dfr['drmc']),fill='true',ax=ax1,levels=lev,color='blue')\n", + "ax1.set_xlabel('Time Distance Origin-Main',size=12)\n", + "ax1.set_ylabel('Spatial Distance Origin-Main',size=12)\n", + "ax1.set_title('2D-PDF for Distance and Time $m\\geq%.0f$ for Origin-Main Earthquakes'%mth,size=20)\n", + "ax1.set_xticks([-3.15,-1.38,0,0.85,1.48,2.56,3.56]) #-3.15 -> 1min -- -4.93 -> 1sec\n", + "ax1.set_xticklabels(['1min','1h','1d','1w','1mth','1y','10y']);\n", + "ax1.set_yticks([-1,0,1,2,3,4])\n", + "ax1.set_yticklabels(['100m','1km','10km','100km','1000km','10000km'])\n", + "ax1.set_xlim(-5,4.5)\n", + "ax1.set_ylim(-1.8,3)\n", + "ax1.grid(which='both',b=1, linestyle='dashed')\n", + "\n", + "mask1=(df['mag']<(mth+1)) & (df['mag']>=mth) \n", + "ax1.scatter(np.log10(dfr['dtmc'][mask1]),np.log10(dfr['drmc'][mask1]),\n", + " color='green',label=r'$%s\\leq m<%s$'%(mth,mth+1),s=20)\n", + "mask1=(df['mag']<(mth+2)) & (df['mag']>=(mth+1)) \n", + "ax1.scatter(np.log10(dfr['dtmc'][mask1]),np.log10(dfr['drmc'][mask1]),\n", + " color='orange',label=r'$%s\\leq m<%s$'%(mth+1,mth+2),s=40)\n", + "mask1=(df['mag']<(mth+3)) & (df['mag']>=(mth+2)) \n", + "ax1.scatter(np.log10(dfr['dtmc'][mask1]),np.log10(dfr['drmc'][mask1]),\n", + " color='red',label=r'$%s\\leq m<%s$'%(mth+2,mth+3),s=60)\n", + "mask1=dfr['mag']>=(mth+3)\n", + "if mth==4:\n", + " ax1.scatter(np.log10(dfr['dtmc'][mask1]),np.log10(dfr['drmc'][mask1]),\n", + " color='m',label='$m\\geq%s$'%(mth+3),s=80)\n", + "ax1.legend()\n", + "ax1.plot();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this graph we can appreciate that the time-space area between a Main and its Origin is extremely wide. There are some earthquakes that are very close in the time scale, but they are very far from each other; on the other hand there are some events that are produced by an extremely long chains of events (1y-10y). The core of the distribution altough, is near 1hour-1month and 1Km-100Km. This study will help us on the idea of correlated and uncorrelated events." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we want to study the frequency, in terms of magnitude, of: originshocks, foreshocks, aftershocks and of all the events without distinction.\n", + "\n", + "#### Recalling notation introduced before:\n", + "\n", + "- Originshocks, all the earthquakes labeled as -1.\n", + "- Mainshocks, all the earthquakes with a magnitude above some threshold $m_{th}~$.\n", + "- Foreshocks, all the earthquakes in between Origins and Mains.\n", + "- Aftershocks, only the first generation of earthquakes produced by the Mains.\n", + "\n", + "\n", + "In particular we calculated the reverse cumulative distribution function (RCDF) for those types of earthquake, because is known that if all earthquakes are considered, then the Gutenberg-Richter Law, which itself is a RCDF, should be obtained. The Law is defined as: $$~N_{>m}=N_{TOT}10^{-bm}=10^{a-bm}~,$$ where $a$ is defined such that $N_{TOT}=10^{a}$ and b, which is the interesting parameter in terms of earthquakes behaviour, has value close to 1 in most analysis.\\\n", + "This law is telling us the relative frequency of earthquakes with magnitude equal or bigger then a given $m$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "#--------------------------------------ORIGINSHOCK---------------\n", + "#fit was done in previous cell\n", + "print('---Fit Parameters with Errors Originshock---')\n", + "print('a=%.2f'%B[0][0],'+-%.2f'%B[1][0,0]**0.5)\n", + "print('b=%.2f'%B[0][1],'+-%.2f'%B[1][1,1]**0.5)\n", + "print('r=%.4f'%resB.rvalue)\n", + "\n", + "#--------------------------------------FORESHOCK---------------\n", + "mth=5 #---CHOOSE MTH HERE-----\n", + "_,originv=bta(mth) #we call the back tree-alg\n", + "counts2, bins2,_ = ax1.hist(df['mag'][originv], bins=binss,visible=0,density=1)\n", + "mask=counts2>0\n", + "logmF=np.log10(counts2[mask])\n", + "xF=((bins2[:-1]+bins2[1:])/2)[mask]\n", + "F=optimize.curve_fit(q,xF,logmF,p0=[0.1,0.2,0]) \n", + "Q=optimize.curve_fit(f,xF,logmF,p0=[6,-1]) \n", + "resF = stats.linregress(xF, logmF)\n", + "\n", + "#Fisher Test\n", + "NFobs=np.sum((logmF-f(xF,Q[0][0],Q[0][1]))**2)-np.sum((logmF-q(xF,F[0][0],F[0][1],F[0][2]))**2)/(3-2)\n", + "DFobs=np.sum((logmF-q(xF,F[0][0],F[0][1],F[0][2]))**2)/(len(logmF)-3)\n", + "Fobs=NFobs/DFobs\n", + "cvF=(1-stats.f.cdf(Fobs,2,len(logmF)-3))*100\n", + "#if cvF is 0, then the model with 3 parameter is necessary at that value of confidence.\n", + "#EX: cvF=1% if we set the threshold at cv=0.5% then we keep the 2p-model, \n", + "#instead if we fix cv=5% we keep the 3p-model\n", + "\n", + "print('---Fit Parameters with Errors 3P Foreshock---')\n", + "print('a=%.2f'%F[0][0],'+-%.2f'%F[1][0,0]**0.5)\n", + "print('b=%.2f'%F[0][1],'+-%.2f'%F[1][1,1]**0.5)\n", + "print('c=%.2f'%F[0][2],'+-%.2f'%F[1][2,2]**0.5)\n", + "print('---Fit Parameters with Errors 2P Foreshock---')\n", + "print('a=%.2f'%Q[0][0],'+-%.2f'%Q[1][0,0]**0.5)\n", + "print('b=%.2f'%Q[0][1],'+-%.2f'%Q[1][1,1]**0.5)\n", + "print('r=%.4f'%resF.rvalue)\n", + "print('Cv=%.4f'%cvF,'% from Fisher-Test')\n", + "\n", + "#--------------------------------------AFTERSHOCK---------------\n", + "mth=5 #---CHOOSE MTH HERE-----\n", + "drf=df['mag'][df['mag']>=mth] \n", + "i1=drf.index\n", + "n=0\n", + "for i in range(len(i1)):\n", + " ind=(df['pointer'][df['pointer']==i1[i]].index).to_numpy() \n", + " n=np.insert(n,-1,ind)\n", + "n=np.delete(n,-1) #finding earthquakes that point to the chosen onescerco tutti i terremoti che puntano ai terremoti scelti.\n", + "\n", + "counts4, bins4,_ = ax1.hist(df['mag'][n], bins=binss,visible=0,density=1,cumulative=-1)\n", + "mask=counts4>0\n", + "logmA=np.log10(counts4[mask])\n", + "xA=((bins4[:-1]+bins4[1:])/2)[mask]\n", + "A=optimize.curve_fit(f,xA,logmA,p0=[6,-1]) \n", + "resA = stats.linregress(xA, logmA)\n", + "print('---Fit Parameters with Errors Aftershock---')\n", + "print('a=%.2f'%A[0][0],'+-%.2f'%A[1][0,0]**0.5)\n", + "print('b=%.2f'%A[0][1],'+-%.2f'%A[1][1,1]**0.5)\n", + "print('r=%.4f'%resA.rvalue)\n", + "\n", + "#---------------------------GUTEMBERG-RICHTER LAW---------------\n", + "counts5, bins5,_ = ax1.hist(df['mag'], bins=binss,visible=0,cumulative=-1,density=1)\n", + "mask=counts5>0\n", + "logmG=np.log10(counts5[mask]) \n", + "xG=((bins5[:-1]+bins5[1:])/2)[mask]\n", + "G=optimize.curve_fit(f,xG,logmG,p0=[6,-1]) \n", + "resG = stats.linregress(xG,logmG)\n", + "print('---Fit Parameters with Errors Gutenberg-Richter Law---')\n", + "print('a=%.2f'%G[0][0],'+-%.2f'%G[1][0,0]**0.5)\n", + "print('b=%.2f'%G[0][1],'+-%.2f'%G[1][1,1]**0.5)\n", + "print('r=%.4f'%resG.rvalue)\n", + "\n", + "#------------------------------FINAL-OUTPUT---------------------\n", + "pandas_bokeh.output_notebook()\n", + "pd.set_option('plotting.backend', 'pandas_bokeh')\n", + "a = figure(title='RCDF & Fit Magnitude',title_location = \"above\",plot_width=800, plot_height=600)\n", + "a.title.text_font_size = \"30px\"\n", + "a.xaxis.axis_label_text_font_size = '15pt'\n", + "a.yaxis.axis_label_text_font_size = '15pt'\n", + "a.yaxis.axis_label = '\\u006C\\u006F\\u0067\\u2081\\u2080[ N\\u2098 /N\\u209C\\u2092\\u209C ]'\n", + "a.xaxis.axis_label = 'Magnitude'\n", + "a.scatter(xB, logmB, color=\"orange\", legend_label=\"RCDF Originshock\", marker = 'square', size=7)\n", + "a.scatter(xA, logmA, color=\"blue\", legend_label=\"RCDF Aftershock\", marker = 'circle', size=7)\n", + "a.scatter(xF, logmF, color=\"green\", legend_label=\"RCDF Foreshock\", marker = 'triangle', size=7)\n", + "a.scatter(xG, logmG, color=\"red\", legend_label=\"RCDF GR-Law\", marker = 'diamond', size=7)\n", + "xx=np.linspace(2,7,100)\n", + "a.line(xx, f(xx,B[0][0],B[0][1]), line_width=2, legend_label=\"Fit Originshock\", color=\"orange\")\n", + "a.line(xx, f(xx,A[0][0],A[0][1]), line_width=2, legend_label=\"Fit Aftershock\", color=\"blue\")\n", + "a.line(xx, q(xx,F[0][0],F[0][1],F[0][2]), line_width=2, legend_label=\"Fit1 Foreshock\", color=\"green\")\n", + "a.line(xx, f(xx,Q[0][0],Q[0][1]), line_width=2, legend_label=\"Fit2 Foreshock\", color=\"green\",line_dash='dashed')\n", + "a.line(xx, f(xx,G[0][0],G[0][1]), line_width=2, legend_label=\"Fit GR-Law\", color=\"red\")\n", + "\n", + "a.legend.location = \"bottom_left\"\n", + "a.legend.click_policy=\"hide\"\n", + "show(a)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### OBSERVATIONS\n", + "\n", + "1. Consider the definition of the Gutenberg-Richter Law, then if we fit the Reverse CDF we expect that our $a$ should be 0. We don't have this result because we have only earthquakes with $m\\geq2$, so in some sense we've translated our fit line.\n", + "2. The most intresting thing is that in GR-Law we've used all the data, instead in the Origin and Aftershock study we've used only a small part of the set, and yet we've found that they have a lot in common with the more general GR-Law. In Aftershocks for example, we see more events for mid values of magnitude and this sounds reasonable because the Mains could induce other strong earthquakes.\n", + "3. Foreshocks don't display the same behaviour, and after performing a Fisher Test, we can say that for $m_{th}=4~$ and $~m_{th}=5$ the model is not linear, so we could use a parabolic fit, maybe because there is some $m^{2}$ contributions, or an activation function, that is constant at the beginning and then falls as a line. In general if we have to use a line, we obtain smaller b-values and this is not understand yet. (see more at https://cds.cern.ch/record/589492/files/0210130.pdf)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ANALYSIS OF THE DISTRIBUTION $P_{m}(r)$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "An important aspect is to find the distribution $P_{m}(r)$ of the spatial distance between two consecutive events (in terms of time) and considering only earthquakes of magnitude m or above. \\\n", + "The model used for the distribution is $K \\cdot r^{-df}$, where $K$ and $d_{f}$ are free parameters. The model is analogue to the Omori Law, that will be introduced and discussed following. \\\n", + "In this section we tried to estimate the value of $d_{f}$, which will be useful in further analysis, and to understand if m affects $P_{m}(r)$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def fit_model (r,p0,p1):\n", + " return p0*(r)**(-p1)\n", + "\n", + "# Pm(r) COMPUTATION \n", + "def Pmr_hist(df,m,n_bins):\n", + "# df = dataframe\n", + "# m = magnitude threshold\n", + "\n", + " # Filtered dataframe\n", + " dr = df[df[\"mag\"]>m]\n", + " \n", + "# Calculation of distaces (km) between two consecutive events\n", + " distance = []\n", + " X = np.array(dr[\"x\"].tolist())\n", + " Y = np.array(dr[\"y\"].tolist())\n", + " Z = np.array(dr[\"z\"].tolist())\n", + " n_events = len(X)\n", + "\n", + " for i in range(0,n_events-1):\n", + " d = np.sqrt( (X[i]-X[i+1])**2. + (Y[i]-Y[i+1])**2. + (Z[i]-Z[i+1])**2. )\n", + " distance.append(d/1000) \n", + " \n", + "# Making histogram with loglog scale, for visualizing data and fit\n", + " log_bins = np.logspace(0,3,n_bins+1)\n", + " histo_range=(0,1000)\n", + " \n", + " log_y_bins, log_x_bins, p = plt.hist(distance,bins=log_bins,range=histo_range,density=True)\n", + " plt.close('all')\n", + " \n", + " #bins centering\n", + " size = len(log_x_bins)\n", + " for i in range (0,size-1): log_x_bins[i]=np.sqrt( log_x_bins[i]*log_x_bins[i+1] )\n", + " #removing right edge of last bin\n", + " log_x_bins = np.delete(log_x_bins, size-1)\n", + " \n", + " # Doing fit; bin with values 0 are excluded\n", + " index = []\n", + " for i in range(n_bins):\n", + " if(log_x_bins[i]>600 or log_y_bins[i]==0): index.append(i)\n", + " log_x = np.delete(log_x_bins,index)\n", + " log_y = np.delete(log_y_bins,index)\n", + " \n", + " #bins poissonian error, necessary for the fit\n", + " y_err = [] \n", + " for yi in log_y: \n", + " y_err.append( np.sqrt(yi/n_events) )\n", + " \n", + " fit_par, fit_cov = optimize.curve_fit(fit_model, xdata=log_x, ydata=log_y, sigma=y_err, p0=[0.15,0.9])\n", + " \n", + " return log_x_bins , log_y_bins , fit_par, fit_cov" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# MAKING BOKEH PLOT\n", + "\n", + "pandas_bokeh.output_notebook()\n", + "pd.set_option('plotting.backend', 'pandas_bokeh')\n", + "\n", + "# getting data point and fit from previous function\n", + "x2, y2, fit_par2, fit_cov2 = Pmr_hist(df,2,30)\n", + "x3, y3, fit_par3, fit_cov3 = Pmr_hist(df,3,30)\n", + "x4, y4, fit_par4, fit_cov4 = Pmr_hist(df,4,30)\n", + "x5, y5, fit_par5, fit_cov5 = Pmr_hist(df,5,30)\n", + "\n", + "fit_par = (fit_par2,fit_par3,fit_par4,fit_par5)\n", + "fit_cov = (fit_cov2,fit_cov3,fit_cov4,fit_cov5)\n", + "\n", + "p = figure(y_axis_type=\"log\",x_axis_type=\"log\",title='Histograms of distance distribution'\n", + " ,title_location = \"above\",plot_width=800, plot_height=600, y_range=(10**(-6), 0.4),x_range=(0.9, 1300))\n", + "p.xaxis.axis_label = 'r (km)'\n", + "p.yaxis.axis_label = 'P\\u2098 (r)'\n", + "p.title.text_font_size = \"30px\"\n", + "p.xaxis.axis_label_text_font_size = '15pt'\n", + "p.yaxis.axis_label_text_font_size = '15pt'\n", + "# scatter points\n", + "p.scatter(x=x2, y=y2, color=\"blue\", legend_label=\"m=2\", marker = 'circle', size=6, alpha=0.75)\n", + "p.scatter(x=x3, y=y3, color=\"orange\", legend_label=\"m=3\", marker = 'square', size=6)\n", + "p.scatter(x=x4, y=y4, color=\"green\", legend_label=\"m=4\", marker = 'triangle', size=6)\n", + "p.scatter(x=x5, y=y5, color=\"red\", legend_label=\"m=5\", marker ='diamond', size=6)\n", + "\n", + "# fit functions\n", + "xf = np.linspace(1,600,1000)\n", + "p.line(xf, fit_model(xf,fit_par[0][0],fit_par[0][1]), line_width=2, legend_label=\"fit m=2\", color=\"blue\",alpha=0.75)\n", + "p.line(xf, fit_model(xf,fit_par[1][0],fit_par[1][1]), line_width=2, legend_label=\"fit m=3\", color=\"orange\")\n", + "p.line(xf, fit_model(xf,fit_par[2][0],fit_par[2][1]), line_width=2, legend_label=\"fit m=4\", color=\"green\")\n", + "p.line(xf, fit_model(xf,fit_par[3][0],fit_par[3][1]), line_width=2, legend_label=\"fit m=5\", color=\"red\")\n", + "\n", + "p.legend.location = \"top_right\"\n", + "p.legend.click_policy=\"hide\"\n", + "show(p)\n", + "\n", + "print(\"Obtained parameters:\")\n", + "print()\n", + "dfv = []\n", + "for m in range(0,4,1):\n", + " cov = fit_cov[m]\n", + " par = fit_par[m]\n", + " print(\"m = \",m+2,\": K =\", round( par[0],3 ), \"+-\", round( np.sqrt(cov[0][0]),3 ),\n", + " \" df =\", round( par[1],2 ), \"+-\", round( np.sqrt(cov[1][1]),2 ) )\n", + " dfv.append(round( par[1],2 ))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### OBSERVATIONS\n", + "\n", + "Because of the spatial limitation of the dataset, $P_{m}(r)$ can not be computed for every value of r, for this reason only r in [1,600] km were considered in fitting the data distribution (histogram). \\\n", + "In fact, as visible in the graph, for value of r > 600 km data distribution decreases very fast.\n", + "\n", + "The obtained values of K and $d_{f}$ show how $P_{m}(r)$ = $K \\cdot r^{-df}$ seems to be weakly dependent on the threshold m. It means that $P_{m}(r)$ is, in a good approximation, a law followed by all earthquakes, at least in the studied region." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ANALYSIS OF THE DISTRIBUTION $P_m(t)$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The _waiting time_ is the time interval between an event and the next one in the sequence, when we examinate these intervals of time we lose information about the real time of each event, which can be seen as a disvantage, but this depends on which relations we consider most interesting. \\\n", + "\\\n", + "In our analysis we calculate the waiting times between time-neighbouring events and make a histogram of the count using successive bins. Because of the range of times is very large, it would have been impratical to use equal-sized bins and so bins that are logarithmically increasing in length were used (the log scale is used also for the y axis). For times below 40 s, earthquakes overlap and it is not easy distinguish separate events, for this reason in the following, we'll be considered only the times for T > 40 s. \\\n", + "It is important to also analyse the data in order to try to identify underlying patterns of behaviour and to compute the distribution $P_m(t)$ of waiting times for events of magnitude m or above, for this reason we have considered three differents models to fit the data based on the Omori law. \\\n", + "\\\n", + "Fusakichi Omori proposed his famous law in 1894, which states that immediately after an earthquake the frequency of a sequence of aftershocks decays with time t as $N(t) \\propto t^{-\\alpha}$. It is relevant to recall that this law is a purely empirical without any complete underlying physical theory and it is therefore likely that it is only an approximation to the reality. The models that we consider are the followings: \\\n", + "\\\n", + "Model 1: $P_m^1 (t) = \\frac{K}{t^\\alpha}$ \\\n", + "\\\n", + "Model 2: $P_m^2 (t) = \\frac{K}{t + c}$ \\\n", + "\\\n", + "Model 3: $P_m^3 (t) = \\frac{K}{(t+c)^\\alpha}$ " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Definition of 3 different models to fit the distribution\n", + "\n", + "def fit_pdf(x, K, p):\n", + " return K*((x)**(-p))\n", + "\n", + "def fit_pdf2(x, K, c):\n", + " return K*((x+c)**(-1))\n", + "\n", + "def fit_pdf3(x, K, c, p):\n", + " return K*((x+abs(c))**(-p))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Tmin = 40 #short times are difficult to detect and distiguish\n", + "#Definition of the function to create histograms and fit them\n", + "def histo(m, nbins, Tdo, Tup):\n", + " plt.figure(figsize=(12,7))\n", + " mask = df['mag'] >= m #select only events with magnitude m or above\n", + " dt = df['t'][mask].diff().dropna(how='any') #compute the waiting times\n", + " dt = dt[dt > Tmin] #select only times above Tmin\n", + " n_events = len(dt)\n", + " \n", + " #histogram of waiting time\n", + " n, bins, patches = plt.hist(dt, bins = np.logspace(0, 8, nbins), \n", + " alpha = 0.4, color = 'r', edgecolor='b', label = 'Data distribution', density=True, \n", + " linewidth=1, zorder=1, visible=0)\n", + " plt.close('all')\n", + " \n", + " Y = n[n>0] #we don't consider empty bins for the fit\n", + " X = np.array([0.5 * (bins[i] + bins[i+1]) for i in range(len(bins)-1)])\n", + " X = X[n>0]\n", + " Yf = Y[X>Tdo]\n", + " Xf = X[X>Tdo]\n", + " \n", + " \n", + " #bins poissonian error, necessary for the fit\n", + " y_err1 = [] \n", + " for yi in Yf[Xf=magth[s(i,j)])\n", + " dr=df[mask]\n", + " dr=dr.drop(columns=['t','pointer','mag','date'])\n", + " dr['dpw']=((dr.diff()**2).sum(axis=1))**0.5 #dpw stay for distance pairs waise\n", + " dr.insert(loc=0, column='t', value=df['t'][dr.index])\n", + " dr.insert(loc=1, column='dt', value=df['t'][dr.index])\n", + " dr['dt']=dr['dt'].diff()\n", + " dr=dr.drop(dr.index[0])\n", + " dr=dr.drop(dr[dr['dpw'] > R[z,s(i,j)]].index)\n", + "\n", + " counts6, bins6, bars =ax6.hist(dr['dt'], visible=0,bins=np.logspace(-1,9,25),density=1) #-1,8 19 // -1,9 26\n", + " yhb=counts6[counts6>0]\n", + " yhb=yhb\n", + " Yu=np.insert(Yu,0,yhb)\n", + " xF=(bins6[:-1]*bins6[1:len(counts6)+1])**0.5\n", + " xF=xF[counts6>0]\n", + " Xu=np.insert(Xu,0,xF)\n", + " ax6.scatter(np.log10(xF),np.log10(yhb),color=colors[z],marker=markers[s(i,j)],\n", + " label='$m\\geq%s, R<%.0f km$'%(magth[s(i,j)],R[z,s(i,j)]/1000))\n", + "\n", + " Yt1=alpha[s(i,j)]*np.log10(xF)+np.log10(yhb)\n", + " Xt1=np.log10(xF) + b*magth[s(i,j)] + dfv[s(i,j)]*np.log10(R[z,s(i,j)]) +const\n", + " ax7.scatter(Xt1,Yt1,color=colors[z],marker=markers[s(i,j)],\n", + " label='$m\\geq%s, R<%.0f km$'%(magth[s(i,j)],R[z,s(i,j)]/1000))\n", + "\n", + "ax7.set_title('Universal Law of Earthquakes',size=20)\n", + "ax7.set_xlabel('$log[~ctS^{-b}R^{d_{f}}~]$',size=15)\n", + "ax7.set_ylabel('$log[~t^{\\\\alpha}P_{m,R}(t)~]$',size=15)\n", + "ax7.plot(np.zeros(100),np.linspace(-3.5,-1,100),'k--')\n", + "ax7.text(0.5,-1, 'Uncorrelated Phase',size=15)\n", + "ax7.text(-2.3,-1, 'Correlated Phase',size=15)\n", + "ax7.legend(loc='best', bbox_to_anchor=(1.05, 1) )\n", + "\n", + "ax6.set_xlim(-1,8)\n", + "ax6.set_ylim(-10.7,-1.5)\n", + "ax6.set_title('$P_{m,R}(t)$',size=20)\n", + "ax6.set_xlabel('$log[~t~]$',size=15)\n", + "ax6.set_ylabel('$log[~P_{m,R}(t)~]$',size=15)\n", + "ax6.legend(loc='best');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### OBSERVATIONS\n", + "The last graph shows the so called \"Universal Law of Earthquakes\", that is a pattern manteined over 8 decades. \n", + "\n", + "A critical part of this analysis was to fix some proper values for the thresholds: with a large $R$ in fact, $P_{m,R}(t)$ becomes the $P_{m}(t)$ studied before, and in some sense choosing $R$ means that we are selecting the area of interest of the study.\n", + "\n", + "As well described in the paper, the constant left part is related to a highly correlated region of power-law behavior, which corresponds to the Omori Law, in which a sequence of earthquakes are temporally correlated; the decaying right part of the graph for large $x$ is consistent with an exponential decay, indicating a random Poisson distribution of interoccurrence times and is indicative of earthquakes that would be classified as uncorrelated events and would be seen as responsible for creating the proceeding sequence of correlated earthquakes. \n", + "\n", + "To understand the Unified Law for Earthquakes, it is essential to see what the value of $x$ represents. The quantity $L^{df}S^{-b}$ in the scaling function represents the average number of earthquakes per unit time, with seismic moment greater than S (where $S=10^{-m}$) occurring inside a region defined by $R$. Therefore, $x$ is a measure of the number of earthquakes happening within a time interval T, so the statistics of earthquakes happening within minutes of an initial earthquake can be related to the statistics of earthquakes separated by 10 years.\n", + "\n", + "The data collapse supports the view that an earthquake is an SOC (Self-Organized Criticality) phenomenon, because only critical processes exhibit data\n", + "collapse, known as scaling in critical phenomena. \n", + "\n", + "#### FINAL RESULT\n", + "Thanks to our results, we have learned that, a plausible consequence of this view, is that the process that creates main shocks also produces aftershocks, and so there is only one relaxation mechanism responsible for all earthquakes.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bibliography\n", + "\n", + "This study is based on \"Unified scaling law for earthquakes\" from Kim Christensen, Leon Danon, Tim Scanlon, and Per Bak. \\\n", + "The paper can be found at the following link: https://www.pnas.org/content/pnas/99/suppl_1/2509.full.pdf" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/June_1992_detailed_earthquake.gif b/June_1992_detailed_earthquake.gif new file mode 100755 index 0000000..2540960 Binary files /dev/null and b/June_1992_detailed_earthquake.gif differ diff --git a/gif_background.png b/gif_background.png new file mode 100644 index 0000000..9ef5fb9 Binary files /dev/null and b/gif_background.png differ