In [1]:
import pandas as pd
import numpy as np
import json
import sys
In [2]:
sys.path.append('../scripts/')
import utils
In [2]:
#with open('../gis/gurla_z_pts.geojson') as f:
# gp = json.load(f)
In [3]:
#for feat in gp['features']:
# coords = feat['geometry']['coordinates']
# z_val = feat['properties']['Z']
# coords.append(z_val)
In [4]:
#with open('../gis/gurla_z_pts.geojson', 'w') as f:
# json.dump(gp, f)
In [5]:
#gp['features'][0]
In [3]:
pdf = pd.read_csv('../data/gis/gurla_z_pts.csv', index_col=3)
pdf.head()
Out[3]:
In [4]:
with open('../data/fault_data/gurla_tris.geojson') as f:
gt = json.load(f)
In [5]:
gt['features'][0]['geometry']['coordinates']#[0][0]
Out[5]:
In [8]:
gc = gt.copy()
In [8]:
for feat in gc['features']:
coors = feat['geometry']['coordinates'][0]
zs = pdf.loc[(int(feat['properties']['POINTA']), int(feat['properties']['POINTB']),
int(feat['properties']['POINTC']), int(feat['properties']['POINTA'])), 'Z'
].values
for i, coor in enumerate(coors):
#coor.append(float(zs[i])) # only do this the first time
pass
In [9]:
gc['features'][-1]
Out[9]:
$S = \frac{|AB \times AC|}{2}$
In [12]:
def area(coords):
ab = np.array(coords[0]) - np.array(coords[1])
ac = np.array(coords[0]) - np.array(coords[2])
return np.linalg.norm(np.cross(ab, ac)) / 2e6
In [13]:
def center(coords):
return ((np.array(coords[0]) + np.array(coords[1])
+ np.array(coords[2])) / 3).tolist()
center(coors)
Out[13]:
In [14]:
for feat in gc['features']:
coords = feat['geometry']['coordinates'][0]
feat['properties']['area_sq_km'] = area(coords)
feat['properties']['center'] = center(coords)
In [15]:
gc['features'][100]
Out[15]:
In [6]:
import halfspace.projections as hsp
In [9]:
for feat in gc['features']:
coords = feat['geometry']['coordinates'][0]
coo = [[co[0], co[1], co[2]] for co in coords[:-1]]
s, d = hsp.strike_dip_from_3_xyz(coo[0], coo[1], coo[2])
feat['properties']['strike'], feat['properties']['dip'] = s, d
In [10]:
gc['features'][100]
Out[10]:
In [11]:
strikes = [feat['properties']['strike'] for feat in gc['features']]
In [12]:
utils.add_strike_dip(gc['features'])
utils.add_rake_from_trend(gc['features'], 277.5, 10)
In [13]:
strikes = [feat['properties']['strike'] for feat in gc['features']]
dips = [feat['properties']['dip'] for feat in gc['features']]
rakes = [feat['properties']['rake'] for feat in gc['features']]
cx = [feat['properties']['center'][0] for feat in gc['features']]
cy = [feat['properties']['center'][1] for feat in gc['features']]
cz = [feat['properties']['center'][2] for feat in gc['features']]
In [14]:
def trend_from_sd_rake(strike=0., dip=0., rake=0., angle='degrees',
aki_richards=False):
if angle == 'degrees':
strike = np.radians(strike)
dip = np.radians(dip)
rake = np.radians(rake)
trend = -1 * (np.arctan( np.cos(dip) * np.tan(rake)) ) + strike
if aki_richards == True:
if np.isscalar(trend):
if rake > np.pi/2:
trend -= np.pi
elif rake < - np.pi/2:
trend += np.pi
else:
trend[rake > np.pi/2] -= np.pi
trend[rake < -np.pi/2] += np.pi
if angle == 'degrees':
trend = np.degrees(trend)# + 180
trend = hsp.unwrap_angle(trend)
return trend
In [15]:
trends = trend_from_sd_rake(strikes, dips, rakes, aki_richards=True)
In [16]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
In [17]:
plt.figure(figsize=(12, 12))
plt.scatter(cx, cy, c=strikes,
s=50, lw=0)
plt.colorbar()
plt.axis('equal')
plt.show()
In [18]:
plt.figure(figsize=(12, 12))
plt.scatter(cx, cy, c=rakes,
s=50, lw=0)
plt.colorbar()
plt.axis('equal')
plt.show()
In [19]:
plt.figure(figsize=(12, 12))
plt.scatter(cx, cy, c=trends,
s=50, lw=0)
plt.colorbar()
plt.axis('equal')
plt.show()
In [20]:
gurla_pts_df = utils.tri_dict_to_df(gc, fault_name='gurla')
In [22]:
gurla_pts_df.to_csv('../data/fault_data/gurla_pts_df.csv')
In [24]:
with open('../data/fault_data/gurla_tris.geojson', 'w') as f:
json.dump(gc, f)