# Topographic and tectonic stress inversions with halfspace: Calculating topographic stresses on a fault

In order to complete the tutorial, please first clone/download the github repo which comes with the DEM, coseismic slip model, and some auxiliary files, all with the correct pathnames.

This part is pretty easy. It's just an interpolation of the topographic stress field $$M(x,y,z)$$ onto the fault, as represented by a coseismic slip model.

In [1]:
import numpy as np
import pandas as pd
import halfspace.projections as hsp
import halfspace.sandbox as hbx
import scipy.ndimage as nd
import h5py
import json


## Set up calculations¶

### Set some paths and filenames¶

In [2]:
data_dir = '../test_data/'
dem_dir = data_dir + 'dem/'
stress_dir = '../stress_arrays/'

fault_slip_file = data_dir + 'baloch_fault_model.csv'
stress_file = stress_dir + 'baloch_topo_stress.h5'


In [3]:
d_meta = json.load(open(dem_metadata_file, 'r')) # DEM metadata dictionary
fault_df = pd.read_csv(fault_slip_file, index_col=0)  # fault slip DataFrame

stress_db = h5py.File(stress_file, 'r')     # open wrapper for stress files
stress_shape = stress_db['xx_MPa'][:,:,0].shape   # size of the stress arrays

In [4]:
x_res = int(d_meta['x_res_m'])  # get info from metadata dict
y_res = int(d_meta['y_res_m'])

z_res = 1000
z_min = int(d_meta['x_res_m'])
z_max = z_min + 26000           # make sure this matches topo stress calcs
z_len = (z_max - z_min) / z_res + 1
z_vec = np.linspace(z_min, z_max, num=z_len)
topo_res = x_res


### Calculate some more important numbers¶

In [5]:
d_meta

Out[5]:
{'north_max': 3329810.0066404096,
'upper_left_y': 3329810.0066404096,
'y_res_m': 445.9049868175547,
'lat_min': 23.993050290002728,
'lon_max': 68.27979759187075,
'east_max': 1009621.1791108495,
'x_rotation': 0.0,
'lat_max': 29.993372267246755,
'x_res_m': 445.9049868175545,
'lon_min': 62.000053658151806,
'east_min': 398285.44218398223,
'y_rotation': 0.0,
'n_cols': 1371,
'north_min': 2653818.0466249967,
'upper_left_x': 398285.44218398223,
'n_rows': 1516}
In [6]:
e_range_dem = np.linspace(d_meta['east_min'], d_meta['east_max'], # orig dem
num=d_meta['n_cols'])
n_range_dem = np.linspace(d_meta['north_min'], d_meta['north_max'],
num=d_meta['n_rows'])

e_range = hs._centered(e_range_dem, stress_shape[1])  # cut down to final size
n_range = hs._centered(n_range_dem, stress_shape[0])

x0 = e_range[0]  # first index
y0 = n_range[0]


## Interpolate stress tensor components onto fault¶

First we have to calculate the coordinates of the fault patches in 'pixel coordinates' of the DEM for interpolation

In [7]:
fault_xyz = hbx.coord_map_inverse_3d([ fault_df.east_utm.values,
fault_df.north_utm.values,
fault_df.z.values],
x_step=x_res, x_shift=x0,
y_step=y_res, y_shift=y0,
z_step=1, z_shift=fault_df.z.min())

In [8]:
# quick check to see if the conversion is set up right.

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt

plt.figure(figsize=(14,12))

plt.pcolormesh(stress_db['zz_MPa'][:,:,0], cmap=plt.cm.gist_earth)
plt.scatter(fault_xyz[0], fault_xyz[1], c=fault_xyz[2],
s=5, lw=0)
plt.colorbar(label='depth (km)')

plt.axis('equal')

plt.show()


Good, now we can do the interpolation.

In [9]:
comp_list = ['zz', 'xy', 'xz', 'yz', 'xx', 'yy']

for comp in comp_list:
print('interpolating {} stresses'.format(comp) )
fault_df['{}_stress'.format(comp)] = nd.map_coordinates(
stress_db['{}_MPa'.format(comp)][:,:,:],
[fault_xyz[1], fault_xyz[0], fault_xyz[2]],
order=1)
stress_db.close()

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


## Project stress tensor components onto fault¶

In [10]:
fault_df['tau_dd'] = 0.
fault_df['tau_ss'] = 0.
fault_df['sig_nn'] = 0.

stress_tens_xyz = {}

for i, fault_patch in enumerate(fault_df.index):
stress_tens_xyz[i] = hsp.make_xyz_stress_tensor(
sig_xx = fault_df.xx_stress.iloc[i],
sig_yy = fault_df.yy_stress.iloc[i],
sig_zz = fault_df.zz_stress.iloc[i],
sig_xy = fault_df.xy_stress.iloc[i],
sig_xz = fault_df.xz_stress.iloc[i],
sig_yz = fault_df.yz_stress.iloc[i])

fault_df['tau_dd'].iloc[i] = hsp.dip_shear_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )

fault_df['tau_ss'].iloc[i] = hsp.strike_shear_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )

fault_df['sig_nn'].iloc[i] = hsp.normal_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )

In [11]:
fault_df.head()

Out[11]:
lat lon rake rise slip_m rupture_time x y z dip ... north_utm zz_stress xy_stress xz_stress yz_stress xx_stress yy_stress tau_dd tau_ss sig_nn
0 27.1152 65.5036 -2 11.5 1.2632 11.494 17.7058 27.2426 0.2171 60 ... 3001667.271426 34.409794 -3.208107 4.297898 -4.019554 25.450635 28.030340 -5.473018 -3.743854 34.502948
1 27.0489 65.4721 7 7.0 1.2632 8.880 14.5800 19.8786 0.2171 60 ... 2994257.958404 30.553797 -1.898560 1.930903 -1.001144 23.029290 24.616245 -3.646218 -1.720015 27.994459
2 26.9826 65.4406 10 4.0 4.0000 6.387 11.4541 12.5145 0.2171 60 ... 2986849.691310 28.060829 -1.194816 1.097168 -1.762246 20.601227 22.731160 -3.566387 -1.978959 24.825576
3 26.9164 65.4090 10 4.0 5.0000 4.239 8.3283 5.1505 0.2171 60 ... 2979453.359426 26.534736 -2.712429 1.860519 -1.841293 20.933961 22.053082 -2.722381 -2.464341 26.031898
4 26.8501 65.3775 13 11.5 6.0000 3.218 5.2024 -2.2136 0.2171 60 ... 2972047.177332 27.472174 -1.913355 3.559284 -1.770538 20.819765 20.779945 -4.271303 -1.258188 26.947084

5 rows × 23 columns

## Done with calculations. Now save and we're through.¶

In [12]:
fault_df.to_csv(fault_slip_file)