Topographic and tectonic stress inversions with halfspace: Preparing a coseismic slip model
Richard Styron
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.
(Or just download this notebook)
We will be using a fault slip model from the Mw 2013 Balochistan, Pakistan earthquake. The model is from Avouac et al., 2014 and hosted at the SRCMOD Finite-Source Rupture Model Database. The model is here. We will use the .fsp file format, which is in the test_data folder.
Go into test_data and open the s2013BALOCH01AVOU.fsp file in your favorite text editor. Then, copy the data for each segment, with the actual header line:
% LAT LON X==EW Y==NS Z SLIP RAKE RISE TRUP
into a separate text file, and delete the % from the header. Name each one something like baloch_s1.fsp where s1 signifies the segment number.
If you want to, you can further format the files (removing heading and trailing spaces, concatenating whitespace, etc.) which is normally the fastest way, but we will illustrate how to do this using pandas, because pandas is awesome.
First, import open ipython (or better yet the ipython notebook so you can see the formatted DataFrames) and import the modules.
import pandas as pd
import numpy as np
import pyproj as pj
Now, we want to read in the .fsp files for each fault segment (which you made with your text editor from hacking the main .fsp file).
data_dir = '../test_data/'
s1 = pd.read_csv(data_dir + 'baloch_s1.fsp', delim_whitespace=True,
skipinitialspace=True)
now check the formatting:
s1.head()
| LAT | LON | X==EW | Y==NS | Z | SLIP | RAKE | RISE | TRUP | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | 17.7058 | 27.2426 | 0.2171 | 1.2632 | -2 | 11.5 | 11.494 |
| 1 | 27.0489 | 65.4721 | 14.5800 | 19.8786 | 0.2171 | 1.2632 | 7 | 7.0 | 8.880 |
| 2 | 26.9826 | 65.4406 | 11.4541 | 12.5145 | 0.2171 | 4.0000 | 10 | 4.0 | 6.387 |
| 3 | 26.9164 | 65.4090 | 8.3283 | 5.1505 | 0.2171 | 5.0000 | 10 | 4.0 | 4.239 |
| 4 | 26.8501 | 65.3775 | 5.2024 | -2.2136 | 0.2171 | 6.0000 | 13 | 11.5 | 3.218 |
Great. We need to do this for the other segments, but first, we need to add some new columns to the DataFrame with important info.
s1['segment'] = 1
s1['strike'] = 203.0
s1['dip'] = 60.0
s1.head()
| LAT | LON | X==EW | Y==NS | Z | SLIP | RAKE | RISE | TRUP | segment | strike | dip | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | 17.7058 | 27.2426 | 0.2171 | 1.2632 | -2 | 11.5 | 11.494 | 1 | 203 | 60 |
| 1 | 27.0489 | 65.4721 | 14.5800 | 19.8786 | 0.2171 | 1.2632 | 7 | 7.0 | 8.880 | 1 | 203 | 60 |
| 2 | 26.9826 | 65.4406 | 11.4541 | 12.5145 | 0.2171 | 4.0000 | 10 | 4.0 | 6.387 | 1 | 203 | 60 |
| 3 | 26.9164 | 65.4090 | 8.3283 | 5.1505 | 0.2171 | 5.0000 | 10 | 4.0 | 4.239 | 1 | 203 | 60 |
| 4 | 26.8501 | 65.3775 | 5.2024 | -2.2136 | 0.2171 | 6.0000 | 13 | 11.5 | 3.218 | 1 | 203 | 60 |
Now we do the same for all of the fault segments.
s2 = pd.read_csv(data_dir + 'baloch_s2.fsp', delim_whitespace=True,
skipinitialspace=True)
s3 = pd.read_csv(data_dir + 'baloch_s3.fsp', delim_whitespace=True,
skipinitialspace=True)
s4 = pd.read_csv(data_dir + 'baloch_s4.fsp', delim_whitespace=True,
skipinitialspace=True)
s5 = pd.read_csv(data_dir + 'baloch_s5.fsp', delim_whitespace=True,
skipinitialspace=True)
s6 = pd.read_csv(data_dir + 'baloch_s6.fsp', delim_whitespace=True,
skipinitialspace=True)
s7 = pd.read_csv(data_dir + 'baloch_s7.fsp', delim_whitespace=True,
skipinitialspace=True)
s2['segment'] = 2
s2['strike'] = 216.0
s2['dip'] = 70.0
s3['segment'] = 3
s3['strike'] = 223.0
s3['dip'] = 47.0
s4['segment'] = 4
s4['strike'] = 230.0
s4['dip'] = 47.0
s5['segment'] = 5
s5['strike'] = 238.0
s5['dip'] = 56.0
s6['segment'] = 6
s6['strike'] = 203.0
s6['dip'] = 56.0
s7['segment'] = 7
s7['strike'] = 251.0
s7['dip'] = 56.0
Combine slip models into single dataframe¶
bdf = pd.concat([s1, s2, s3, s4, s5, s6, s7])
bdf.head()
| LAT | LON | RAKE | RISE | SLIP | TRUP | X==EW | Y==NS | Z | dip | segment | strike | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | -2 | 11.5 | 1.2632 | 11.494 | 17.7058 | 27.2426 | 0.2171 | 60 | 1 | 203 |
| 1 | 27.0489 | 65.4721 | 7 | 7.0 | 1.2632 | 8.880 | 14.5800 | 19.8786 | 0.2171 | 60 | 1 | 203 |
| 2 | 26.9826 | 65.4406 | 10 | 4.0 | 4.0000 | 6.387 | 11.4541 | 12.5145 | 0.2171 | 60 | 1 | 203 |
| 3 | 26.9164 | 65.4090 | 10 | 4.0 | 5.0000 | 4.239 | 8.3283 | 5.1505 | 0.2171 | 60 | 1 | 203 |
| 4 | 26.8501 | 65.3775 | 13 | 11.5 | 6.0000 | 3.218 | 5.2024 | -2.2136 | 0.2171 | 60 | 1 | 203 |
Rename columns so they're easier to type. Also, pay attention because the column order got shifted in transit somehow.
bdf.columns = ['lat', 'lon', 'rake', 'rise', 'slip_m', 'rupture_time',
'x', 'y', 'z', 'dip', 'segment', 'strike']
bdf.head(10)
| lat | lon | rake | rise | slip_m | rupture_time | x | y | z | dip | segment | strike | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | -2 | 11.5 | 1.2632 | 11.494 | 17.7058 | 27.2426 | 0.2171 | 60 | 1 | 203 |
| 1 | 27.0489 | 65.4721 | 7 | 7.0 | 1.2632 | 8.880 | 14.5800 | 19.8786 | 0.2171 | 60 | 1 | 203 |
| 2 | 26.9826 | 65.4406 | 10 | 4.0 | 4.0000 | 6.387 | 11.4541 | 12.5145 | 0.2171 | 60 | 1 | 203 |
| 3 | 26.9164 | 65.4090 | 10 | 4.0 | 5.0000 | 4.239 | 8.3283 | 5.1505 | 0.2171 | 60 | 1 | 203 |
| 4 | 26.8501 | 65.3775 | 13 | 11.5 | 6.0000 | 3.218 | 5.2024 | -2.2136 | 0.2171 | 60 | 1 | 203 |
| 5 | 26.7838 | 65.3459 | -8 | 11.5 | 6.0000 | 4.239 | 2.0766 | -9.5776 | 0.2171 | 60 | 1 | 203 |
| 6 | 27.1222 | 65.4851 | 31 | 5.5 | 1.2632 | 9.831 | 15.8648 | 28.0241 | 3.6812 | 60 | 1 | 203 |
| 7 | 27.0559 | 65.4535 | 25 | 8.5 | 1.1579 | 7.450 | 12.7390 | 20.6600 | 3.6812 | 60 | 1 | 203 |
| 8 | 26.9897 | 65.4220 | 16 | 5.5 | 3.4737 | 5.111 | 9.6131 | 13.2960 | 3.6812 | 60 | 1 | 203 |
| 9 | 26.9234 | 65.3904 | 19 | 4.0 | 5.2105 | 2.914 | 6.4873 | 5.9319 | 3.6812 | 60 | 1 | 203 |
Project lat, lon into UTM¶
Because all of our calculations are done in meters, we need to project the geographic coordinates into a projection. We're going to use UTM here (UTM Zone 41 N) because it's easy. It is possible to make a custom projection that might be slightly more accurate, but that might not be necessary, or worth the time. YMMV.
We will use the PyProj library to do the projections. It can use custom projections, but it can also simply use the EPSG number, which greatly simplifies things on our end.
# First, initiate (or instantiate, more accurately) the projection classes.
wgs84 = pj.Proj(init='epsg:4326')
utm41 = pj.Proj(init='epsg:32641')
# Now, transform coordinates
bdf['east_utm'], bdf['north_utm'] = pj.transform(wgs84, utm41, # old_proj, new_proj
bdf.lon.values, bdf.lat.values)
bdf.head()
| lat | lon | rake | rise | slip_m | rupture_time | x | y | z | dip | segment | strike | east_utm | north_utm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | -2 | 11.5 | 1.2632 | 11.494 | 17.7058 | 27.2426 | 0.2171 | 60 | 1 | 203 | 748188.019649 | 3001667.271426 |
| 1 | 27.0489 | 65.4721 | 7 | 7.0 | 1.2632 | 8.880 | 14.5800 | 19.8786 | 0.2171 | 60 | 1 | 203 | 745208.616285 | 2994257.958404 |
| 2 | 26.9826 | 65.4406 | 10 | 4.0 | 4.0000 | 6.387 | 11.4541 | 12.5145 | 0.2171 | 60 | 1 | 203 | 742225.243808 | 2986849.691310 |
| 3 | 26.9164 | 65.4090 | 10 | 4.0 | 5.0000 | 4.239 | 8.3283 | 5.1505 | 0.2171 | 60 | 1 | 203 | 739227.769404 | 2979453.359426 |
| 4 | 26.8501 | 65.3775 | 13 | 11.5 | 6.0000 | 3.218 | 5.2024 | -2.2136 | 0.2171 | 60 | 1 | 203 | 736236.493174 | 2972047.177332 |
Done making fault model. Let's save to CSV.¶
bdf.to_csv(data_dir + 'baloch_fault_model.csv')
Table of Contents for the halfspace tutorial
- Installation and configuration
- Preparing a Digital Elevation Model file
- Preparing the coseismic slip model
- Calculating the topographic stress field
- Calculating the topographic stresses on the fault
- Bayesian inversion of tectonic stress pt. 1: Rake inversion
- Bayesian inversion of tectonic stress pt. 2: Conditions at failure