Skip to content

Instantly share code, notes, and snippets.

@cossatot

cossatot/rupture_probabilities.ipynb

Created April 22, 2017 10:16
Show Gist options
  • Save cossatot/f5d301bfeb359366533e6633157d6c77 to your computer and use it in GitHub Desktop.
Save cossatot/f5d301bfeb359366533e6633157d6c77 to your computer and use it in GitHub Desktop.
Download ZIP
jupyter notebook for calculating fault rupture probabilites (or likelihoods)
Display the source blob
Display the rendered blob
Raw
rupture_probabilities.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib notebook"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import json\n",
"import numba\n",
"from numba import jit, float64\n",
"import time\n",
"import shapely\n",
"from shapely.geometry import LineString\n",
"from multiprocess import Process, Manager, Pool\n",
"from scipy import sparse\n",
"import copy\n",
"from ast import literal_eval"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Rupture probability calculations\n",
"\n",
"This is a notebook demonstrating how to take a GIS file of fault traces,\n",
"and calculate potential single- and multifault-ruptures, and assign some\n",
"probability to each rupture based on geometrical criteria. \n",
"\n",
"The results can be used for hazard calculations\n",
"or for evaluating which ruptures are worth more \n",
"computationally-expensive mechanical or dynamic modeling."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### First we load geojson data\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"with open('../../faulted_earth/c_am_car/geojson/central_am_caribbean.geojson', 'r') as f:\n",
" gj = json.load(f)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"187"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(gj['features'])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"faults = gj['features']\n",
"\n",
"coords = [feat['geometry']['coordinates'] for feat in faults]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def gc_dist(p0, p1, R=6371):\n",
" '''\n",
" Returns the (approximate) great-circle distance\n",
" between two points (in long-lat). Returned\n",
" distance in kilometers.\n",
" '''\n",
" p0 = np.radians(p0)\n",
" p1 = np.radians(p1)\n",
" \n",
" x = (p0[0]-p1[0]) * np.cos((p0[1]+p1[1]) /2)\n",
" y = p0[1] - p1[1]\n",
" \n",
" return R * np.sqrt((x**2 + y**2))\n",
"\n",
"\n",
"# this is a just-in-time compiled version using Numba\n",
"\n",
"R = 6371 # earth radius\n",
"@jit(float64[:](float64[:]))\n",
"def pairwise_dist(X):\n",
" '''\n",
" Pairwise distance function (using same calculations as gc_dist)\n",
" but accelerated for making a pairwise distance matrix for all\n",
" points in list\n",
" '''\n",
" M = X.shape[0]\n",
" N = X.shape[1]\n",
" D = np.empty((M, M), dtype=np.float)\n",
" for i in range(M):\n",
" for j in range(M):\n",
" d = ((X[i,0] - X[j,0]) * np.cos((X[i,1] - X[j,1])/2))**2\n",
" d += (X[i,1]-X[j,1])**2\n",
" D[i,j] = np.sqrt(d) * R\n",
" \n",
" return D"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compile pairwise distance function\n",
"_ = pairwise_dist(np.radians(coords[0]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Resample coordinates along each fault to 0.1° (~10 km)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Copy GeoJSON dict to new file for modification\n",
"# Only consider first 100 features for computational quickness\n",
"res_faults = copy.deepcopy(gj)\n",
"res_faults['features'] = res_faults['features'][:100]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def resample(feat, distance):\n",
" coords = feat['geometry']['coordinates']\n",
" ls = LineString(coords)\n",
" \n",
" n_pts = int(round(ls.length / distance))\n",
" if n_pts == 0:\n",
" n_pts = 1\n",
" nls = LineString([ls.interpolate(float(n)/n_pts, normalized=True)\n",
" for n in range(n_pts + 1)])\n",
" \n",
" newfeat = copy.deepcopy(feat)\n",
" newfeat['geometry']['coordinates'] = list(nls.coords)\n",
" \n",
" return newfeat"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"res_faults['features'] = [resample(feat, 0.1) for feat in res_faults['features']]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# write results to geojson to inspect\n",
"\n",
"with open('../../faulted_earth/c_am_car/geojson/'\n",
" +'central_am_caribbean_resampled.geojson', 'w') as f:\n",
" json.dump(res_faults, f)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"492"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"num_nodes = sum(len(f['geometry']['coordinates']) for f in res_faults['features'])\n",
"num_nodes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Node DataFrame\n",
"\n",
"Make dataframe of nodes on each fault, with id (the feature id from the GIS file)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def fault_to_df(fault):\n",
" df = pd.DataFrame(data=np.array(fault['geometry']['coordinates']),\n",
" columns=['lon', 'lat'])\n",
" df['fid'] = fault['properties']['ogc_fid'] #fault ID\n",
" return df"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"node_df = pd.concat((fault_to_df(fault) for fault in res_faults['features']),\n",
" ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(492, 3)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"node_df.shape"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>lon</th>\n",
" <th>lat</th>\n",
" <th>fid</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-93.647518</td>\n",
" <td>16.983588</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-93.549913</td>\n",
" <td>16.964436</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-93.455689</td>\n",
" <td>16.932247</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-93.360046</td>\n",
" <td>16.904757</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-93.263945</td>\n",
" <td>16.879527</td>\n",
" <td>1</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" lon lat fid\n",
"0 -93.647518 16.983588 1\n",
"1 -93.549913 16.964436 1\n",
"2 -93.455689 16.932247 1\n",
"3 -93.360046 16.904757 1\n",
"4 -93.263945 16.879527 1"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"node_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>lon</th>\n",
" <th>lat</th>\n",
" <th>fid</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>487</th>\n",
" <td>-86.702478</td>\n",
" <td>12.521014</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>488</th>\n",
" <td>-86.615644</td>\n",
" <td>12.465571</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>489</th>\n",
" <td>-86.528504</td>\n",
" <td>12.410605</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>490</th>\n",
" <td>-87.647444</td>\n",
" <td>12.998417</td>\n",
" <td>100</td>\n",
" </tr>\n",
" <tr>\n",
" <th>491</th>\n",
" <td>-87.587523</td>\n",
" <td>12.943162</td>\n",
" <td>100</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" lon lat fid\n",
"487 -86.702478 12.521014 99\n",
"488 -86.615644 12.465571 99\n",
"489 -86.528504 12.410605 99\n",
"490 -87.647444 12.998417 100\n",
"491 -87.587523 12.943162 100"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"node_df.tail()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Matrix of distances between each node\n",
"node_dists = pd.DataFrame(index=node_df.index, columns=node_df.index,\n",
" data=pairwise_dist(np.radians(node_df[['lon', 'lat']].values)))"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# List of fault IDs\n",
"fids = [f['properties']['ogc_fid'] for f in res_faults['features']]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Segment DataFrame\n",
"\n",
"This is a dataframe with segments (the section of fault in between\n",
"each pair of nodes). It will have attributes such as strike, dip and rake\n",
"that the nodes themselves do not have."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Make dataframe of all fault segments (i.e. the fault section in between two nodes)\n",
"\n",
"def make_segs(fid, df=node_df):\n",
" fdf = df[df.fid==fid]\n",
" \n",
" inds = list(fdf.index)\n",
" starts = inds[:-1]\n",
" stops = inds[1:]\n",
" \n",
" pairs = list(zip(starts, stops))\n",
" \n",
" lens = [gc_dist(node_df.iloc[s0][['lon', 'lat']].values, \n",
" node_df.iloc[s1][['lon', 'lat']].values) \n",
" for (s0, s1) in pairs]\n",
" \n",
" seg_df = pd.DataFrame(columns=['start', 'stop'],\n",
" data=pairs)\n",
" seg_df['fid'] = fid\n",
" seg_df['seg_length'] = lens\n",
" return seg_df"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"seg_df = pd.concat((make_segs(fid) for fid in fids), ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def calc_strike(seg):\n",
"\n",
" p0 = node_df.ix[seg.start, ['lon', 'lat']]\n",
" p1 = node_df.ix[ seg.stop, ['lon', 'lat']]\n",
" \n",
" lon0, lat0 = np.radians(p0[['lon', 'lat']])\n",
" lon1, lat1 = np.radians(p1[['lon', 'lat']])\n",
" \n",
" y = np.sin(lon1-lon0) * np.cos(lat1)\n",
" x = np.cos(lat0) * np.sin(lat1) - np.sin(lat0) * np.cos(lat1) * np.cos(lon1-lon0)\n",
" \n",
" angle = np.degrees(np.arctan2(y,x))\n",
" \n",
" az = 90 - angle #-(angle - 90)\n",
" \n",
" while az < 0:\n",
" az += 360\n",
" while az > 360:\n",
" az -= 360\n",
" \n",
" return az\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"seg_df['strike'] = seg_df.apply(calc_strike, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>start</th>\n",
" <th>stop</th>\n",
" <th>fid</th>\n",
" <th>seg_length</th>\n",
" <th>strike</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>10.596599</td>\n",
" <td>348.420283</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" <td>10.642168</td>\n",
" <td>340.360652</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>3</td>\n",
" <td>1</td>\n",
" <td>10.624013</td>\n",
" <td>343.292410</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>1</td>\n",
" <td>10.602795</td>\n",
" <td>344.670753</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>5</td>\n",
" <td>1</td>\n",
" <td>10.525126</td>\n",
" <td>334.825527</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop fid seg_length strike\n",
"0 0 1 1 10.596599 348.420283\n",
"1 1 2 1 10.642168 340.360652\n",
"2 2 3 1 10.624013 343.292410\n",
"3 3 4 1 10.602795 344.670753\n",
"4 4 5 1 10.525126 334.825527"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"seg_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"100"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"seg_df.fid.max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get fault dip, dip dir, rake, and slip type from GeoJSON"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"{'Sinistral',\n",
" 'Normal-Sinistral',\n",
" 'Normal',\n",
" 'Thrust',\n",
" 'Sinistral-Normal',\n",
" 'Dextral',\n",
" None,\n",
" 'Dextral Normal'}"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"set(f['properties']['slip_type'] for f in res_faults['features'])"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dip_d = {'Normal-Sinistral' : 70.,\n",
" 'Sinistral' : 90.,\n",
" 'Dextral Normal' : 70.,\n",
" 'Normal' : 50.,\n",
" 'Thrust' : 25.,\n",
" 'Dextral' : 90.,\n",
" # None, not sure how to handle this\n",
" 'Sinistral-Normal' : 70.}"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def get_dip(f):\n",
" \n",
" props = gj['features'][int(f.fid)]['properties']\n",
" \n",
" dip = props['ns_average_dip']\n",
" if dip is not None:\n",
" if type(dip) == str:\n",
" try:\n",
" dip = literal_eval(dip)\n",
" except SyntaxError: # frequently '(90,,)' etc.\n",
" dip = dip.strip('()')\n",
" dip = literal_eval(dip.split(',')[0])\n",
" \n",
" if np.isscalar(dip):\n",
" dip = float(dip)\n",
" else:\n",
" dip = float(dip[0])\n",
" \n",
" else:\n",
" try:\n",
" dip = dip_d[props['slip_type']]\n",
" except KeyError: # in case slip_type == None:\n",
" dip = 60. # whatever\n",
" \n",
" return dip\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"seg_df['dip'] = seg_df.apply(get_dip, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rake_d = {'Normal-Sinistral' : -45.,\n",
" 'Sinistral' : 0.,\n",
" 'Dextral Normal' : -135.,\n",
" 'Normal' : -90.,\n",
" 'Thrust' : 90.,\n",
" 'Dextral' : 180.,\n",
" # None, not sure how to handle this\n",
" 'Sinistral-Normal' : -45.}"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def get_rake(f):\n",
" \n",
" props = gj['features'][int(f.fid)]['properties']\n",
" \n",
" rake = props['ns_average_rake']\n",
" if rake is not None:\n",
" if type(rake) == str:\n",
" try:\n",
" rake = literal_eval(rake)\n",
" except SyntaxError: # frequently '(90,,)' etc.\n",
" rake = rake.strip('()')\n",
" rake = literal_eval(rake.split(',')[0])\n",
" \n",
" if np.isscalar(rake):\n",
" rake = float(rake)\n",
" else:\n",
" rake = float(rake[0])\n",
" \n",
" else:\n",
" try:\n",
" rake = rake_d[props['slip_type']]\n",
" except KeyError: # in case slip_type == None:\n",
" rake = 60. # whatever\n",
" \n",
" return rake"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"seg_df['rake'] = seg_df.apply(get_rake, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Do right-hand-rule dip thing\n",
"# Do rake"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>start</th>\n",
" <th>stop</th>\n",
" <th>fid</th>\n",
" <th>seg_length</th>\n",
" <th>strike</th>\n",
" <th>dip</th>\n",
" <th>rake</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>387</th>\n",
" <td>485</td>\n",
" <td>486</td>\n",
" <td>99</td>\n",
" <td>11.258960</td>\n",
" <td>327.193535</td>\n",
" <td>80.0</td>\n",
" <td>-135.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>388</th>\n",
" <td>486</td>\n",
" <td>487</td>\n",
" <td>99</td>\n",
" <td>11.263397</td>\n",
" <td>326.562061</td>\n",
" <td>80.0</td>\n",
" <td>-135.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>389</th>\n",
" <td>487</td>\n",
" <td>488</td>\n",
" <td>99</td>\n",
" <td>11.263858</td>\n",
" <td>326.825157</td>\n",
" <td>80.0</td>\n",
" <td>-135.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>390</th>\n",
" <td>488</td>\n",
" <td>489</td>\n",
" <td>99</td>\n",
" <td>11.264340</td>\n",
" <td>327.149775</td>\n",
" <td>80.0</td>\n",
" <td>-135.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>391</th>\n",
" <td>490</td>\n",
" <td>491</td>\n",
" <td>100</td>\n",
" <td>8.939188</td>\n",
" <td>316.587898</td>\n",
" <td>80.0</td>\n",
" <td>-45.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop fid seg_length strike dip rake\n",
"387 485 486 99 11.258960 327.193535 80.0 -135.0\n",
"388 486 487 99 11.263397 326.562061 80.0 -135.0\n",
"389 487 488 99 11.263858 326.825157 80.0 -135.0\n",
"390 488 489 99 11.264340 327.149775 80.0 -135.0\n",
"391 490 491 100 8.939188 316.587898 80.0 -45.0"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"seg_df.tail()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define floating ruptures\n",
"\n",
"Floating ruptures are defined as any possible sequence of contiguous segments\n",
"on the fault. For a fault with $n$ segments, there are $\\frac{n(n+1)}{2}$ possible\n",
"floating ruptures. For a fault with 4 segments, this is $\\frac{4\\cdot5}{2}$ or 10 ruptures.\n",
"\n",
"```\n",
" a,b,c,d == fault\n",
" -------\n",
"1 a,b,c,d | = full rupture\n",
"2 a,b,c + = 3 segment ruptures\n",
"3 b,c,d +\n",
"4 a,b |\n",
"5 b,c | = 2 segment ruptures\n",
"5 c,d |\n",
"7 a + \n",
"8 b + = 1 segment rupture\n",
"9 c +\n",
"10 d +\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def subsegs(seg_list, n):\n",
" '''\n",
" Returns all contiguous segments of length n\n",
" from segment list\n",
" '''\n",
" N = len(seg_list)\n",
" \n",
" if n == N-1:\n",
" return [((s)) for s in seg_list]\n",
" else:\n",
" return [seg_list[i : N-n+i] for i in range(n+1)]\n",
"\n",
" \n",
"def segs(seg_list):\n",
" '''\n",
" Returns all possible contiguous segments\n",
" from segment list.\n",
" '''\n",
" \n",
" segz = [seg_list]\n",
" n = len(seg_list)\n",
" \n",
" subsegz = [subsegs(seg_list, i) for i in range(n)]\n",
" \n",
" return [s for sub in subsegz for s in sub]"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Make dataframe of all possible continuous floating ruptures on the faults\n",
"\n",
"def make_rupts(fid, df=seg_df):\n",
" fdf = df[df.fid == fid]\n",
" \n",
" inds = tuple(fdf.index)\n",
" \n",
" rupts = segs(inds)\n",
" nrups = len(rupts)\n",
" \n",
" rupt_df = pd.DataFrame(pd.Series(rupts), columns=['segment_ids'])\n",
" rupt_df['fid'] = fid\n",
" \n",
" return rupt_df"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rupture_df = pd.concat((make_rupts(fid) for fid in fids), ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(1766, 2)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rupture_df.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ruptures contain tuples with all the segment IDs, as well as the fault IDs."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>segment_ids</th>\n",
" <th>fid</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>1746</th>\n",
" <td>(386, 387, 388, 389, 390)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1747</th>\n",
" <td>(385, 386, 387, 388)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1748</th>\n",
" <td>(386, 387, 388, 389)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1749</th>\n",
" <td>(387, 388, 389, 390)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1750</th>\n",
" <td>(385, 386, 387)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1751</th>\n",
" <td>(386, 387, 388)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1752</th>\n",
" <td>(387, 388, 389)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1753</th>\n",
" <td>(388, 389, 390)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1754</th>\n",
" <td>(385, 386)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1755</th>\n",
" <td>(386, 387)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1756</th>\n",
" <td>(387, 388)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1757</th>\n",
" <td>(388, 389)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1758</th>\n",
" <td>(389, 390)</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1759</th>\n",
" <td>385</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1760</th>\n",
" <td>386</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1761</th>\n",
" <td>387</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1762</th>\n",
" <td>388</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1763</th>\n",
" <td>389</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1764</th>\n",
" <td>390</td>\n",
" <td>99</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1765</th>\n",
" <td>391</td>\n",
" <td>100</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" segment_ids fid\n",
"1746 (386, 387, 388, 389, 390) 99\n",
"1747 (385, 386, 387, 388) 99\n",
"1748 (386, 387, 388, 389) 99\n",
"1749 (387, 388, 389, 390) 99\n",
"1750 (385, 386, 387) 99\n",
"1751 (386, 387, 388) 99\n",
"1752 (387, 388, 389) 99\n",
"1753 (388, 389, 390) 99\n",
"1754 (385, 386) 99\n",
"1755 (386, 387) 99\n",
"1756 (387, 388) 99\n",
"1757 (388, 389) 99\n",
"1758 (389, 390) 99\n",
"1759 385 99\n",
"1760 386 99\n",
"1761 387 99\n",
"1762 388 99\n",
"1763 389 99\n",
"1764 390 99\n",
"1765 391 100"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rupture_df.tail(20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define probability for each rupture based on geometric characteristics\n",
"\n",
"Geometric characteristics that define floating rupture probability\n",
"potentially include changes in the strike, dip and rake of ruptures.\n",
"At this point in the fault database, the strike is the only characteristic\n",
"that is different for each segment of a fault: Dip and rake are defined\n",
"per fault, while strike is calculated directly from the node geometry.\n",
"\n",
"\n",
"We define the probability of rupture $r$ based on changes in strike as \n",
"$$ p_{strike}(r)= \\prod_{i=1}\\cos(strike_i - strike_{i+1}) $$\n",
"where $i$ and $i+1$ are adjacent fault segments in $r$.\n"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def rupt_strike_probs(rupture):\n",
" strikes = seg_df.ix[rupture.segment_ids, 'strike']\n",
" if np.isscalar(strikes):\n",
" probs = 1.\n",
" else:\n",
" strike_rads = np.radians(strikes)\n",
" probs = np.product(np.cos(np.diff(strike_rads)))\n",
" return probs\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rupture_df['strike_probs'] = rupture_df.apply(rupt_strike_probs, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"count 1766.000000\n",
"mean 0.926696\n",
"std 0.115041\n",
"min 0.390364\n",
"25% 0.910718\n",
"50% 0.982654\n",
"75% 0.999679\n",
"max 1.000000\n",
"Name: strike_probs, dtype: float64"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rupture_df.strike_probs.describe()"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mean_strike(rupture):\n",
" strikes = seg_df.ix[rupture.segment_ids, 'strike']\n",
" if np.isscalar(strikes):\n",
" return strikes\n",
" else:\n",
" return np.mean(strikes)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rupture_df['mean_strike'] = rupture_df.apply(mean_strike, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>segment_ids</th>\n",
" <th>fid</th>\n",
" <th>strike_probs</th>\n",
" <th>mean_strike</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>1761</th>\n",
" <td>387</td>\n",
" <td>99</td>\n",
" <td>1.0</td>\n",
" <td>327.193535</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1762</th>\n",
" <td>388</td>\n",
" <td>99</td>\n",
" <td>1.0</td>\n",
" <td>326.562061</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1763</th>\n",
" <td>389</td>\n",
" <td>99</td>\n",
" <td>1.0</td>\n",
" <td>326.825157</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1764</th>\n",
" <td>390</td>\n",
" <td>99</td>\n",
" <td>1.0</td>\n",
" <td>327.149775</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1765</th>\n",
" <td>391</td>\n",
" <td>100</td>\n",
" <td>1.0</td>\n",
" <td>316.587898</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" segment_ids fid strike_probs mean_strike\n",
"1761 387 99 1.0 327.193535\n",
"1762 388 99 1.0 326.562061\n",
"1763 389 99 1.0 326.825157\n",
"1764 390 99 1.0 327.149775\n",
"1765 391 100 1.0 316.587898"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rupture_df.tail()"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>start</th>\n",
" <th>stop</th>\n",
" <th>fid</th>\n",
" <th>seg_length</th>\n",
" <th>strike</th>\n",
" <th>dip</th>\n",
" <th>rake</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>10.596599</td>\n",
" <td>348.420283</td>\n",
" <td>90.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" <td>10.642168</td>\n",
" <td>340.360652</td>\n",
" <td>90.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>3</td>\n",
" <td>1</td>\n",
" <td>10.624013</td>\n",
" <td>343.292410</td>\n",
" <td>90.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>1</td>\n",
" <td>10.602795</td>\n",
" <td>344.670753</td>\n",
" <td>90.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>5</td>\n",
" <td>1</td>\n",
" <td>10.525126</td>\n",
" <td>334.825527</td>\n",
" <td>90.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop fid seg_length strike dip rake\n",
"0 0 1 1 10.596599 348.420283 90.0 0.0\n",
"1 1 2 1 10.642168 340.360652 90.0 0.0\n",
"2 2 3 1 10.624013 343.292410 90.0 0.0\n",
"3 3 4 1 10.602795 344.670753 90.0 0.0\n",
"4 4 5 1 10.525126 334.825527 90.0 0.0"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"seg_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#rupture_df['dip'] = \n",
"#seg_df.ix[rupture_df.fid, 'dip']\n",
"\n",
"def get_seg_dip(rupture):\n",
" dips = seg_df[seg_df.fid == rupture.fid]['dip']\n",
" if np.isscalar(dips):\n",
" return dips\n",
" else:\n",
" return dips.values[0]\n",
" \n",
"def get_seg_rake(rupture):\n",
" rakes = seg_df[seg_df.fid == rupture.fid]['rake']\n",
" if np.isscalar(rakes):\n",
" return rakes\n",
" else:\n",
" return rakes.values[0]"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rupture_df['rake'] = rupture_df.apply(get_seg_rake, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rupture_df['dip'] = rupture_df.apply(get_seg_dip, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>segment_ids</th>\n",
" <th>fid</th>\n",
" <th>strike_probs</th>\n",
" <th>mean_strike</th>\n",
" <th>rake</th>\n",
" <th>dip</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...</td>\n",
" <td>1</td>\n",
" <td>0.698518</td>\n",
" <td>332.808503</td>\n",
" <td>0.0</td>\n",
" <td>90.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...</td>\n",
" <td>1</td>\n",
" <td>0.717178</td>\n",
" <td>331.669365</td>\n",
" <td>0.0</td>\n",
" <td>90.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14...</td>\n",
" <td>1</td>\n",
" <td>0.705486</td>\n",
" <td>332.098877</td>\n",
" <td>0.0</td>\n",
" <td>90.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...</td>\n",
" <td>1</td>\n",
" <td>0.742712</td>\n",
" <td>331.045487</td>\n",
" <td>0.0</td>\n",
" <td>90.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14...</td>\n",
" <td>1</td>\n",
" <td>0.724333</td>\n",
" <td>330.871703</td>\n",
" <td>0.0</td>\n",
" <td>90.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" segment_ids fid strike_probs \\\n",
"0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 0.698518 \n",
"1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 0.717178 \n",
"2 (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... 1 0.705486 \n",
"3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 0.742712 \n",
"4 (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... 1 0.724333 \n",
"\n",
" mean_strike rake dip \n",
"0 332.808503 0.0 90.0 \n",
"1 331.669365 0.0 90.0 \n",
"2 332.098877 0.0 90.0 \n",
"3 331.045487 0.0 90.0 \n",
"4 330.871703 0.0 90.0 "
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rupture_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def rupture_length(rupture):\n",
" seg_lengths = seg_df.ix[rupture.segment_ids, 'seg_length']\n",
" \n",
" return np.sum(seg_lengths)"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rupture_df['rup_length'] = rupture_df.apply(rupture_length, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1766"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(rupture_df)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multi-fault ruptures\n",
"\n",
"Multi-fault ruptures are those that rupture multiple segments on different\n",
"faults. This is currently defined as any combination of two single-fault\n",
"ruptures, though it can easily be expanded to any number of ruptures once\n",
"the feasible 2-fault ruptures are selected from the millions of possible\n",
"2 fault ruptures."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate possible multi-fault ruptures\n",
"\n",
"The first step is to calculate the minimum distance between all possible\n",
"2-fault ruptures. This is used to select the feasible multi-fault ruptures."
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Calculate minimum distance between nodes in each rupture\n",
"\n",
"def min_rup_dist(rup0, rup1):\n",
" if rup0.fid == rup1.fid:\n",
" return -1. \n",
" segs0 = rup0.segment_ids\n",
" segs1 = rup1.segment_ids\n",
" \n",
" if np.isscalar(segs0):\n",
" if np.isscalar(segs1):\n",
" return node_dists.ix[segs0, segs1]\n",
" else:\n",
" return min(node_dists.ix[segs0, j] for j in segs1)\n",
" \n",
" elif np.isscalar(segs1):\n",
" return min(node_dists.ix[i, segs1] for i in segs0)\n",
" \n",
" else: \n",
" return min( node_dists.ix[i, j] for i in segs0 for j in segs1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The maximum allowable distance is set at 30 km."
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"shared_d = {}\n",
"\n",
"def calc_rup_dists(ij, d=shared_d, max_dist=30.):\n",
" i, j = ij\n",
" \n",
" if i == j:\n",
" return False\n",
" \n",
" if i != j:\n",
" #rup_dists[i,j] = 0.\n",
" #return\n",
" #return False\n",
" rup0 = rupture_df.ix[i]\n",
" rup1 = rupture_df.ix[j]\n",
" \n",
" md = min_rup_dist(rup0, rup1) \n",
" if md <= max_dist:\n",
" if md > 0.:\n",
" #print(md)\n",
" d[(i, j)] = md\n",
" return True\n",
" else: return False\n",
" else:\n",
" return False"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we do the calculations. Because there are millions to do,\n",
"we'll parallelize it using the `multiprocess` module (this is\n",
"a fork from the standard `multiprocessing` library that allows\n",
"for better shared-memory data structures)."
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"54525 changes, done in 269 s\n"
]
}
],
"source": [
"manager = Manager() # Multiprocess manager\n",
"\n",
"shared_d = manager.dict() # shared-memory dictionary\n",
"\n",
"t0 = time.time()\n",
"\n",
"p = Pool(7)\n",
"\n",
"ijs = ((i, j) for i in rupture_df.index[:] # Calculate upper-diagonal distance matrix\n",
" for j in rupture_df.index[i:])# indices\n",
"\n",
"z = sum(p.imap_unordered(calc_rup_dists, ijs, chunksize=100)) #do calcs; sum=kept ruptures\n",
"\n",
"t1 = time.time()\n",
"print('{0} changes, done in {1} s'.format(z, int(t1-t0)))"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(shared_d.keys()) # checking to see that this is the same as above"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"mult_rups = shared_d.items() # make list of results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make multi-rupture DataFrame"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def make_multifault_df(item):\n",
" key, val = item\n",
" r1, r2 = key\n",
" \n",
" return pd.Series((r1, r2, val), index=['r1', 'r2', 'dist_km'])\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {
"collapsed": false
},
"outputs": [
{
"ename": "ValueError",
"evalue": "No objects to concatenate",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-60-e1ceefdd8880>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m multifault_rupture_df = pd.concat((make_multifault_df(item) for item in mult_rups),\n\u001b[1;32m 2\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m ignore_index=True).transpose()\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mmultifault_rupture_df\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'r1'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r2'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmultifault_rupture_df\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'r1'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r2'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/itchy/src/anaconda/lib/python3.5/site-packages/pandas/tools/merge.py\u001b[0m in \u001b[0;36mconcat\u001b[0;34m(objs, axis, join, join_axes, ignore_index, keys, levels, names, verify_integrity, copy)\u001b[0m\n\u001b[1;32m 843\u001b[0m \u001b[0mkeys\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevels\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlevels\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnames\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnames\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 844\u001b[0m \u001b[0mverify_integrity\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverify_integrity\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 845\u001b[0;31m copy=copy)\n\u001b[0m\u001b[1;32m 846\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mop\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_result\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 847\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/itchy/src/anaconda/lib/python3.5/site-packages/pandas/tools/merge.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, objs, axis, join, join_axes, keys, levels, names, ignore_index, verify_integrity, copy)\u001b[0m\n\u001b[1;32m 876\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 877\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobjs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 878\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'No objects to concatenate'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 879\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 880\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkeys\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: No objects to concatenate"
]
}
],
"source": [
"multifault_rupture_df = pd.concat((make_multifault_df(item) for item in mult_rups),\n",
" axis=1,\n",
" ignore_index=True).transpose()\n",
"multifault_rupture_df[['r1', 'r2']] = multifault_rupture_df[['r1', 'r2']].astype(int)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df['fid1'] = rupture_df.ix[multifault_rupture_df['r1'], 'fid'].values\n",
"multifault_rupture_df['fid2'] = rupture_df.ix[multifault_rupture_df['r2'], 'fid'].values"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`r1` and `r2` are the indices of the single-fault ruptures that are a part of this multi-segment fault.\n",
"\n",
"`fid1` and `fid2` are the indices of the faults themselves."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.dist_km.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate multi-rupture probability based on gap distance\n",
"\n",
"*Shaw and Dieterich [2007]* studied multi-fault ruptures using simulations\n",
"and found that the likelihood that a rupture jumps a gap can be\n",
"approximated by the formula \n",
"$$ p(jump) = \\exp(-d / d_0) $$\n",
"where $d$ is the distance and $d_0$ is a characteristic distance.\n",
"\n",
"Based on an unofficial survey of multi-fault ruptures, I have set\n",
"$d_0=5$ so that it is essentially 0 at 30 km. There is about a\n",
"1% probability at 20 km:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.plot(np.arange(31), np.exp(-np.arange(31)/5))\n",
"plt.xlabel('distance (km)')\n",
"plt.ylabel('probability')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df['dist_prob'] = np.exp(-multifault_rupture_df.dist_km/5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.dist_prob.describe()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multifault rupture probability based on strike, dip and rake compatibility\n",
"\n",
"There are several geometric characteristics of the ruptures that we can use\n",
"to evaluate the likelihood of each potential multifault rupture, such as the\n",
"compatibilities between strike, dip and rake of each of the faults in the \n",
"\n",
"It's important to note that a mechanical (rather than geometric) characterization\n",
"of these ruptures would take the place of this analysis. This is quick, however,\n",
"and can serve as a basic example."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"'''\n",
"Some basic equations to check geometric compatibility.\n",
"\n",
"Note that the angle differences may need to be checked for\n",
"bugs, i.e. differences in strike angles should not be\n",
"greater than 90 (strike is not a vector), while\n",
"a difference of 180 in rake is more meaningful.\n",
"\n",
"This is only an example.\n",
"'''\n",
"\n",
"\n",
"def strike_compatibility(multirupt):\n",
"\n",
" strike1 = rupture_df.ix[multirupt.r1, 'mean_strike']\n",
" strike2 = rupture_df.ix[multirupt.r2, 'mean_strike']\n",
" \n",
" return np.abs(np.cos(np.radians(strike2-strike1)))\n",
" \n",
" \n",
"def rake_compatibility(multirupt):\n",
"\n",
" rake1 = rupture_df.ix[multirupt.r1, 'rake']\n",
" rake2 = rupture_df.ix[multirupt.r2, 'rake']\n",
" \n",
" return np.abs(np.cos(np.radians(rake2-rake1)))\n",
"\n",
"\n",
"def dip_compatibility(multirupt):\n",
"\n",
" dip1 = rupture_df.ix[multirupt.r1, 'dip']\n",
" dip2 = rupture_df.ix[multirupt.r2, 'dip']\n",
" \n",
" return np.abs(np.cos(np.radians(dip2-dip1)))\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"multifault_rupture_df['strike_probs'] = multifault_rupture_df.apply(strike_compatibility,\n",
" axis=1)\n",
"multifault_rupture_df['dip_probs'] = multifault_rupture_df.apply(dip_compatibility,\n",
" axis=1)\n",
"multifault_rupture_df['rake_probs'] = multifault_rupture_df.apply(rake_compatibility,\n",
" axis=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculating final probabilities for multifault ruptures\n",
"\n",
"This is done by multiplying the probabilities for each category\n",
"across each rupture."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"multifault_rupture_df['final_probs'] = ( multifault_rupture_df.dist_prob\n",
" * multifault_rupture_df.strike_probs\n",
" * multifault_rupture_df.dip_probs\n",
" * multifault_rupture_df.rake_probs)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"multifault_rupture_df.final_probs.describe()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Number of ruptures with a probability of > 1%, out of total possible ruptures\n",
"\n",
"n_01 = sum((multifault_rupture_df.final_probs > 0.01)) \n",
"n_poss = int((len(rupture_df)**2 - len(rupture_df))/2)\n",
"\n",
"print('{0} possible ruptures out of {1} ({2:1f}%)'.format(n_01, n_poss, (100 * n_01 /n_poss)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Results\n",
"\n",
"These results illustrate how it is possible to calculate rupture probabilities\n",
"for single and multiple faults based on geometric characteristics.\n",
"\n",
"The results are not final or comprehensive, but lay out a computationally-efficient\n",
"framework for evaluating a huge number of possible ruptures without having to\n",
"compute expensive mechanical (much less dynamic) ruptures.\n",
"\n",
"This framework can be easily modified, and is very useful primarily for limiting\n",
"the set of all possible ruptures down to a much smaller set that is more possible\n",
"to model dynamically."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Future work\n",
"\n",
"The next step is to perform mechanical modeling of ruptures.\n",
"\n",
"The general plan is to calculate the shear and normal stress change\n",
"on each fault segment from a unit dislocation on each other segment.\n",
"These unit dislocations can then be scaled for slip (based on rupture\n",
"length) and summed for each rupture. This can be used to both\n",
"evaluate ruptures on a single fault, and multi-fault rupture."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$ p(r) = p_{strike}(r) \\dot p_{dip}(r) \\dot \\ldots \\dot p_{azimuth}(r) $"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment