RectangularPatches class

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

Classes implementing a fault made of rectangular patches. Inherits from Fault

Args:

• name : Name of the fault.

Kwargs:

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

• lon0 : Longitude of the center of the UTM zone

• lat0 : Latitude of the center of the UTM zone

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

• verbose : Speak to me (default=True)

AverageAlongStrikeOffsets(name, insars, filename, discretized=True, smooth=None)

If the profiles have the lon lat vectors as the fault, This routines averages it and write it to an output file. Weird method, I don’t know what it does…

ExtractAlongStrikeAllDepths(filename=None, discret=0.5)

Extracts the Along Strike Variations of the creep at all depths for the discretized version.

Kwargs:

• filename : save in this file

• discret : Discretize the fault

Returns:

• None

ExtractAlongStrikeVariations(depth=0.5, origin=None, filename=None, orientation=0.0)

Extract the Along Strike Variations of slip at a given depth

Args:

• depth : Depth at which we extract the along strike variations of slip.

• origin : Computes a distance from origin. Give [lon, lat].

• filename : Saves to a file.

• orientation : defines the direction of positive distances.

Returns:

• None

ExtractAlongStrikeVariationsOnDiscretizedFault(depth=0.5, filename=None, discret=0.5, interpolation='linear')

Extracts the Along Strike variations of the slip at a given depth, resampled along the discretized fault trace.

Kwargs:

• depth : Depth at which we extract the along strike variations of slip.

• discret : Discretization length

• filename : Saves to a file.

• interpolation : Interpolation method

Returns:

• None

addPatchesFromOtherFault(fault, indexes=None)

The name of this method is pretty self-explanatory.

Args:

• fault : Another fault instance with rectangular patches.

Kwargs:

• indexes : List of indices to consider

Returns:

• None

addpatch(patch, slip=[0, 0, 0])

Adds a patch to the list.

Args:

• patch : Geometry of the patch to add

Kwargs:

• slip : List of the strike, dip and tensile slip.

Returns:

• None

buildLaplacian(verbose=False, method=None, irregular=False)

Build normalized Laplacian smoothing array. This routine is not designed for unevenly paved faults. It does not account for the variations in size of the patches.

Kwargs:

• verbose : speak to me

• method : Useless argument only here for compatibility reason

• irregular : Can be used if the parametrisation is not regular along dip (N patches above or below one patch). If True, the Laplacian takes into account that there is not only one patch above (or below)

chCoordinates(p, p_ref, angle_rad)

Returns a coordinate system change.

Args:

• p : Translation vector

• p_ref : Translation reference

• angle_rad : rotation angle

Returns:

• r_p: Rotation matrix in the form of list

Note:

This is a weird method. Who wrote that?

computeAdjacencyMat(verbose=False, patchinc='alongstrike')

Computes the adjacency matrix for the fault geometry provided by ndip x nstrike. Values of 0 indicate no adjacency while values of 1 indicate patches share an edge.

Kwargs:

• verbose : Speak to me

• patchinc : For a patch N, if patch N+1 is located along-strike, patchinc should be set to ‘alongstrike’ (default). If patch N+1 is located along-dip, patchinc should be set to ‘alongdip’.

computeArea()

Computes the area of all patches. Stores that in {self.area}

Returns:

• None

computeCoulombOnPatches(friction=0.6, sign='strike')

Computes the Coulomb failure stress change on patches. Normal stress is positive away from the center of the patch.

Kwargs:

• friction: Standard coefficient of friction

• ‘strike’ or ‘dip’

Note:

Has not been tested in a loooong time…

computeEquivRectangle()

In the case where patches are not exactly rectangular, this method computes the best rectangle that fits within patches. Stores all the equivalent rectangles in equivpatch and equivpatchll.

Returns:

• None

computeImposedTractionOnEachFaultPatch(factor=0.001, mu=30000000000.0, nu=0.25)

Uses okada92 to compute the traction change on each patch from slip on the other patches.

Note:

Has not been tested in a loooong time…

Kwargs:

• factor : Conversion fator between slip units and distance units. In a regular case, distance units are km. If the slip is in mm, then factor is 10e-6.

• mu : Shear Modulus (default is 30e9 Pa).

• nu : Poisson’s ratio (default is 0.25).

computeSlipDirection(scale=1.0, factor=1.0, ellipse=False, flipstrike=False, nsigma=1.0)

Computes the segment indicating the slip direction. Direclty stores it in self.slipdirection

Kwargs:

• scale can be a real number or a string in ‘total’, ‘strikeslip’, ‘dipslip’ or ‘tensile’

• factor is a scaling factor

• ellipse: if True: design an ellipse for each slip vector

• flipstrike: if True will flip strike direction

• nsigma: nsigma for error ellipses

computeTractionOnEachFaultPatch(factor=0.001, mu=30000000000.0, nu=0.25)

Uses okada92 to compute the traction change on each patch.

Kwargs:

• factor : Conversion fator between slip units and distance units. In a regular case, distance units are km. If the slip is in mm, then factor is 10e-6.

• mu : Shear Modulus (default is 30e9 Pa).

• nu : Poisson’s ratio (default is 0.25).

Note:

This has not been tested in a while

computetotalslip()

Computes the total slip. Stores is in self.totalslip

Returns:

• None

deletepatch(patch)

Deletes a patch.

Args:

• patch : index of the patch to remove.

Returns:

• None

deletepatches(titi)

Deletes a list of patches.

Args:

• titi : List of indexes

Returns:

• None

distanceMatrix(distance='center', lim=None, nproc=1, verbose=False)

Returns a matrix of the distances between patches.

Kwargs:

• distance : has to be ‘center’ for now. No other method is implemented for now.

• lim : if not None, list of two float, the first one is the distance above which d=lim[1].

distancePatchToPatch(patch1, patch2, distance='center', lim=None)

Measures the distance between two patches.

Args:

• patch1 : geometry of the first patch.

• patch2 : geometry of the second patch.

Kwargs:

• distance : has to be ‘center’ for now. Distance between the centers of the patches.

• lim : if not None, list of two float, the first one is the distance above which d=lim[1].

distancesMatrix(distance='center', lim=None, nproc=1, verbose=False)

Returns two matrices of the distances between patches. One along the horizontal dimension, the other for the vertical.

Kwargs:

• distance : has to be ‘center’ for now. No other method is implemented for now.

• lim : if not None, list of two float, the first one is the distance above which d=lim[1].

extrap1d(interpolator)

Linear extrapolation routine. Found on StackOverflow by sastanin.

Args:

• interpolator : An instance of scipy.interpolation.interp1d

Returns:

• ufunc : An extrapolating method

extrapolate(length_added=50, tol=2.0, fracstep=5.0, extrap='ud')

Extrapolates the surface trace. This is usefull when building deep patches for interseismic loading.

Args:

• length_added : Length to add when extrapolating.

• tol : Tolerance to find the good length.

• fracstep : control each jump size.

• extrap : combination of ‘u’ (extrapolate from the beginning) and ‘d’ (extrapolate from the end). Default is ‘ud’

geometry2patch(Lon, Lat, Depth, Strike, Dip, Length, Width, initializeSlip=True)

Builds the list of patches from lists of lon, lat, depth, strike, dip, length and width

Args:

• Lon : List of longitudes

• Lat : List of Latitudes

• Depth : List of depths (km)

• Strike : List of strike angles (degree)

• Dip : List of dip angles (degree)

• Length : List of length (km)

• Width : List of width (km)

Kwargs:

• initializeSlip : Set slip values to zero

getDips()

Returns an array of dip angles for each patch (radians)

Returns:

• dip : Array of angles in radians

getEllipse(patch, ellipseCenter=None, Npoints=100, factor=1.0, nsigma=1.0)

Compute the ellipse error given Cm for a given patch

Args:

• patch : Patch

Kwargs:

• center : center of the ellipse

• Npoints : number of points on the ellipse

• factor : scaling factor

• nsigma : will design a nsigma*sigma error ellipse

Returns:

• RE : Ellipse

getFaultVector(i, discretized=True, normal=False)

Returns the vector tangential to the fault at the i-th point of the fault (discretized if True). if normal is True, returns a normal vector (+90 counter-clockwise wrt. tangente)

Args:

• i : index of the point og the fault trace

Kwargs:

• discretized: Use the discretized fault trace

• normal: Returns a fault normal vector (default is False)

Returns:

• vector: 3D vector

getPatchPositionAlongStrike(p, discretized=True)

Returns the position of a patch along strike (distance to the first point of the fault).

Args:

• p: patch index or the patch.

Returns:

• position: Float

getPatchesThatAreUnder(xf, yf, ref='utm', tolerance=0.5)

For a list of positions, returns the patches that are directly underneath. Only works if you have a vertical fault.

Args:

• xf : fault trace x coordinate

• yf : fault trace y coordinates

Kwargs:

• ref: ‘utm’ or ‘lonlat’

• tolerance: Tolerance for finding the patches (in km)

Note:

Has not been tested in a looooon time….

getStrikes()

Returns an array of strike angle for each patch (radians).

Returns:

• strike : Array of angles in radians

getcenter(p)

Get the center of one rectangular patch.

Args:

• p : Patch geometry.

Returns:

• x,y,z : Coordinates of the center

getcenters()

Get the center of the patches.

Returns:

• centers : list of centers [x,y,z]

getpatchgeometry(patch, center=False, checkindex=True)

Returns the patch geometry as needed for okada92.

Args:

• patch : index of the wanted patch or patch;

Kwargs:

• center : if true, returns the coordinates of the center of the patch. if False, returns the UL corner.

• checkindex : Checks the index of the patch

Returns:

• x, y, depth, width, length, strike, dip

When we build the fault, the patches are not exactly rectangular. Therefore, this routine will return the rectangle that matches with the two shallowest points and that has an average dip angle to match with the other corners.

horizshrink1patch(ipatch, fixedside='south', finallength=25.0)

Takes an existing patch and shrinks its size in the horizontal direction.

Args:

• ipatch : Index of the patch of concern.

Kwargs:

• fixedside : One side has to be fixed, takes the southernmost if ‘south’, takes the northernmost if ‘north’

• finallength : Length of the final patch.

importPatches(filename, origin=[45.0, 45.0])

Builds a patch geometry and the corresponding files from a relax co-seismic file type.

Args:

• filename : Input from Relax (See Barbot and Cie on the CIG website).

Kwargs:

• origin : Origin of the reference frame used by relax. [lon, lat]

Returns:

• None

lonlat2patch(lon, lat, depth, strike, dip, length, width)

Builds a patch from its longitude {lon}, latitude {lat}, depths {depth}, strike angles {strike}, dip angles {dip}, patch length {length} and patch width {width}

Args:

• lon : Longitude of the center of the patch

• lat : Latitude of the center of the patch

• depth : Depth of the center of the patch (km)

• strike : Strike of the patch (degree)

• dip : Dip of the patch (degree)

• length : Length of the patch (km)

• width : Width of the patch (km)

Returns:

• patch : a list for patch corners

• patchll : a list patch corners in lonlat

mapFault2Fault(Map, fault)

User provides a Mapping function np.array((len(self.patch), len(fault.patch))) and a fault and the slip from the argument fault is mapped into self.slip.

Args:

• Map : 2D array for mapping function

• fault : fault object

Returns:

• None

mapSlipPlane2Plane(fault, interpolation='linear', verbose=False, addlimits=True, smooth=0.1)

Maps the slip distribution from fault onto self. Mapping is built by computing the best plane between the two faults, projecting the center of patches on that plane and doing a simple resampling. The closer the two faults, the better…

Args:

• fault : Fault that has a slip distribution

Kwargs:

• interpolation : Type of interpolation method. Can be ‘linear’, ‘cubic’ or ‘quintic’

• verbose : if True, the routine says a few things.

• addlimits : Adds the upper and lower bounds of the fault in the interpolation scheme.

• smooth : Some interpolation smoothing factor

Returns:

• None

mapUnder2Above(deepfault)

This routine is very very particular. It only works with 2 vertical faults. It Builds the mapping function from one fault to another, when these are vertical. These two faults must have the same surface trace. If the deep fault has more than one raw of patches, it might go wrong and give some unexpected results.

Args:

• deepfault : Deep section of the fault.

Returns:

• None

mergePatches(p1, p2, eps=1e-06, verbose=True)

Merges 2 patches that have common corners. This modifies directly the attribute patch

Args:

• p1 : index of the patch #1.

• p2 : index of the patch #2.

Kwargs:

• eps : tolerance value for the patch corners (in km)

• verbose : Speak to me (default is True)

patchArea(p)

Computes the area of one patch.

Args:

• p : One item of self.patch

Returns:

• area : The area of the patch

patchesUtm2LonLat()

Perform the utm to lonlat conversion for all patches.

Returns:

• None

plot(figure=134, slip='total', Fault=True, Map=True, show=True, shadedtopo=None, norm=None, linewidth=1.0, plot_on_2d=True, view=None, colorbar=True, cbaxis=[0.1, 0.2, 0.1, 0.02], cborientation='horizontal', cblabel='', drawCoastlines=True, expand=0.2, figsize=(None, None))

Plot the available elements of the fault.

Args:

• figure : Number of the figure.

• slip : which slip to plot

• Fault : True to plot the fault, False otherwise

• Map : True to plot the map, False otherwise

• equiv : plot the equivalent patches

• show : True to show the plot, False otherwise

• axesscaling : True to perform axes scaling, False otherwise

• shadedtopo : True to plot shaded topography, False otherwise

• norm : Colorbar limits for slip

• linewidth : width of the lines

• plot_on_2d : True to make a map plot of the fault, False otherwise

• colorbar : True to show colorbar, False otherwise

• cbaxis : Colorbar axis position [left, bottom, width, height]

• cborientation : Colorbar orientation (‘horizontal’ or ‘vertical’)

• cblabel : Colorbar label

• drawCoastlines: True to draw coastlines, False otherwise

• expand : How much to extend the map around the fault (degrees)

• figsize : Figure size (width, height)

pointRotation3D(iPatch, iPoint, theta, p_axis1, p_axis2)

Rotate a point with an arbitrary axis (fault tip). Used in rotatePatch.

Args:

• iPatch: index of the patch to be rotated

• iPoint: index of the patch corner (point) to be rotated

• theta : angle of rotation in degrees

• p_axis1 : first point of axis (ex: one side of a fault)

• p_axis2 : second point to define the axis (ex: the other side of a fault)

Returns:

• rotated point

Reference: ‘Rotate A Point About An Arbitrary Axis (3D)’ (Paul Bourke)

read3DrectangularGrid(filename, aggregatePatchNodes=None, square=False)

This routine read the rectangular geometry.

Format:

lon	lat	E[km]	N[km]	Dep[km]	strike	dip	length	Area	ID

Args:

• filename : Name of the text file

Kwargs:

• aggregatePatchNodes : Aggregates patche nodes that are closer than a distance (float)

• square : If square == True, length = width and format becomes

lon	lat	E[km]	N[km]	Dep[km]	strike	dip	Area	ID

readPatchesFromFile(filename, Cm=None, readpatchindex=True, inputCoordinates='lonlat', donotreadslip=False, increasingy=True)

Read patches from a GMT formatted file. This means the file is a list of patches separated by ‘>’.

Args:

• filename : Name of the file.

Kwargs:

• Cm : Posterior covariances (array nslip x nslip)

• readpatchindex : Read the index of the patch and organize them

• inputCoordinates : lonlat or utm

• donotreadslip : Do not read slip values in the file

• increasingy : if you don’t want csi to set your patches corners according to increasing y, set increasingy = False

replacePatch(patch, iPatch)

Replaces one patch by the given geometry.

Args:

• patch : Patch geometry.

• iPatch : index of the patch to replace.

Returns:

• None

rotatePatch(iPatch, theta, p_axis1, p_axis2)

Rotate a patch with an arbitrary axis (fault tip) Used by class uncertainties

Args:

• iPatch: index of the patch to be rotated

• theta : angle of rotation in degrees

• p_axis1 : first point of axis (ex: one side of a fault)

• p_axis2 : second point to define the axis (ex: the other side of a fault)

Returns:

• None

setdepth(nump=None, top=None, width=None)

Set depth patch attributes

Args:

• nump : Number of fault patches at depth.

• top : Depth of the top row

• width : Width of the patches

slip2dis(data, patch, slip=None)

Computes the surface displacement at the data location using okada.

Args:

• data : data object from gps or insar.

• patch : number of the patch that slips

Kwargs:

• slip : if a number is given, that is the amount of slip along strike. If three numbers are given, that is the amount of slip along strike, along dip and opening. If None, values from slip are taken.

splitPatch(patch)

Splits a patch in 4 patches and returns 4 new patches.

Args:

• patch : item of the list of patch

Returns:

• p1, p2, p3, p4: Four patches

splitPatches()

Split all patches in 4 patches.

Returns:

• None

splitPatchesHoriz(nPatches, equiv=False, indices=None, keepSlip=False)

Splits all the patches in nPatches Horizontally. Directly modifies the patch attribute. Default behavior re-intializes slip to zero. If keepSlip, then slip values are mapped directly.

Args:

• nPatches : Number of new patches per patch (can be a list of the size of indices or an integer).

Kwargs:

• equiv : Do it on the equivalentPatches (default False)

• indices : Specify which patches to split (list of int)

• keepSlip : Map current slip values to new patches

strikedip2normal(strike, dip)

Returns a vector normal to a plane with a given strike and dip.

Args:

• strike : Strike angle in radians

• dip : Dip angle in radians

Returns:

• list a 3 components

surfacePatches2Trace()

Takes the shallowest patches of the fault and use them to build a fault trace. Direclty modifies attributes xf, yf, lonf and latf

Returns:

• None

surfacesimulation(box=None, disk=None, err=None, lonlat=None, npoints=10, slipVec=None, verbose=False)

Takes the slip vector and computes the surface displacement that corresponds on a regular grid. Returns a gps object Has not been tested in a long time…

Kwargs:

• box : Can be a list of [minlon, maxlon, minlat, maxlat, n].

• disk : list of [xcenter, ycenter, radius, n]

• err : Errors are set randomly using a uniform distribution multiplied by {err}

• lonlat : Arrays of lat and lon. [lon, lat]

• slipVec : Specify slip

translatePatch(iPatch, tr_vector)

Translate a patch Used by class uncertainties

Args:

• iPatch: index of the patch to be rotated

• tr_vector: array, translation vector in 3D

Returns:

• None

vertshrink1patch(ipatch, finalwidth, fixedside='up')

Takes an existing patch and shrinks its size in the vertical direction. This method does not account for the dip angle and just moves the row of nodes up or down.

Args:

• ipatch : Index of the patch

• finalwidth : Final width of the patch along depth.

Kwargs:

• fixedside : Which side is fixed

writePatches2File(filename, add_slip=None, scale=1.0, patch='normal', stdh5=None, decim=1)

Writes the patch corners in a file that can be used in psxyz.

Args:

• filename : Name of the file.

Kwargs:

• add_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.

• patch : Can be ‘normal’ or ‘equiv’

• stdh5 : Get standard deviation from an h5 file

• decim : How to decimate the file to get the standard dev.

writeSlipCenter2File(filename, add_slip=None, scale=1.0, neg_depth=False, addCovar=False)

Write a psxyz, surface or greenspline compatible file with the center of each patch and magnitude slip. Scale can be a real number or a string in ‘total’, ‘strikeslip’, ‘dipslip’ or ‘tensile’

Args:

• filename : Name of the output file

Kwargs:

• add_slipPut the slip as a value at the center

  Can be None, strikeslip, dipslip, total, coupling

• scale : Multiply the slip value by a factor.

• neg_depth : if True, depth is a negative nmber

• addCovar : add the covariance Cm values to the slip for each patch

Returns:

• None

writeSlipDirection2File(filename, scale=1.0, factor=1.0, neg_depth=False, ellipse=False, flipstrike=False, nsigma=1.0)

Write a psxyz compatible file to draw lines starting from the center of each patch, indicating the direction of slip. Tensile slip is not used…

Args:

• filename : Name of the output file

Kwargs:

• scale: a real number or a string in ‘total’, ‘strikeslip’, ‘dipslip’ or ‘tensile’

• factor: scaling factor

• neg_depth: use True if depth is negative

• ellipse: if True, design error ellipse for each slip vector

• flipstrike: if True, flip strike

• nsigma: if ellipse==True, design nsigma*sigma error ellipses

Returns:

• None