startTime=AbsoluteTime[];

Seismo-Lineament Geomorphology

This code is copyright © 2010 by Vincent S. Cronin.  It may be used for educational or scientific purposes with appropriate attribution.

Cronin, V.S., Millard, M.A., Seidman, L.E., and Bayliss, B.G., 2008, The Seismo-Lineament Analysis Method [SLAM] -- A reconnaissance tool to help find seismogenic faults: Environmental and Engineering Geoscience, v. 14, no. 3, p. 199-219.

Begun May 28, 2010

Introduction

The purpose of this code is to produce two hillshade images of a seismo-lineament, each illuminated perpendicular to the trend of the lineament.

The output from this application is a pair of JPEG images, whose names are defined by the user in a gray box near the end of this code.

Explanation of changes to this code that are necessary to evaluate another fault plane solution

Places in this notebook where the user must imput data are marked with a gray box.
(1) In the "Description of input data" section, the user must supply values for the following variables:  focalLat, focalLong, focalDepthKm, reportedVertError, reportedHorizError, firstFltDipTrend, firstFltDipPlunge, secondFltDipTrend, secondFltDipPlunge, zoneMeridian.
(2) In the "Import the DEM data from a DAT file" subsection of the "Digital elevation model (DEM)" section, the user may need to modify the path information so this notebook can find the input DEM/DAT file.
(3) In the "Export data and image files" section at the end of this code, the user must supply a different name for the output JPEG graphics files to avoid over-writing the results of previous runs.

Description of input data

Data source for focal mechanism solution used in this run

Brumbaugh, David S. 2008. Seismicity and active faulting of the Kanad-Fredoni are of the southern Colorado Plateau. Journal of Geophysical Research, v. 113, B05309, doi:10.1029/2007JB005278   M5.75 event 19590721 UTC 1239:29

From focal mechanism solution

Epicenter location:  latitude (focalLat) and longitude (focalLong) in decimal degrees.

focalLat=36.800;

focalLong=-112.370;

Focal depth (focalDepthKm) in kilometers

focalDepthKm=11;

Reported vertical (reportedVertError) and horizontal (horizError) uncertainties in the location of the earthquake focus, in kilometers.

reportedVertError=3;

reportedHorizError=2;

We use the trend and plunge of the nodal plane's dip vector to define the orientation of the nodal plane in this notebook.  This avoids the potential ambiguity in using the traditional "strike, dip and dip direction" used by geologists, who may report data in a quadrant, azimuth, right-hand rule, left-hand rule, or some other format.

For example, the focal mechanism solution data were given by Brumbaugh (2008, Table 2) in terms of strike, dip and dip direction.  The orientation of nodal plane 1 was listed as 347/47NE, which has a dip vector plunging 47° and trending 347°-270°=77°.  The orientation of nodal plane 2 was listed as 353/43SW, which has a dip vector plunging 43° and trending 353°-90°=263°.

faultDipTrend=263;

faultDipPlunge=43;

The central meridian of a UTM zone is the longitude of the middle of the UTM zone.  Some common UTM central meridians are listed in the table below.

({{({{Zone}, {10}, {11}})}, {12}})     ({{({{Longitude Range}, {-126° to -120°}, {-120° to -114°}})}, {-114° to -108°}})     ({{({{Central Meridian}, {-123°}, {-117°}})}, {-111°}})

In the appropriate space below, input the numerical value of the central meridian, without adding the degree symbol.  If the study area spans two zones, use the more westerly zone's central meridian.

zoneMeridian=-111;

Some user-defined functions

makeVector[plunge_,trend_]:={Cos[plunge Degree] Sin[trend Degree],Cos[plunge Degree] Cos[trend Degree],-Sin[plunge Degree]};

vectorNorm[x_]:=(Sqrt[x.x]);

unitVector[x_]:={x[[1]]/vectorNorm[x],x[[2]]/vectorNorm[x],x[[3]]/vectorNorm[x]};

vectorAngle[a_,b_]:=ArcCos[a.b/(vectorNorm[a] vectorNorm[b])];

The lat/long to UTM conversion is after Snyder (1982), and assumes use of the North American Datum of 1927 (NAD27).  Another excellent resource for converting from lat/long to UTM and back, using a variety of datums (e.g., NAD27, WGS84 and so on) is provided by Steven Dutch at http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm, and by an online calculator at http://www.rcn.montana.edu/resources/tools/coordinates.aspx

convertToUTM[inLat_,inLong_,centMerid_]:=Module[{c1,c2,c3,c4,c5,v1,v2,v3,v4,v5,v6,utmX,utmY},c1=6378206.4;c2=0.00676866;c3=0;c4=centMerid;c5=0.9996;v1=c2/(1-c2);v2=c1/Sqrt[1-(c2*(((Sin[inLat Degree])^2)))];v3=((Tan[inLat Degree])^2);v4=v1*(((Cos[inLat Degree])^2));v5=(Cos[inLat Degree])*((inLong-c4) ((π/180)));v6=(111132.0894*inLat)-(16216.94*Sin[2*(inLat Degree)])+(17.21*Sin[4*(inLat Degree)])-(0.02*Sin[6*(inLat Degree)]);utmX=(c5*v2*(v5+(((1-v3+v4)*(v5^3))/6)+((5-(18*v3)+((v3^2))+(72*v4)-(58*v1))*(v5^5))/120))+500000;utmY=c5*(v6 - 0+(v2*Tan[inLat Degree]*(((v5^2)/2)+(((5-v3+(9*v4)+(4*((v4^2))))*(v5^4))/24)+(((61-(58*v3)+((v3^2))+(600*v4)-(330*v1))*(v5^6))/720))));{utmX,utmY}];

The following module differentiates between points that are within "width" meters from the fault plane and those that are further away.  Given unit vector N that is normal to the fault plane that passes through the origin of the coordinate system, the distance from an arbitrary point (whose position vector is P) to that plane is given by |N·P|.

pointEvaluator[xCoord_,yCoord_,zCoord_,width_,fltUNrml_,nulData_]:=Module[{locVect,distToFlt,result},locVect={xCoord,yCoord,zCoord};distToFlt=If[(zCoord<(-1000)),width,Abs[Dot[fltUNrml,locVect]]];result=If[(distToFlt≤width),zCoord,nulData];result];

Change the units of some input data

The UTM coordinates that specify the horizontal position of the earthquake focus (i.e., the epicenter location) are given in meters.  The variable focalDepth gives the depth of the earthquake focus, in meters.

focalDepth=focalDepthKm*(-1000);

The reported vertical and horizontal uncertainties in the location of the earthquake focus are given in kilometers.  The variables horizError and vertError convert these to meters.

horizError=reportedHorizError*1000;

vertError = reportedVertError*1000;

Find the unit vector normal to the fault surface at the earthquake focus

Convert the fault dip vector from trend and plunge to a unit location vector

faultDipVect=unitVector[makeVector[faultDipPlunge,faultDipTrend]];

Find the strike vector defined using the right-hand rule

faultStrike=faultDipTrend-90;

faultStrikeVect={Sin[faultStrike Degree],Cos[faultStrike Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

faultNormalVect=unitVector[faultDipVect\[Cross]faultStrikeVect];

Digital elevation model (DEM)

The digital elevation model is presented as a set of z data in meters, with x and y coordinates explicit from the position of the z datum in the data file.  The data file is derived from the U.S. Geological Survey's DEM data for the Pt. Dume, Malibu Beach and Topanga 7.5 minute quadrangles (30 m resolution).
DEM files can be acquired from the National Elevation Dataset, produced by the U.S. Geological Survey (USGS) and obtained from USGS web servers (http://seamless.usgs.gov), the Geocommunity web server (http://data.geocomm.com) or from www.gisdatadepot.com.  The file can be in ASCII (DEM) format or Spatial Data Transfer Standard (SDTS) format.

Import the DEM data from a DAT file

The file pathway to the DEM file must be specified.

mydata=Import["/Users/vince_cronin/Desktop/Lancaster/thin1991.dat"];

Read and interpret the DEM header information

Read the header information from the first six records of the DEM file.

headerData=Table[mydata[[i,j]],{i,6},{j,2}];

"ncols" is the number of columns;  label in position [[1,1]] and value in [[1,2]]

ncols=headerData[[1,2]];

"nrows" is the number of rows;  label in position [[2,1]] and value in [[2,2]]

nrows=headerData[[2,2]];

"xllcorner" is the UTM zone 11 "x" coordinate, in meters, of the lower left corner of the DEM;  label in position [[3,1]] and value in [[3,2]].

xllcorner=headerData[[3,2]];

"yllcorner" is the UTM zone 11 "y" coordinate, in meters, of the lower left corner of the DEM;  label in position [[4,1]] and value in [[4,2]]

yllcorner=headerData[[4,2]];

"cellsize" a.k.a. "grid spacing" is the distance between grid nodes in the DEM file, in meters;  label in position [[5,1]] and value in [[5,2]].

cellsize=headerData[[5,2]];

"nodata_value" is the value used in the data file to indicate that there is no data in a specific location in the data file;  label in position [[6,1]] and value in [[6,2]]

nodataValue=headerData[[6,2]];

The seventh row of the data set is the first data row.

The true UTM coordinates of the datum in position [[i,j]] of the data set (from row 7 to the end) are as follows:

x coordinate = xllcorner + (cellsize * (j - 1))          y coordinate = yllcorner + (cellsize * (nrows - i))          
           z coordinate = the datum in position [[i,j]]

nullData=0.;

Convert the input coordinates of the earthquake focus to the same coordinate system as the DEM data

utmCoordinates=convertToUTM[focalLat,focalLong,zoneMeridian];

The "focus" is the (x, y, z) coordinates in meters of the reported earthquake focus, in a coordinate system in which the origin is the lower-left (southwest) corner of the DEM at sea level.

focus={(utmCoordinates[[1]]-xllcorner),(utmCoordinates[[2]]-yllcorner),(focalDepth)};

Define the threshold for resolving whether a point in the DEM is within the fault swath

The "zoneHalfWidth" is half the distance between the parallel planes that bound the uncertainty region for this fault-plane solution.  The units are meters.

zoneHalfWidth=(vertError*Cos[faultDipPlunge Degree])+(horizError*Sin[faultDipPlunge Degree]);

Define the location of the seismo-lineament

Determine which points on the DEM lie within the uncertainty swath for the fault location, ± the width of the "zoneHalfWidth," and map the corresponding boundaries

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focus[[1]]),((cellsize*(nrows-i))-focus[[2]]),(mydata[[i+6,j]]-focus[[3]]),zoneHalfWidth,faultNormalVect,nullData],{i,nrows},{j,ncols}];

elev21=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

lightDirection1=(faultDipTrend*(π/180))+(π/2);

lightDirection2=lightDirection1+π;

Several options for color gradients exist for use with "ColorFunction" including "GreenBrownTerrain", "AlpineColors", "ArmyColors", "GrayTones", "SandyTerrain", "DarkTerrain", and  "CoffeeTones".  These color schemes can be recovered using the function ColorData.

finalImage3=ReliefPlot[elev21,ColorFunction->"GreenBrownTerrain",LightingAngle->{lightDirection1,Pi/12}]

finalImage4=ReliefPlot[elev21,ColorFunction->"GreenBrownTerrain",LightingAngle->{lightDirection2,Pi/12}]

Output

Export data and image files

The file created in the next line is JPEG image that contains the graphic showing the hillshade image within the seismo-lineament, illuminated from the dip direction and 180° to the dip direction of the fault plane.

IMPORTANT NOTE:  The user must supply a name for the output file that is different from existing file names, or else the existing files may be over-written.

Export["Ariz1959Swath1a.jpg", finalImage3,ImageResolution->300];

Export["Ariz1959Swath1b.jpg", finalImage4,ImageResolution->300];

Before we end the program, we do some housecleaning by clearing the data values from the variable names (or list names).

ClearAll[mergedData,mergedSwathData,traceImageFile,traceImageFile1,traceImageFile2,finalImage1,finalImage2,mydata,inFaultTrace,uncertSwath,elev21];

How long did this program take to run, in minutes?

minutesForProcessing=(AbsoluteTime[]-startTime)/60

References

Colley, S.J., 2002, Vector calculus [2nd edition]:  Upper Saddle River, New Jersey, Prentice-Hall, 558 p., (pp. 43-50), ISBN 0-13-041531-6.

Cronin, V.S., and Sverdrup, K.A., 1998, Preliminary assessment of the seismicity of the Malibu Coast Fault Zone, southern California, and related issues of philosophy and practice, in Welby, C.W., and Gowan, M.E. [editors], A Paradox of Power--Voices of Warning and Reason in the Geosciences:  Geological Society of America, Reviews in Engineering Geology, p. 123-155.

Cronin, V.S., Byars, B.W., and Gammill, T., 2003, Developing techniques for regional structural interpretation using GIS and DEM-based terrain analysis:  Geological Society of America, Abstracts with Programs, v. 34, no. 7, p. 261.  http://gsa.confex.com/gsa/2003AM/finalprogram/abstract_61851.htm

Cronin, V.S., Millard, M.A., Seidman, L.E., and Bayliss, B.G., 2008, The Seismo-Lineament Analysis Method [SLAM] -- A reconnaissance tool to help find seismogenic faults: Environmental and Engineering Geoscience, v. 14, no. 3, p. 199-219.

Davis, H.F., and Snider, A.D., 1987, Introduction to vector analysis [5th edition]:  Boston, Allyn and Bacon, 365 p. (pp. 34-37), ISBN 0-205-10263-8.

Gammill, T., Cronin, V.S., and Byars, B.W., 2004, Combining earthquake focal data and digital map analysis in reconnaissance for active faults, central Santa Monica Mountains and northern Santa Monica Bay, California:  Geological Society of America, Abstracts with Programs,v. 36, no. 5, http://gsa.confex.com/gsa/2004AM/finalprogram/abstract_80352.htm

Haneberg, W.C., 2004, Computational geosciences with Mathematica:  Berlin, Springer-Verlag, 381 p., CD, ISBN 3-540-40245-4.

Snyder, J.P., 1982, Map projections used by the U.S. Geological Survey:  U.S. Geological Survey Bulletin 1532, pp. 63-69 and 233-235.