Rocks and Water

Topographic and tectonic stress inversions with halfspace: Preparing a coseismic slip model

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)

Fault Slip Data Preparation

This section will cover taking a text file with a fault slip model and making a pandas DataFrame that is correctly formatted for doing the topographic and tectonic stress calculations. We will take a multi-segment model and make new files for each segment, and then make a concatenated file. We have to do it this way so that we can associate each segment with its parameters that are not already included as columns, such as strike and dip.

Some editing will need to be done in the user's text editor of choice.

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.

In [1]:
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).

In [2]:
data_dir = '../test_data/'

s1 = pd.read_csv(data_dir + 'baloch_s1.fsp', delim_whitespace=True, 
                 skipinitialspace=True)

now check the formatting:

In [3]:
s1.head()
Out[3]:
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.

In [4]:
s1['segment'] = 1
s1['strike'] = 203.0
s1['dip'] = 60.0
In [5]:
s1.head()
Out[5]:
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.

In [6]:
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)
In [7]:
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

In [8]:
bdf = pd.concat([s1, s2, s3, s4, s5, s6, s7])

bdf.head()
Out[8]:
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.

In [9]:
bdf.columns = ['lat', 'lon', 'rake', 'rise', 'slip_m', 'rupture_time',
               'x', 'y', 'z', 'dip', 'segment', 'strike']
In [10]:
bdf.head(10)
Out[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.

In [11]:
# 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()
Out[11]:
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.

In [12]:
bdf.to_csv(data_dir + 'baloch_fault_model.csv')

Table of Contents for the halfspace tutorial

Comments