Fault class

class csi.Fault(name, utmzone=None, ellps='WGS84', lon0=None, lat0=None, verbose=True)

Parent class implementing what is common in all fault objects.

You can specify either an official utm zone number or provide longitude and latitude for a custom zone.

Args:

• name : Name of the fault.

• utmzone : UTM zone (optional, default=None)

• lon0 : Longitude defining the center of the custom utm zone

• lat0 : Latitude defining the center of the custom utm zone

• ellps : ellipsoid (optional, default=’WGS84’)

addfaults(filename)

Add some other faults to plot with the modeled one.

Args:

• filename : Name of the file. File is ascii format. First column is longitude. Second column is latitude. Separator between faults is > as in GMT style

Return:

• None

assembleCd(datas, add_prediction=None, verbose=False)

Assembles the data covariance matrices that have been built for each data structure.

Args:

• datas : List of data instances or one data instance

Kwargs:

• add_prediction: Precentage of displacement to add to the Cd diagonal to simulate a Cp (dirty version of a prediction error covariance, see Duputel et al 2013, GJI).

• verbose : Talk to me (overwrites self.verbose)

Returns:

• None

assembleGFs(datas, polys=None, slipdir='sdt', verbose=True, custom=False, computeNormFact=True)

Assemble the Green’s functions corresponding to the data in datas. This method allows to specify which transformation is going to be estimated in the data sets, through the polys argument.

Assembled Green’s function matrix is stored in self.Gassembled

Args:

• datas : list of data sets. If only one data set is used, can be a data instance only.

Kwargs:

• polys : None, nothing additional is estimated

  • For InSAR, Optical, GPS:

    • 1: estimate a constant offset

    • 3: estimate z = ax + by + c

    • 4: estimate z = axy + bx + cy + d

  • For GPS only:

    • ‘full’ : Estimates a rotation, translation and scaling (Helmert transform).

    • ‘strain’ : Estimates the full strain tensor (Rotation + Translation + Internal strain)

    • ‘strainnorotation’ : Estimates the strain tensor and a translation

    • ‘strainonly’ : Estimates the strain tensor

    • ‘strainnotranslation’ : Estimates the strain tensor and a rotation

    • ‘translation’ : Estimates the translation

    • ‘translationrotation : Estimates the translation and a rotation

• slipdir : Directions of slip to include. Can be any combination of s (strike slip), d (dip slip), t (tensile), c (coupling)

• custom : If True, gets the additional Green’s function from the dictionary self.G[data.name][‘custom’]

• computeNormFact : bool. if True, compute new OrbNormalizingFactor. if False, uses parameters in self.OrbNormalizingFactor

• verbose : Talk to me (overwrites self.verbose)

Returns:

• None

assembled(datas, verbose=True)

Assembles a data vector for inversion using the list datas Assembled vector is stored in self.dassembled

Args:

• datas : list of data objects

Returns:

• None

buildCm(sigma, lam, lam0=None, extra_params=None, lim=None, verbose=True)

Builds a model covariance matrix using the equation described in Radiguet et al 2010. We use

\(C_m(i,j) = \frac{\sigma \lambda_0}{ \lambda }^2 e^{-\frac{||i,j||_2}{ \lambda }}\)

extra_params allows to add some diagonal terms and expand the size of the matrix, in case the fault object is also hosting the estimation of transformation parameters.

Model covariance is stored in self.Cm

Args:

• sigma : Amplitude of the correlation.

• lam : Characteristic length scale.

Kwargs:

• lam0 : Normalizing distance. if None, lam0=min(distance between patches)

• extra_params : A list of extra parameters.

• lim : Limit distance parameter (see self.distancePatchToPatch)

• verbose : Talk to me (overwrites self.verrbose)

Returns:

• None

buildCmGaussian(sigma, extra_params=None)

Builds a diagonal Cm with sigma values on the diagonal. Sigma is a list of numbers, as long as you have components of slip (1, 2 or 3). extra_params allows to add some diagonal terms and expand the size of the matrix, in case the fault object is also hosting the estimation of transformation parameters.

Model covariance is hold in self.Cm

Args:

• sigma : List of numbers the size of the slip components requried for the modeling

Kwargs:

• extra_params : a list of extra parameters.

Returns:

• None

buildCmLaplacian(lam, diagFact=None, extra_params=None, sensitivity=True, method='distance', sensitivityNormalizing=False, irregular=False)

Implements the Laplacian smoothing with sensitivity (optional) into a model covariance matrix. Description can be found in F. Ortega-Culaciati’s PhD thesis.

extra_params allows to add some diagonal terms and expand the size of the matrix, in case the fault object is also hosting the estimation of transformation parameters.

Model covariance is hold in self.Cm

Args:

• lam : Damping factor (list of size of slipdirections)

Kwargs:

• extra_params : a list of extra parameters.

• sensitivity : Weights the Laplacian by Sensitivity (default True)

• sensitivityNormalizing : Normalizing the Sensitivity?

• method : which method to use to build the Laplacian operator

• irregular : Only used for rectangular patches. Allows to account for irregular meshing along dip.

Returns:

• None

buildCmSensitivity(sigma, lam, lam0=None, extra_params=None, lim=None, verbose=True)

Builds a model covariance matrix using the equation described in Radiguet et al 2010. We use

\(C_m(i,j) = \frac{\sigma\lambda_0}{\lambda}^2 e^{-\frac{||i,j||_2}{\lambda}}\)

Then correlation length is weighted by the sensitivity matrix described in Ortega’s PhD thesis: \(S = diag(G'G)\)

Here, Sigma and Lambda are lists specifying values for the slip directions

extra_params allows to add some diagonal terms and expand the size of the matrix, in case the fault object is also hosting the estimation of transformation parameters.

Model covariance is stored in self.Cm

Args:

• sigma : Amplitude of the correlation.

• lam : Characteristic length scale.

Kwargs:

• lam0 : Normalizing distance. if None, lam0=min(distance between patches)

• extra_params : a list of extra parameters.

• lim : Limit distance parameter (see self.distancePatchToPatch)

• verbose : Talk to me (overwrites self.verrbose)

Returns:

• None

buildCmSlipDirs(sigma, lam, lam0=None, extra_params=None, lim=None, verbose=True)

Builds a model covariance matrix using the equation described in Radiguet et al 2010. Here, Sigma and Lambda are lists specifying values for the slip directions. We use

\(C_m(i,j) = \frac{\sigma\lambda_0}{\lambda}^2 e^{-\frac{||i,j||_2}{\lambda}}\)

Args:

• sigma : Amplitude of the correlation.

• lam : Characteristic length scale.

Kwargs:

• lam0 : Normalizing distance. If None, lam0=min(distance between patches)

• extra_params : A list of extra parameters.

• lim : Limit distance parameter (see self.distancePatchToPatch)

• verbose : Talk to me (overwrites self.verrbose)

Returns:

• None

buildCmXY(sigma, lam, lam0=None, extra_params=None, lim=None, verbose=True)

Builds a model covariance matrix using the equation described in Radiguet et al 2010 with a different characteristic lengthscale along the horizontal and vertical directions. We use

\(C_m(i,j) = \frac{\sigma \lambda_0}{ \lambda_x }^2 e^{-\frac{||i,j||_{x2}}{ \lambda_x }} \frac{\sigma \lambda_0}{ \lambda_z }^2 e^{-\frac{||i,j||_{z2}}{ \lambda_z }}\)

extra_params allows to add some diagonal terms and expand the size of the matrix, in case the fault object is also hosting the estimation of transformation parameters.

Model covariance is stored in self.Cm

Args:

• sigma : Amplitude of the correlation.

• lam : Characteristic length scale (lamx, lamz)

Kwargs:

• lam0 : Normalizing distance. if None, lam0=min(distance between patches)

• extra_params : A list of extra parameters.

• lim : Limit distance parameter (see self.distancePatchToPatch)

• verbose : Talk to me (overwrites self.verrbose)

Returns:

• None

buildGFs(data, vertical=True, slipdir='sd', method='homogeneous', verbose=True, convergence=None)

Builds the Green’s function matrix based on the discretized fault.

The Green’s function matrix is stored in a dictionary. Each entry of the dictionary is named after the corresponding dataset. Each of these entry is a dictionary that contains ‘strikeslip’, ‘dipslip’, ‘tensile’ and/or ‘coupling’

Args:

• data : Data object (gps, insar, optical, …)

Kwargs:

• vertical : If True, will produce green’s functions for the vertical displacements in a gps object.

• slipdir : Direction of slip along the patches. Can be any combination of s (strikeslip), d (dipslip), t (tensile) and c (coupling)

• method : Can be ‘okada’ (Okada, 1982) (rectangular patches only), ‘meade’ (Meade 2007) (triangular patches only), ‘edks’ (Zhao & Rivera, 2002), ‘homogeneous’ (Okada for rectangles, Meade for triangles)

• verbose : Writes stuff to the screen (overwrites self.verbose)

• convergence : If coupling case, needs convergence azimuth and rate [azimuth in deg, rate]

Returns:

• None

TODO: Implement the homogeneous case for the Node-based triangular GFs

calcGFsCp(datasets, edks=False, **edks_params)

Calculate Green’s Functions using Okada or EDKS Used in class uncertainties

Args:

• datasetsList of data objects

  ex: dataset=[gps]+[insar1,insar2]

Kwargs:

• edksIf True, GFs calculated using a layered Earth model calculated with EDKS.

  If False, GFs with Okada

If edks is True, please specify in **edks_params:

ex: Cp_dip(fault,datasets,[40,50],multi_segments=2,edks=True,edksdir=’PATH’,modelname=’CIA’,sourceSpacing=0.5) * modelname : xxx.edks = Filename of the EDKS kernels * sourceSpacing : source spacing to calculate the Green’s Functions

OR sourceNumber : Number of sources per patches. OR sourceArea : Maximum Area of the sources.

Returns:

• Gassembled

cumdis2xy(distance, recompute=True, mode='lonlat', discretized=False)

For a given {distance}, returns the x and y position along strike.

Args:

• distance : Along strike distance

Kwargs:

• recompute : recompute the interpolator

• mode : ‘lonlat’ returns lon/lat while ‘xy’ returns x and y

• discretized : use the discretized fault trace

Returns:

• xy : tuple of floats

cumdistance(discretized=False, recompute=True)

Computes the distance between the first point of the fault and every other point. The distance is cumulative along the fault.

Args:

• discretized : if True, use the discretized fault trace (default False)

• recompute : if False, just returns the attribute cumdis

Returns:

• dis : Cumulative distance array

discretize(every=2.0, tol=0.01, fracstep=0.2, xaxis='x', cum_error=True)

Refine the surface fault trace by setting a constant distance between each point. Pay attention, the fault cannot be exactly a straight line north-south. Descretized fault trace is stored in self.xi and self.yi

Kwargs:

• every : Spacing between each point (in km)

• tol : Tolerance in the spacing (in km)

• fracstep : fractional step in the chosen direction for the discretization optimization

• xaxis : Axis for the discretization. ‘x’= use x as the axis, ‘y’= use y as the axis

• cum_error : if True, accounts for cumulated error to define the axis bound for the last patch

Returns:

• None

distance2trace(lon, lat, discretized=False, coord='ll', recompute=True)

Computes the distance between a point and the trace of a fault. This is a slow method, so it has been recoded in a few places throughout the whole library.

Args:

• lon : Longitude of the point.

• lat : Latitude of the point.

Kwargs:

• discretized : Uses the discretized trace.

• recompute : recompute the cumulative distance

• coord : if ‘ll’ or ‘lonlat’, input in degree. If ‘xy’ or ‘utm’, input in km

Returns:

• dalong : Distance to the first point of the fault along the fault

• dacross : Shortest distance between the point and the fault

dropPointSources()

Drops point sources along the fault. Point sources can then be used to compute GFs using the EDKS software.

The process is controlled by the attributes:

• self.sourceSpacing : Distance between sources

• self.sourceArea : Area of the sources

• self.sourceNumber : Number of sources per patch

One needs to set at least one of those three attributes.

Sources are saved in self.plotSources and self.edksSources

Returns:

• None

duplicateFault()

Returns a full copy (copy.deepcopy) of the fault object.

Return:

• fault : fault object

edksGFs(data, vertical=True, slipdir='sd', verbose=True, convergence=None, method='fortran')

Builds the Green’s functions based on the solution by Zhao & Rivera 2002. The corresponding functions are in the EDKS code that needs to be installed and the executables should be found in the directory set by the environment variable EDKS_BIN.

A few variables need to be set in before running this method

Required:

• self.kernelsEDKS : Filename of the EDKS kernels.

One of the Three:

• self.sourceSpacing : Spacing between the sources in each patch.

• self.sourceNumber : Number of sources per patches.

• self.sourceArea : Maximum Area of the sources.

Args:

• data : Data object

Kwargs:

• vertical : If True, will produce green’s functions for the vertical displacements in a gps object.

• slipdir : Direction of slip along the patches. Can be any combination of s (strikeslip), d (dipslip), t (tensile) and c (coupling)

• verbose : Writes stuff to the screen (overwrites self.verbose)

• convergence : If coupling case, needs convergence azimuth and rate [azimuth in deg, rate]

Returns:

• G : Dictionary of the built Green’s functions

emptyGFs(data, vertical=True, slipdir='sd', verbose=True)

Build zero GFs.

Args:

• data : Data object (gps, insar, optical, …)

Kwargs:

• vertical : If True, will produce green’s functions for the vertical displacements in a gps object.

• slipdir : Direction of slip along the patches. Can be any combination of s (strikeslip), d (dipslip), t (tensile) and c (coupling)

• verbose : Writes stuff to the screen (overwrites self.verbose)

Returns:

• G : Dictionnary of GFs

estimateSeismicityRate(earthquake, extra_div=1.0, epsilon=1e-05)

Counts the number of earthquakes per patches and divides by the area of the patches. Sets the results in

self.earthquakeInPatch (Number of earthquakes per patch) and self.seismicityRate (Seismicity rate for this patch)

Args:

• earthquake : seismiclocation object

Kwargs:

• extra_div : Extra divider to get the seismicity rate.

• epsilon : Epsilon value for precision of earthquake location.

Returns:

• None

file2trace(filename, utm=False, header=0)

Reads the fault trace from a text file (ascii 2 columns)

• If utm is False, format is Lon Lat

• If utm is True, format is X Y (in km)

Args:

• filename : Name of the fault file.

Kwargs:

• utm : Specify nature of coordinates

• header : Number of lines to skip at the beginning of the file

Returns:

• None

gaussianSlipSmoothing(length)

Smoothes the slip distribution using a Gaussian filter. Smooth slip distribution is in self.slip

Args:

• length : Correlation length.

Returns:

• None

getindex(p)

Returns the index of a patch.

Args:

• p : Patch from a fault object.

Returns:

• iout : index of the patch

getslip(p)

Returns the slip vector for a patch or tent

Args:

• p : patch or tent

Returns:

• iout : Index of the patch or tent

homogeneousGFs(data, vertical=True, slipdir='sd', verbose=True, convergence=None)

Builds the Green’s functions for a homogeneous half-space.

If your patches are rectangular, Okada’s formulation is used (Okada, 1982) If your patches are triangular, Meade’s formulation is used (Meade, 2007)

Args:

• data : Data object (gps, insar, optical, …)

Kwargs:

• vertical : If True, will produce green’s functions for the vertical displacements in a gps object.

• slipdir : Direction of slip along the patches. Can be any combination of s (strikeslip), d (dipslip), t (tensile) and c (coupling)

• verbose : Writes stuff to the screen (overwrites self.verbose)

• convergence : If coupling case, needs convergence azimuth and rate [azimuth in deg, rate]

Returns:

• G : Dictionary of the built Green’s functions

initializeEmptyFault()

Initializes what is required for a fualt with no patches

Returns:

• None

initializeslip(n=None, values=None)

Re-initializes the fault slip array to zero values. Slip array will be the size of the number of patches/tents times the 3 components of slip (strike-slip, dip slip and tensile).

• 1st Column is strike slip

• 2nd Column is dip slip

• 3rd Column is tensile

Kwargs:

• n : Number of slip values. If None, it’ll take the number of patches.

• values : Can be ‘depth’, ‘strike’, ‘dip’, ‘length’, ‘width’, ‘area’, ‘index’ or a numpy array. The array can be of size (n,3) or (n,1)

Returns:

• None

interpolateSlip(lon, lat, slip='strikeslip', rebuild=False, coord='LL')

Creates an interpolator and interpolates the slip values to the position given in {lon} and {lat}.

Args:

• lon : Longitude

• lat : Latitude

Kwargs:

• slip : Which slip value to take

• rebuild : Rebuild the interpolator

• coord : LL or utm (in km)

patch2ll()

Takes all the patches in self.patch and convert them to lonlat. Patches are stored in self.patchll

Returns:

• None

patch2xy()

Takes all the patches in self.patchll and convert them to xy Patches are stored in self.patch

Returns:

• None

pickle(filename)

Pickle myself.

Args:

• filename : Name of the pickle fault.

readPointSourcesFromPickle(filename)

Reads the point sources for computing Green’s functions with EDKS from a pickle file. Sets the sources in self.edksSources

Args:

• filename : Name of the pickle file

Returns:

• None

rotateGFs(data, azimuth)

For the data set data, returns the rotated GFs so that dip slip motion is aligned with the azimuth. It uses the Greens functions stored in self.G[data.name].

Args:

• data : Name of the data set.

• azimuth : Direction in which to rotate the GFs

Returns:

• rotatedGar : GFs along the azimuth direction

• rotatedGrp : GFs in the direction perpendicular to the azimuth direction

saveData(dtype='d', outputDir='.')

Saves the Data in binary files.

Kwargs:

• dtype : Format of the binary data saved. ‘d’ for double. ‘f’ for np.float32

• outputDir : Directory to save binary data

Returns:

• None

saveGFs(dtype='d', outputDir='.', suffix={'coupling': 'Coupling', 'dipslip': 'DS', 'strikeslip': 'SS', 'tensile': 'TS'})

Saves the Green’s functions in different files.

Kwargs:

• dtype : Format of the binary data saved. ‘d’ for double. ‘f’ for np.float32

• outputDir : Directory to save binary data.

• suffix : suffix for GFs name (dictionary)

Returns:

• None

setCustomGFs(data, G)

Sets a custom Green’s Functions matrix in the G dictionary.

Args:

• data : Data concerned by the Green’s function

• G : Green’s function matrix

Returns:

• None

setGFs(data, strikeslip=[None, None, None], dipslip=[None, None, None], tensile=[None, None, None], coupling=[None, None, None], vertical=False, synthetic=False)

Stores the input Green’s functions matrices into the fault structure.

These GFs are organized in a dictionary structure in self.G Entries of self.G are the data set names (data.name). Entries of self.G[data.name] are ‘strikeslip’, ‘dipslip’, ‘tensile’ and/or ‘coupling’

If you provide GPS GFs, those are organised with E, N and U in lines

If you provide Optical GFs, those are organised with E and N in lines

If you provide InSAR GFs, these need to be projected onto the LOS direction already.

Args:

• data : Data structure

Kwargs:

• strikeslip : List of matrices of the Strikeslip Green’s functions

• dipslip : List of matrices of the dipslip Green’s functions

• tensile : List of matrices of the tensile Green’s functions

• coupling : List of matrices of the coupling Green’s function

Returns:

• None

setGFsFromFile(data, strikeslip=None, dipslip=None, tensile=None, coupling=None, custom=None, vertical=False, dtype='d', inDir='.')

Sets the Green’s functions reading binary files. Be carefull, these have to be in the good format (i.e. if it is GPS, then GF are E, then N, then U, optional, and if insar, GF are projected already). Basically, it will work better if you have computed the GFs using csi…

Args:

• data : Data object

kwargs:

• strikeslip : File containing the Green’s functions for strikeslip related displacements.

• dipslip : File containing the Green’s functions for dipslip related displacements.

• tensile : File containing the Green’s functions for tensile related displacements.

• coupling : File containing the Green’s functions for coupling related displacements.

• vertical : Deal with the UP component (gps: default is false, insar: it will be true anyway).

• dtype : Type of binary data. ‘d’ for double/float64. ‘f’ for np.float32

Returns:

• None

setTrace(delta_depth=0.0, sort='y')

Uses the patches to build a fault trace. Fault trace is made of the vertices that are shallower than fault top + delta_depth Fault trace is in self.xf and self.yf

Args:

• delta_depth : Depth extension below top of the fault

setmu(model_file, tents=False)

Gets the shear modulus corresponding to each patch using a model file from the EDKS software. Shear moduli are set in self.mu

The model file format is as follows:

N	F
RHO_1	VP_1	VS_1	TH_1
RHO_2	VP_2	VS_2	TH_2
…	…	…	…
RHO_N	VP_N	VS_N	TH_N

where N is the number of layers, F a conversion factor to SI units RHO_i is the density of the i-th layer VP_i is the P-wave velocity in the i-th layer VS_i is the S-wave velocity in the i-th layer TH_i is the thickness of the i-th layer

Args:

• model_file : path to model file

• tents : if True, set mu values every point source in patches

Returns:

• None

slipIntegrate(slip=None)

Integrates slip on the patch by simply multiplying slip by the patch area. Sets the results in self.volume

Kwargs:

• slip : Can be strikeslip, dipslip, tensile, coupling or a list/array of floats.

Returns:

• None

smoothTrace(winsize=10, std=1.0, discretized=True)

Smoothes the trace of the fault using a median filter (scipy.signal.medfilt)

Kwargs:

• std : Standard deviation of the Gaussian kernel in km

• discretized : Use the discretized fault trace (self.xi and self.yi) (default True)

Returns:

• None

strikeOfTrace(discretized=True, npoints=4)

Computes the strike of the fault trace from the discretized (default) fault trace.

Kwargs:

• discretized : Use the discretized fault trace (self.xi and self.yi) (default True)

• npoints : Number of points to average strike

Return:

• None. Stores the strike in self.strike

sumPatches(iPatches, finalPatch)

Takes a list of indexes of patches, sums the corresponding GFs and replace the corresponding patches by the finalPatch in self.patch

Args:

• patches : List of the patche indexes to sum

• finalPatch : Geometry of the final patch.

Returns:

• None

surfaceGFs(data, slipdir='sd', verbose=True)

Build the GFs for the surface slip case. We assume the data are within the bounds of the fault.

Args:

• data : surfaceslip data ojbect

Kwargs:

• slipdir : any combinatino of s and d. default: ‘sd’

• verbose : Default True

trace(x, y, utm=False)

Set the surface fault trace from Lat/Lon or UTM coordinates Surface fault trace is stored in self.xf, self.yf (UTM) and self.lon, self.lat (Lon/lat)

Args:

• Lon : Array/List containing the Lon points.

• Lat : Array/List containing the Lat points.

Kwargs:

• utm : If False, considers x and y are lon/lat. If True, considers x and y are utm in km

Returns:

• None

trace2ll()

Transpose the fault trace UTM coordinates into lat/lon. Lon/Lat coordinates are stored in self.lon and self.lat in degrees

Returns:

• None

trace2xy()

Transpose the fault trace lat/lon into the UTM reference. UTM coordinates are stored in self.xf and self.yf in km

Returns:

• None

writePatchesCenters2File(filename, slip=None, scale=1.0)

Write the patch center coordinates in an ascii file the file format is so that it can by used directly in psxyz (GMT).

Args:

• filename : Name of the file.

Kwargs:

• slip : Put the slip as a value for the color. Can be None, strikeslip, dipslip, total, coupling

• scale : Multiply the slip value by a factor.

Retunrs:

• None

writePointSources2Pickle(filename)

Writes the point sources to a pickle file. Always writes the Facet based point sources.

Args:

• filename : Name of the pickle file.

Returns:

• None

writeTrace2File(filename, ref='lonlat')

Writes the trace to a file. Format is ascii with two columns with either lon/lat (in degrees) or x/y (utm in km).

Args:

• filename : Name of the file

Kwargs:

• ref : can be lonlat or utm.

Returns:

• None