In [1]:
import numpy as np
import pandas as pd
import json
import halfspace.projections as hsp
from halfspace.scripts import *
  • load point dataframes
  • load geojson tris
  • get point stresses
  • average
  • make new df of average stresses, using tri center coords
  • merge w/ other fault df(s)

In [2]:
tib_dog_pts = pd.read_csv('../data/fault_data/tib_dog_pts_df.csv')
gurla_pts = pd.read_csv('../data/fault_data/gurla_pts_df.csv')

In [3]:
tdg_pts = pd.concat((tib_dog_pts, gurla_pts), axis=0, ignore_index=True)

tdg_pts.columns = ['tri', 'vertex', 'east', 'north', 'depth', 'fault_name']

In [4]:
tdg_pts.tail()


Out[4]:
tri vertex east north depth fault_name
2931 340 m -54140.576926 3406463.683422 -18788.333333 gurla
2932 341 a -57886.265887 3407087.521804 -20455.000000 gurla
2933 341 b -54899.206517 3407362.546041 -20455.000000 gurla
2934 341 c -55140.956705 3405082.406588 -15455.000000 gurla
2935 341 m -55975.476370 3406510.824811 -18788.333333 gurla

In [5]:
with open('../data/fault_data/tib_dog_tris.geojson') as f:
    tdt = json.load(f)
    
with open('../data/fault_data/gurla_tris.geojson') as f:
    gt = json.load(f)

In [6]:
stress_file = '../../nepal_2015/stress_arrays/nepal_topo_stress.h5'
stress_meta_file = '../../nepal_2015/stress_arrays/stress_calcs_meta.json'
dem_meta_file = '../../nepal_2015/data/dem/s_tibet_dem_meta.json'

In [7]:
tdg = calc_topo_stresses_at_pts(tdg_pts, stress_file=stress_file,
                                dem_meta_file=dem_meta_file,
                                stress_meta_file=stress_meta_file,
                                resolve_stress=False)


interpolating zz stresses
interpolating xy stresses
interpolating xz stresses
interpolating yz stresses
interpolating xx stresses
interpolating yy stresses

In [8]:
tdg.head()


Out[8]:
tri vertex east north depth fault_name zz_stress xy_stress xz_stress yz_stress xx_stress yy_stress
0 0 a 124496.100000 3172901.700000 4000.000000 tib_dog 0 0 0 0 0 0
1 0 b 125070.300000 3172690.800000 3743.000000 tib_dog 0 0 0 0 0 0
2 0 c 124256.800000 3172817.900000 4156.000000 tib_dog 0 0 0 0 0 0
3 0 m 124607.733333 3172803.466667 3966.333333 tib_dog 0 0 0 0 0 0
4 1 a 124058.100000 3173298.000000 4000.000000 tib_dog 0 0 0 0 0 0

In [9]:
tdg_avg = tdg.groupby(('fault_name', 'tri')).mean()[['zz_stress',
                                                     'xy_stress',
                                           'xz_stress','yz_stress',
                                           'xx_stress','yy_stress']]

In [10]:
tdg_avg.head()


Out[10]:
zz_stress xy_stress xz_stress yz_stress xx_stress yy_stress
fault_name tri
gurla 0 122.365046 1.262838 -3.130255 -7.833799 67.215418 67.963298
1 122.802498 1.562838 -2.546067 -7.169587 68.833418 69.845652
2 123.452559 2.246759 -1.578494 -5.545806 72.404882 74.020897
3 123.926082 2.731499 -1.128569 -4.596917 74.340048 76.356889
4 122.374934 1.156617 -3.483587 -7.801336 66.988821 67.718253

In [11]:
tdg_avg['i'] = range(tdg_avg.shape[0])
tdg_avg['fault_name'],tdg_avg['tri'] = zip(*tdg_avg.index.values)

In [12]:
tdg_avg.set_index('i', inplace=True)

In [13]:
tdg_avg.head(), tdg_avg.tail()


Out[13]:
(    zz_stress  xy_stress  xz_stress  yz_stress  xx_stress  yy_stress  \
 i                                                                      
 0  122.365046   1.262838  -3.130255  -7.833799  67.215418  67.963298   
 1  122.802498   1.562838  -2.546067  -7.169587  68.833418  69.845652   
 2  123.452559   2.246759  -1.578494  -5.545806  72.404882  74.020897   
 3  123.926082   2.731499  -1.128569  -4.596917  74.340048  76.356889   
 4  122.374934   1.156617  -3.483587  -7.801336  66.988821  67.718253   
 
   fault_name  tri  
 i                  
 0      gurla    0  
 1      gurla    1  
 2      gurla    2  
 3      gurla    3  
 4      gurla    4  ,
       zz_stress  xy_stress  xz_stress  yz_stress  xx_stress  yy_stress  \
 i                                                                        
 729  102.925521  -0.030489  -7.277902 -10.392833  72.365598  74.368759   
 730  101.775826  -0.468165  -2.945500 -10.891380  73.363035  77.416506   
 731  100.918067  -0.009567  -4.248177 -10.401864  73.603741  76.640471   
 732  101.976276   0.130320  -5.909874 -10.411796  72.501793  74.848722   
 733  103.119852   0.064570  -6.424231 -10.678737  71.582790  73.881284   
 
     fault_name  tri  
 i                    
 729    tib_dog  387  
 730    tib_dog  388  
 731    tib_dog  389  
 732    tib_dog  390  
 733    tib_dog  391  )

In [14]:
fd = {'gurla': gt, 'tib_dog':tdt} # fault dict

In [15]:
gt['features'][0]


Out[15]:
{'geometry': {'coordinates': [[[-55783.583006836, 3342716.71580536, -14985.0],
    [-54677.8965490291, 3339938.7432105276, -19768.0],
    [-56725.8969340146, 3341331.705430983, -19747.0],
    [-55783.583006836, 3342716.71580536, -14985.0]]],
  'type': 'Polygon'},
 'properties': {'POINTA': 130.0,
  'POINTB': 204.0,
  'POINTC': 203.0,
  'area_sq_km': 6.251594281640157,
  'center': [-55729.12549662657, 3341329.0548156234, -18166.666666666668],
  'dip': 70.61925626388545,
  'rake': -123.58798416497373,
  'rake_err': 10,
  'strike': 124.05097108466613},
 'type': 'Feature'}

In [16]:
def get_row(row):
    g = fd[row.fault_name]
    p = g['features'][row.tri]['properties']
    o = {}
    o['east'], o['north'], o['depth'] = p['center']
    o['strike'] = p['strike']
    o['dip'] = p['dip']
    o['rake'] = p['rake']
    o['area_sq_km'] = p['area_sq_km']
    
    return pd.Series(o)

In [17]:
get_row(tdg_avg.iloc[0])


Out[17]:
area_sq_km          6.251594
depth          -18166.666667
dip                70.619256
east           -55729.125497
north         3341329.054816
rake             -123.587984
strike            124.050971
dtype: float64

In [18]:
tdg_avg.apply(get_row, axis=1).head()


Out[18]:
area_sq_km depth dip east north rake strike
i
0 6.251594 -18166.666667 70.619256 -55729.125497 3341329.054816 -123.587984 124.050971
1 7.592646 -16527.333333 73.934283 -54657.964479 3341155.562365 -112.795196 130.864597
2 7.633754 -13288.333333 68.460230 -53918.655475 3341974.752790 -118.845817 131.186418
3 7.399779 -11725.000000 69.749617 -54321.288308 3343089.493691 -119.903116 128.541681
4 7.278628 -18510.000000 70.715283 -57216.150551 3342311.810013 -117.256531 130.162511

In [19]:
tri_stress = pd.merge(tdg_avg, tdg_avg.apply(get_row, axis=1),
                       left_index=True, right_index=True)

In [20]:
tri_stress.head()


Out[20]:
zz_stress xy_stress xz_stress yz_stress xx_stress yy_stress fault_name tri area_sq_km depth dip east north rake strike
i
0 122.365046 1.262838 -3.130255 -7.833799 67.215418 67.963298 gurla 0 6.251594 -18166.666667 70.619256 -55729.125497 3341329.054816 -123.587984 124.050971
1 122.802498 1.562838 -2.546067 -7.169587 68.833418 69.845652 gurla 1 7.592646 -16527.333333 73.934283 -54657.964479 3341155.562365 -112.795196 130.864597
2 123.452559 2.246759 -1.578494 -5.545806 72.404882 74.020897 gurla 2 7.633754 -13288.333333 68.460230 -53918.655475 3341974.752790 -118.845817 131.186418
3 123.926082 2.731499 -1.128569 -4.596917 74.340048 76.356889 gurla 3 7.399779 -11725.000000 69.749617 -54321.288308 3343089.493691 -119.903116 128.541681
4 122.374934 1.156617 -3.483587 -7.801336 66.988821 67.718253 gurla 4 7.278628 -18510.000000 70.715283 -57216.150551 3342311.810013 -117.256531 130.162511

In [21]:
tri_stress.rename(columns={'tri':'point_index'}, inplace=True)

In [22]:
tri_stress = resolve_stresses(tri_stress)

In [23]:
tri_stress.shape


Out[23]:
(734, 18)

In [24]:
ts = tri_stress.set_index(['fault_name', 'point_index'])

In [25]:
ts.head()


Out[25]:
zz_stress xy_stress xz_stress yz_stress xx_stress yy_stress area_sq_km depth dip east north rake strike tau_dd tau_ss sig_nn
fault_name point_index
gurla 0 122.365046 1.262838 -3.130255 -7.833799 67.215418 67.963298 6.251594 -18166.666667 70.619256 -55729.125497 3341329.054816 -123.587984 124.050971 -10.308706 -0.711939 69.627094
1 122.802498 1.562838 -2.546067 -7.169587 68.833418 69.845652 7.592646 -16527.333333 73.934283 -54657.964479 3341155.562365 -112.795196 130.864597 -7.784614 -0.500037 71.159574
2 123.452559 2.246759 -1.578494 -5.545806 72.404882 74.020897 7.633754 -13288.333333 68.460230 -53918.655475 3341974.752790 -118.845817 131.186418 -12.552396 -0.437116 78.443939
3 123.926082 2.731499 -1.128569 -4.596917 74.340048 76.356889 7.399779 -11725.000000 69.749617 -54321.288308 3343089.493691 -119.903116 128.541681 -11.568097 -0.336551 80.918145
4 122.374934 1.156617 -3.483587 -7.801336 66.988821 67.718253 7.278628 -18510.000000 70.715283 -57216.150551 3342311.810013 -117.256531 130.162511 -10.359369 -0.626613 69.307510

In [26]:
for i, feat in enumerate(tdt['features']):
    stresses = ts.loc[('tib_dog', i), :]
    feat['properties']['tau_d'] = stresses['tau_dd']
    feat['properties']['tau_s'] = stresses['tau_ss']
    feat['properties']['sigma_n'] = stresses['sig_nn']
    feat['properties']['topo_stress_rake'] = hsp.get_rake_from_shear_components(
                                                stresses['tau_dd'], stresses['tau_ss'])
    feat['properties']['tau'] = np.sqrt((stresses['tau_dd']**2 + stresses['tau_ss']**2))
    
for i, feat in enumerate(gt['features']):
    stresses = ts.loc[('gurla', i), :]
    feat['properties']['tau_d'] = stresses['tau_dd']
    feat['properties']['tau_s'] = stresses['tau_ss']
    feat['properties']['sigma_n'] = stresses['sig_nn']
    feat['properties']['topo_stress_rake'] = hsp.get_rake_from_shear_components(
                                                stresses['tau_dd'], stresses['tau_ss'])
    feat['properties']['tau'] = np.sqrt((stresses['tau_dd']**2 + stresses['tau_ss']**2))

In [27]:
tri_stress = tri_stress[tri_stress.tau_dd != 0.]

In [28]:
tri_stress.shape


Out[28]:
(610, 18)

In [29]:
tri_stress.slip_m = 0.1

In [30]:
tri_stress.to_csv('../data/fault_data/wnfs_tris.csv', index=False)

In [31]:
tdt['features'][0]


Out[31]:
{'geometry': {'coordinates': [[[124496.1, 3172901.7, 4000.0],
    [125070.3, 3172690.8, 3743.0],
    [124256.8, 3172817.9, 4156.0],
    [124496.1, 3172901.7, 4000.0]]],
  'type': 'Polygon'},
 'properties': {'POINTA': 195.0,
  'POINTB': 11.0,
  'POINTC': 12.0,
  'area_sq_km': 0.058032022780406034,
  'center': [124607.73333333334, 3172803.466666667, 3966.3333333333335],
  'dip': 31.852099159915703,
  'rake': -121.26474086089195,
  'rake_err': 15,
  'sigma_n': 0.0,
  'strike': 332.71822991430054,
  'tau': 0.0,
  'tau_d': 0.0,
  'tau_s': 0.0,
  'topo_stress_rake': 180.0},
 'type': 'Feature'}

In [32]:
with open('../data/fault_data/tib_dog_tri_stresses.geojson', 'w') as f:
    json.dump(tdt, f)
    
with open('../data/fault_data/gurla_tri_stresses.geojson', 'w') as f:
    json.dump(gt, f)

In [33]:
import matplotlib.pyplot as plt

In [34]:
tri_stress[tri_stress.fault_name=='tib_dog'].rake.hist(bins=20)
plt.show()

In [35]:
tri_stress[tri_stress.fault_name=='gurla'].tau_dd.hist(bins=20)
plt.show()

In [36]:
ts[['xz_stress', 'yz_stress']].describe()


Out[36]:
xz_stress yz_stress
count 734.000000 734.000000
mean -4.389081 -3.304068
std 5.835849 6.720918
min -26.264985 -34.373000
25% -8.448389 -6.257592
50% -1.874759 -1.698347
75% 0.000000 0.000000
max 4.991560 28.835635