| Title: | Graphics for Spherical Distributions and Earthquake Focal Mechanisms |
|---|---|
| Description: | Graphics for statistics on a sphere, as applied to geological fault data, crystallography, earthquake focal mechanisms, radiation patterns, ternary plots and geographical/geological maps. Non-double couple plotting of focal spheres and source type maps are included for statistical analysis of moment tensors. |
| Authors: | Jonathan M. Lees [aut, cre], Keehoon Kim [ctb] |
| Maintainer: | Jonathan M. Lees <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 3.4-10 |
| Built: | 2025-01-30 03:36:37 UTC |
| Source: | https://github.com/cran/RFOC |
Graphics for statistics on a sphere, as applied to geological fault data, crystallography, earthquake focal mechanisms, radiation patterns, ternary plots and geographical/geological maps. Given strike-dip-rake or a set of fault planes, focal planes, RFOC creates structures for manipulating and plotting earthquake focal mechanisms as individual plots or distributed spatially maps.
RFOC can be used for analysis of plane orientation, geologic structure, distribution of stress and strain analyses.
Visualize focal mechanisms in a number of modes, including: beachball plots, radiation plots, fault planes and ternary diagrams. Shows spatial distribution of spherically distributed data.
Jonathan M. Lees Maintainer: Jonathan M. Lees <[email protected]>
J. M. Lees. Geotouch: Software for three and four dimensional GIS in the earth sciences. Computers and Geosciences , 26(7):751–761, 2000.
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
Snyder, John P., 1987, Map Projections-a working manual, USGS-Professional Paper, 383p.
C. Frohlich. Triangle diagrams: ternary graphs to display similarity and diversity of earthquake focal mechanisms. Physics of the Earth and Planetary Interiors, 75:193-198, 1992.
RSEIS, GEOmap, zoeppritz
############# plot one focal mechanism: M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=TRUE) ############# plot many P-axes: paz = rnorm(100, mean=297, sd=100) pdip = rnorm(100, mean=52, sd=20) net() focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) ############# #### Show many Focal mechanisms on a plot: Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) plot(MZ[,1], MZ[,2], type='n', xlab="LON", ylab="LAT", asp=1) for(i in 1:length(MZ[,1])) { paste(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = SDRfoc(MZ[i,3], MZ[i,4], MZ[i,5], u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) fcol = foc.color(foc.icolor(MEC$rake1), pal=1) justfocXY(MEC, x=MZ[i,1], y =MZ[i,2] , focsiz=0.5, fcol =fcol , fcolback = "white", xpd = TRUE) }############# plot one focal mechanism: M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=TRUE) ############# plot many P-axes: paz = rnorm(100, mean=297, sd=100) pdip = rnorm(100, mean=52, sd=20) net() focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) ############# #### Show many Focal mechanisms on a plot: Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) plot(MZ[,1], MZ[,2], type='n', xlab="LON", ylab="LAT", asp=1) for(i in 1:length(MZ[,1])) { paste(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = SDRfoc(MZ[i,3], MZ[i,4], MZ[i,5], u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) fcol = foc.color(foc.icolor(MEC$rake1), pal=1) justfocXY(MEC, x=MZ[i,1], y =MZ[i,2] , focsiz=0.5, fcol =fcol , fcolback = "white", xpd = TRUE) }
Add a standard set of points to a Focal Mechanism
addmecpoints(MEC, pch = 5)addmecpoints(MEC, pch = 5)
MEC |
MEC structure list |
pch |
plotting character |
Graphical Side effects
Jonathan M. Lees<[email protected]>
SDRfoc, focpoint
MEC= SDRfoc(12,34,-120) addmecpoints(MEC)MEC= SDRfoc(12,34,-120) addmecpoints(MEC)
Add Pressure and tension Axes to focal mechanism
addPT(MEC, pch = 5)addPT(MEC, pch = 5)
MEC |
MEC structure |
pch |
plotting character |
Graphical Side Effect
Jonathan M. Lees<[email protected]>
addPTarrows
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Beachfoc(MEC) addPT(MEC, pch = 5)MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Beachfoc(MEC) addPT(MEC, pch = 5)
Illustrate Pressure and Tension axis on Focal Plot using 3D arrows
addPTarrows(MEC)addPTarrows(MEC)
MEC |
Mechanism Structure |
Graphical Side Effects
This function looks better when plotting the upper hemisphere
Jonathan M. Lees<[email protected]>
focpoint, BOXarrows3D,Z3Darrow
MEC = SDRfoc(65,25,13, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=TRUE) addPTarrows(MEC)MEC = SDRfoc(65,25,13, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=TRUE) addPTarrows(MEC)
Calculate and plot small circle on Stereo net at arbitrary azimuth, orientation and conical angle
addsmallcirc(az, iang, alphadeg, BALL.radius = 1, N = 100, add = TRUE, ...)addsmallcirc(az, iang, alphadeg, BALL.radius = 1, N = 100, add = TRUE, ...)
az |
Azimuth of axis |
iang |
angle of dip, degrees |
alphadeg |
width of cone in degrees |
BALL.radius |
size of sphere |
N |
NUmber of points to calculate |
add |
logical, TRUE=add to existing plot |
... |
graphical parameters |
Given the azimuth and dip of a vector, plot the small circle around the pole with conical angle alphadeg
LIST:
x |
x-coordinates |
y |
y-coordinates |
alphadeg is the radius of the conic projection
Jonathan M. Lees<[email protected]>
net
net() addsmallcirc(65, 13, 20, BALL.radius = 1, N = 100, add = TRUE) addsmallcirc(165, 73, 5.6, BALL.radius = 1, N = 100, add = TRUE)net() addsmallcirc(65, 13, 20, BALL.radius = 1, N = 100, add = TRUE) addsmallcirc(165, 73, 5.6, BALL.radius = 1, N = 100, add = TRUE)
Using a Starting LAT-LON, return points along an azimuth
AlongGreat(LON1, LAT1, km1, ang, EARTHRAD= 6371)AlongGreat(LON1, LAT1, km1, ang, EARTHRAD= 6371)
LON1 |
Longitude, point |
LAT1 |
Latitude, point |
km1 |
Kilometers in direction ang |
ang |
Direction from North |
EARTHRAD |
optional earth radius, default = 6371 |
Returns LAT-LON points along a great circle, so many kilometers away in a specified direction
LIST:
lat |
Latitude, destination point |
lon |
Longitude, destination point |
distdeg |
distance in degrees |
distkm |
distance in km |
Jonathan M. Lees <[email protected]>
london = c(51.53333, -0.08333333) AlongGreat(london[2], london[1], 450, 56)london = c(51.53333, -0.08333333) AlongGreat(london[2], london[1], 450, 56)
Calculates conical projection angle for 95% confidence bounds for mean of spherically distributed data.
alpha95(az, iang)alpha95(az, iang)
az |
vector of azimuths, degrees |
iang |
vector of dips, degrees |
Program calculates the cartesian coordinates of all poles, sums and returns the resultant vector, its azimuth and length (R). For N points, statistics include:
where 's are the relevant eigenvalues of matrix MAT and
angles are in degrees.
LIST:
Ir |
resultant inclination, degrees |
Dr |
resultant declination, degrees |
R |
resultant sum of vectors, normalized |
K |
K-dispersion value |
S |
spherical variance |
Alph95 |
95% confidence angle, degrees |
Kappa |
log ratio of eignevectors |
E |
Eigenvactors |
MAT |
matrix of cartesian vectors |
Jonathan M. Lees<[email protected]>
Davis, John C., 2002, Statistics and data analysis in geology, Wiley, New York, 637p.
addsmallcirc
paz = rnorm(100, mean=297, sd=10) pdip = rnorm(100, mean=52, sd=8) ALPH = alpha95(paz, pdip) ######### draw stereonet net() ############ add points focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) ############### add 95 percent confidence bounds addsmallcirc(ALPH$Dr, ALPH$Ir, ALPH$Alph95, BALL.radius = 1, N = 25, add = TRUE, lwd=1, col='blue') ############ second example: paz = rnorm(100, mean=297, sd=100) pdip = rnorm(100, mean=52, sd=20) ALPH = alpha95(paz, pdip) net() focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) addsmallcirc(ALPH$Dr, 90-ALPH$Ir, ALPH$Alph95, BALL.radius = 1, N = 25, add = TRUE, lwd=1, col='blue')paz = rnorm(100, mean=297, sd=10) pdip = rnorm(100, mean=52, sd=8) ALPH = alpha95(paz, pdip) ######### draw stereonet net() ############ add points focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) ############### add 95 percent confidence bounds addsmallcirc(ALPH$Dr, ALPH$Ir, ALPH$Alph95, BALL.radius = 1, N = 25, add = TRUE, lwd=1, col='blue') ############ second example: paz = rnorm(100, mean=297, sd=100) pdip = rnorm(100, mean=52, sd=20) ALPH = alpha95(paz, pdip) net() focpoint(paz, pdip, col='red', pch=3, lab="", UP=FALSE) addsmallcirc(ALPH$Dr, 90-ALPH$Ir, ALPH$Alph95, BALL.radius = 1, N = 25, add = TRUE, lwd=1, col='blue')
Interactive extract axis point on Stereonet
AXpoint(UP = TRUE, col=2, n=1)AXpoint(UP = TRUE, col=2, n=1)
UP |
logical, TRUE=upper hemisphere |
col |
plotting color |
n |
maximum number to locate, default=unlimited |
Program uses locator to create a vector of poles. Points outside the focal sphere (r>1) are ignored. If n is missing, locator continues until stopped (middle mouse in linux, stop in windows).
phiang |
azimuth angle, degrees |
dip |
dip angle, degrees |
x |
x-coordinate of cartesian vector |
y |
y-coordinate of cartesian vector |
z |
z-coordinate of cartesian vector |
gx |
x-coordinate of prjection |
gy |
y-coordinate of prjection |
Jonathan M. Lees<[email protected]>
locator, qpoint, EApoint
#################### this is interactive ## Not run: net() Z = AXpoint(UP = TRUE) ## click in steronet Z ## End(Not run)#################### this is interactive ## Not run: net() Z = AXpoint(UP = TRUE) ## click in steronet Z ## End(Not run)
Calculates the angle between two 2D normalized vectors using dot and cross product
bang(x1, y1, x2, y2)bang(x1, y1, x2, y2)
x1 |
x coordinate of first normalized vector |
y1 |
y coordinate of first normalized vector |
x2 |
x coordinate of second normalized vector |
y2 |
y coordinate of second normalized vector |
The sign of angle is determined by the sign of the cross product of the two vectors.
angle in radians
Vectors must be normalized prior to calling this routine. Used mainly for vectors on the unit sphere.
Jonathan M. Lees <[email protected]>
v1 = c(5,3) v2 = c(6,1) a1 = c(5,3)/sqrt(v1[1]^2+v1[2]^2) a2 = c(6,1)/sqrt(v2[1]^2+v2[2]^2) plot(c(0, v1[1],v2[1] ) , c(0, v1[2],v2[2]), type='n', xlab="x", ylab="y" ) text(c(v1[1],v2[1]) , c(v1[2],v2[2]), labels=c("v1", "v2"), pos=3, xpd=TRUE) arrows(0, 0, c(v1[1],v2[1] ), c(v1[2],v2[2])) B = 180*bang(a1[1], a1[2], a2[1], a2[2])/pi title(paste(sep=" ", "Angle from V1 to V2=",format(B, digits=2)) )v1 = c(5,3) v2 = c(6,1) a1 = c(5,3)/sqrt(v1[1]^2+v1[2]^2) a2 = c(6,1)/sqrt(v2[1]^2+v2[2]^2) plot(c(0, v1[1],v2[1] ) , c(0, v1[2],v2[2]), type='n', xlab="x", ylab="y" ) text(c(v1[1],v2[1]) , c(v1[2],v2[2]), labels=c("v1", "v2"), pos=3, xpd=TRUE) arrows(0, 0, c(v1[1],v2[1] ), c(v1[2],v2[2])) B = 180*bang(a1[1], a1[2], a2[1], a2[2])/pi title(paste(sep=" ", "Angle from V1 to V2=",format(B, digits=2)) )
Plots a focal mechanism in beachball style
Beachfoc(MEC, fcol = gray(0.9), fcolback = "white", ALIM = c(-1, -1, +1, +1))Beachfoc(MEC, fcol = gray(0.9), fcolback = "white", ALIM = c(-1, -1, +1, +1))
MEC |
Mechanism Structure |
fcol |
color for the filled portion of the beachball |
fcolback |
color for the background portion of the beachball, default='white' |
ALIM |
Bounding box for beachball, default=c(-1, -1, +1, +1) |
Beachfoc is run after MEC is set using SDRfoc. Options for plotting the beachball in various modes are controlled by flags set in MEC
Used for its graphical side effect
Jonathan M. Lees <[email protected]>
K. Aki and P. G. Richards. Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002. Keiiti Aki, Paul G. Richards. ill. ; 26 cm.
CONVERTSDR, SDRfoc, justfocXY
MEC = SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=TRUE) Beachfoc(MEC, fcol=MEC$fcol, fcolback="white")MEC = SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=TRUE) Beachfoc(MEC, fcol=MEC$fcol, fcolback="white")
Calculates Angles for determining ternary distribution of faults based on P-T axis orientation.
Bfocvec(Paz, Pdip, Taz, Tdip)Bfocvec(Paz, Pdip, Taz, Tdip)
Paz |
vector of azimuths, degrees |
Pdip |
vector of dips, degrees |
Taz |
vector of azimuths, degrees |
Tdip |
vector of dips, degrees |
This calculation is based on Froelich's paper.
LIST:
Bdip |
azimuths, degrees |
Baz |
dips, degrees |
Jonathan M. Lees<[email protected]>
C. Frohlich. Triangle diagrams: ternary graphs to display similarity and diversity of earthquake focal mechanisms. Physics of the Earth and Planetary Interiors, 75:193-198, 1992.
ternfoc.point
Msdr = CONVERTSDR(55.01, 165.65, 29.2 ) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td )Msdr = CONVERTSDR(55.01, 165.65, 29.2 ) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td )
Create and project and plot 3D arrows with viewing Matrix.
BOXarrows3D(x1, y1, z1, x2, y2, z2, aglyph = NULL, Rview = ROTX(0), col = grey(0.5), border = "black", len = 0.7, basethick = 0.05, headlen = 0.3, headlip = 0.02)BOXarrows3D(x1, y1, z1, x2, y2, z2, aglyph = NULL, Rview = ROTX(0), col = grey(0.5), border = "black", len = 0.7, basethick = 0.05, headlen = 0.3, headlip = 0.02)
x1 |
x-coordinates of base of arrows |
y1 |
y-coordinates of base of arrows |
z1 |
z-coordinates of base of arrows |
x2 |
x-coordinates of head of arrows |
y2 |
y-coordinates of head of arrows |
z2 |
z-coordinates of head of arrows |
aglyph |
glyph structure, default is Z3Darrow |
Rview |
Viewing matrix |
col |
fill color |
border |
Border color |
len |
Length |
basethick |
thickness of the base |
headlen |
thickness of the head |
headlip |
width of the overhanging lip |
Arrows point from base to head.
Used for graphical side effects.
Any 3D glyph strucutre can be used
Jonathan M. Lees <[email protected]>
Z3Darrow
## Not run: #### animate 10 random arrow vectors L = list(x1 = runif(10, min=-2, max=2), y1 = runif(10, min=-2, max=2), z1=runif(10, min=-4, max=4), x2 = runif(10, min=-2, max=2), y2 = runif(10, min=-2, max=2), z2=runif(10, min=-4, max=4) ) headlen = .3 len = .7 basethick = 0.05 headlip = .02 aglyph = Z3Darrow(len = len , basethick =basethick , headlen =headlen , headlip=headlip ) r1 = 8 theta = seq(from=0, to=2*360, length=200) mex = r1*cos(theta*pi/180) mey = r1*sin(theta*pi/180) mez = seq(from=r1, to =0 , length=length(mex)) ## mez=rep(r1, length=length(mex)) angz = atan2(mey, mex)*180/pi angx = atan2(sqrt(mex^2+mey^2), mez)*180/pi pal=c("red", "blue", "green") ## aglyph = gblock for(j in 1:length(angz)) { Rview = ROTZ(angz[j]) plot(c(-4,4), c(-4,4), type='n', asp=1); grid() BOXarrows3D(L$x1,L$y1,L$z1, L$x2,L$y2,L$z2, aglyph=aglyph, Rview=Rview, col=pal) Sys.sleep(.1) } ## End(Not run)## Not run: #### animate 10 random arrow vectors L = list(x1 = runif(10, min=-2, max=2), y1 = runif(10, min=-2, max=2), z1=runif(10, min=-4, max=4), x2 = runif(10, min=-2, max=2), y2 = runif(10, min=-2, max=2), z2=runif(10, min=-4, max=4) ) headlen = .3 len = .7 basethick = 0.05 headlip = .02 aglyph = Z3Darrow(len = len , basethick =basethick , headlen =headlen , headlip=headlip ) r1 = 8 theta = seq(from=0, to=2*360, length=200) mex = r1*cos(theta*pi/180) mey = r1*sin(theta*pi/180) mez = seq(from=r1, to =0 , length=length(mex)) ## mez=rep(r1, length=length(mex)) angz = atan2(mey, mex)*180/pi angx = atan2(sqrt(mex^2+mey^2), mez)*180/pi pal=c("red", "blue", "green") ## aglyph = gblock for(j in 1:length(angz)) { Rview = ROTZ(angz[j]) plot(c(-4,4), c(-4,4), type='n', asp=1); grid() BOXarrows3D(L$x1,L$y1,L$z1, L$x2,L$y2,L$z2, aglyph=aglyph, Rview=Rview, col=pal) Sys.sleep(.1) } ## End(Not run)
Draw circular ticmarks
circtics(r = 1, dr = 0.02, dang = 10, ...)circtics(r = 1, dr = 0.02, dang = 10, ...)
r |
radius |
dr |
length of tics |
dang |
angle between tics |
... |
graphical parameters |
graphical side effects
Jonathan M. Lees <[email protected]>
phi = seq(from =0, to = 2 * pi, length=360) x = cos(phi) y = sin(phi) plot(x, y, col = 'blue', asp=1, type='l') circtics(r = 1, dr = 0.02, dang = 10, col='red')phi = seq(from =0, to = 2 * pi, length=360) x = cos(phi) y = sin(phi) plot(x, y, col = 'blue', asp=1, type='l') circtics(r = 1, dr = 0.02, dang = 10, col='red')
Takes Strike-Dip-Rake and creates planes and pole locations for MEC structure
CONVERTSDR(strike, dip, rake)CONVERTSDR(strike, dip, rake)
strike |
angle, degrees, strike of down dip directin |
dip |
angle, degrees, dip is measured from the horizontal NOT from the NADIR |
rake |
angle, degrees |
input is strike dip and rake in degrees
LIST:
strike |
strike |
dipdir |
dip |
rake |
rake |
F |
list(az, dip) of F-pole |
G |
list(az, dip) of G-pole |
U |
list(az, dip) of U-pole |
V |
list(az, dip) of V-pole |
P |
list(az, dip) of P-pole |
T |
list(az, dip) of T-pole |
M |
list( az1=0, d1=0, az2=0, d2=0, uaz=0, ud=0, vaz=0, vd=0, paz=0, pd =0, taz=0, td=0) |
Jonathan M. Lees <[email protected]>
BeachFoc
s=65 d=25 r=13 mc = CONVERTSDR(s,d,r )s=65 d=25 r=13 mc = CONVERTSDR(s,d,r )
Vector Cross Product with list as arguments and list as values
cross.prod(B, A)cross.prod(B, A)
B |
list of x,y,z |
A |
list of x,y,z |
LIST
x, y, z
|
vector of cross product |
Jonathan M. Lees<[email protected]>
RSEIS::xprod
B1 = list(x=4, y=9, z=2) B2 = list(x=2,y=-5,z=4) cross.prod(B1, B2)B1 = list(x=4, y=9, z=2) B2 = list(x=2,y=-5,z=4) cross.prod(B1, B2)
returns cross product of two vectors in list format
CROSSL(A1, A2)CROSSL(A1, A2)
A1 |
list x,y,z |
A2 |
list x,y,z |
List
x, y, z
|
input vector |
az |
azimuth, degrees |
dip |
dip, degrees |
Jonathan M. Lees<[email protected]>
RSEIS::xprod
A1 = list(x=1,y=2, z=3) A2 = list(x=12,y=-2, z=-5) N = CROSSL(A1, A2)A1 = list(x=1,y=2, z=3) A2 = list(x=12,y=-2, z=-5) N = CROSSL(A1, A2)
Plot Non-double Couple Moment
doNonDouble(moments, sel = 1, col=rgb(1, .75, .75))doNonDouble(moments, sel = 1, col=rgb(1, .75, .75))
moments |
list of moments: seven elements. See details. |
sel |
integer vector, index of moments to plot |
col |
color, either a single color, rgb, or a color palette |
Plot, sequentially the moments using the CLVD (non-double couple component. The first element of the list is the integer index of the event. The next six elements are the moments in the following order, c(Mxx, Myy, Mzz, Mzy, Mxz, Mxy) .
If the data is in spherical coordinates, one must switch the
sign of the Mrp and Mtp components, so:
Mrr = Mzz
Mtt = Mxx
Mpp = Myy
Mrt = Mxz
Mrp = -Myz
Mtp = -Mxy
Side effects
If events are read in using spherical rather than cartesian
coordinates
need a conversion:
Mrr = Mzz
Mtt = Mxx
Mpp = Myy
Mrt = Mxz
Mrp = -Myz
Mtp = -Mxy
Jonathan M. Lees<[email protected]>
Ekstrom, G.; Nettles, M. and DziewoDski, A. The Global CMT Project 2004-2010: centroid-moment tensors for 13,017 earthquakes Physics of the Earth and Planetary Interiors, 2012.
MapNonDouble, ShadowCLVD, angles, nodalLines, PTaxes
mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)
Interactive locator to calculate x,y orientation, dip coordinates and plots on an equalarea stereonet
EApoint()EApoint()
Used for returning a set of strike/dip angles on Equal-area stereonet plot.
LIST:
phi |
orientation, degrees |
iang |
angle of dip, degrees |
x |
x-coordinate |
y |
y-coordinate |
Jonathan M. Lees<[email protected]>
qpoint, focpoint
#################### this is interactive ### collect points with locator() ## Not run: net() eps = EApoint() ### plot results net() qpoint(eps$phi , eps$iang) ## End(Not run)#################### this is interactive ### collect points with locator() ## Not run: net() eps = EApoint() ### plot results net() qpoint(eps$phi , eps$iang) ## End(Not run)
Cartesian moment tensors from Tungurahua Volcano, Ecuador
data(egl)data(egl)
A list of 84 moment tensors, each elelment consists of: lam1, lam2, lam3, vec1, vec2,vec3, ratio, force.
See below
Kim, K., Lees, J.M. and Ruiz, M., (2014) Source mechanism of Vulcanian eruption at Tungurahua Volcano, Ecuador, derived from seismic moment tensor inversions, J. Geophys. Res., February, 2014. Vol. 119(2): pp. 1145-1164.
data(egl) typl1=c(2,4,7,12,13,16,17,18,19,20,24,25,26,27, 28,29,30,31,33,34,35,36,37,38,40,43,50, 59,62,73,74,77,8,79,80,81,83,84) typl2=c(5,6,8,9,10,11,14,15,22,42,46,47,48,49, 51,52,53,54,55,56,57,58,60,61,63,72,82) evtns=1:84 par(mfrow=c(1,2)) T1 = TapeBase() TapePlot(T1) for(i in 1:length(egl)) { i1 = egl[[i]] E1 = list(values=c(i1$lam1, i1$lam2, i1$lam3), vectors = cbind(i1$vec1, i1$vec2, i1$vec3)) testrightHAND(E1$vectors) E1$vectors = forcerighthand(E1$vectors) mo=sort(E1$values,decreasing=TRUE) # M=sum(mo)/3 # Md=mo-M h = SourceType(mo) h$dip = 90-h$phi h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) if(i %in% typl1) { col="red" }else{col="blue" } points(h1$x, h1$y, pch=21, bg=col ) } par(mai=c(0,0,0,0)) hudson.net() for(i in 1:length(typl1)) { egv=egl[[typl1[i]]] m=c(egv$lam1,egv$lam2,egv$lam3) col='red' hudson.plot(m=m,col=col) } for(i in 1:length(typl2)) { egv=egl[[typl2[i]]] m=c(egv$lam1,egv$lam2,egv$lam3) col='blue' hudson.plot(m=m,col=col,lwd=2) }data(egl) typl1=c(2,4,7,12,13,16,17,18,19,20,24,25,26,27, 28,29,30,31,33,34,35,36,37,38,40,43,50, 59,62,73,74,77,8,79,80,81,83,84) typl2=c(5,6,8,9,10,11,14,15,22,42,46,47,48,49, 51,52,53,54,55,56,57,58,60,61,63,72,82) evtns=1:84 par(mfrow=c(1,2)) T1 = TapeBase() TapePlot(T1) for(i in 1:length(egl)) { i1 = egl[[i]] E1 = list(values=c(i1$lam1, i1$lam2, i1$lam3), vectors = cbind(i1$vec1, i1$vec2, i1$vec3)) testrightHAND(E1$vectors) E1$vectors = forcerighthand(E1$vectors) mo=sort(E1$values,decreasing=TRUE) # M=sum(mo)/3 # Md=mo-M h = SourceType(mo) h$dip = 90-h$phi h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) if(i %in% typl1) { col="red" }else{col="blue" } points(h1$x, h1$y, pch=21, bg=col ) } par(mai=c(0,0,0,0)) hudson.net() for(i in 1:length(typl1)) { egv=egl[[typl1[i]]] m=c(egv$lam1,egv$lam2,egv$lam3) col='red' hudson.plot(m=m,col=col) } for(i in 1:length(typl2)) { egv=egl[[typl2[i]]] m=c(egv$lam1,egv$lam2,egv$lam3) col='blue' hudson.plot(m=m,col=col,lwd=2) }
Create and plot fancy arrows. Aspect ratio must be set to 1-1 for these arrows to plot correctly.
fancyarrows(x1, y1, x2, y2, thick = 0.08, headlength = 0.4, headthick = 0.2, col = grey(0.5), border = "black")fancyarrows(x1, y1, x2, y2, thick = 0.08, headlength = 0.4, headthick = 0.2, col = grey(0.5), border = "black")
x1 |
x tail coordinate |
y1 |
y tail coordinate |
x2 |
x head coordinate |
y2 |
y head coordinate |
thick |
thickness of arrow |
headlength |
length of head |
headthick |
thickness of head |
col |
fill color |
border |
color of border |
Graphical side effects.
fancyarrows only work if te aspect ratio is set to 1. See example below.
Jonathan M. Lees<[email protected]>
TEACHFOC
thick = 0.01; headlength = 0.2; headthick = 0.1 x = runif(10, -1, 1) y = runif(10, -1, 1) ############ MUST set asp=1 here plot(x,y, asp=1) fancyarrows(rep(0, 10) , rep(0, 10) ,x, y, thick =thick , headlength = headlength, headthick =headthick)thick = 0.01; headlength = 0.2; headthick = 0.1 x = runif(10, -1, 1) y = runif(10, -1, 1) ############ MUST set asp=1 here plot(x,y, asp=1) fancyarrows(rep(0, 10) , rep(0, 10) ,x, y, thick =thick , headlength = headlength, headthick =headthick)
given azimuth and dip of fault mechanism, calculate and plot the fault plane.
faultplane(az, dip, col = par("col"), PLOT = TRUE, UP = FALSE,lwd=2, lty=1, ...)faultplane(az, dip, col = par("col"), PLOT = TRUE, UP = FALSE,lwd=2, lty=1, ...)
az |
degrees, strike of the plane (NOT down dip azimuth) |
dip |
degrees, dip from horizontal |
col |
color for line |
PLOT |
option for adding to plot |
UP |
upper or lower hemisphere |
lwd |
Line Width |
lty |
Line Type |
... |
graphical parameters |
Azimuth is the strike in degrees, not the down dip azimuth as described in other routines.
list of points along fault plane
x |
coordinates on focal sphere |
y |
coordinates on focal sphere |
Jonathan M. Lees <[email protected]>
Beachfoc
gcol='black' border='black' ndiv=36 phi = seq(0,2*pi, by=2*pi/ndiv); x = cos(phi); y = sin(phi); plot(x,y, type='n', asp=1) lines(x,y, col=border) lines(c(-1,1), c(0,0), col=gcol) lines(c(0,0), c(-1,1), col=gcol) faultplane(65, 34)gcol='black' border='black' ndiv=36 phi = seq(0,2*pi, by=2*pi/ndiv); x = cos(phi); y = sin(phi); plot(x,y, type='n', asp=1) lines(x,y, col=border) lines(c(-1,1), c(0,0), col=gcol) lines(c(0,0), c(-1,1), col=gcol) faultplane(65, 34)
Fix az, dip angles so they fall in correct quadrant.
FixDip(A)FixDip(A)
List:
A |
|
Quadrants are determined by the sine and cosine of the dip angle:
co = cos(dip)
si = sin(dip)
quad[co>=0 & si>=0] = 1
quad[co<0 & si>=0] = 2
quad[co<0 & si<0] = 3
quad[co>=0 & si<0] = 4
List:
az |
azimuthm angle, degrees |
dip |
dip angle, degrees |
Jonathan M. Lees<[email protected]>
RPMG::fmod
B = list(az=231, dip = -65) FixDip(B)B = list(az=231, dip = -65) FixDip(B)
Switch a focal mechanism so the auxilliary plane is the nodal plane.
flipnodal(s1, d1, r1)flipnodal(s1, d1, r1)
s1 |
Strike |
d1 |
Dip |
r1 |
Rake |
Fuunction is used for orienting a set of fault planes to line up according to a geologic interpretation.
List:
s1 |
Strike |
d1 |
Dip |
r1 |
Rake |
Jonathan M. Lees<[email protected]>
s=65 d=25 r=13 mc = CONVERTSDR(s,d,r ) mc2 = flipnodal(s, d, r)s=65 d=25 r=13 mc = CONVERTSDR(s,d,r ) mc2 = flipnodal(s, d, r)
Based on the rake angle, focal styles are assigned an index and assigned a color by foc.color
foc.color(i, pal = 0)foc.color(i, pal = 0)
i |
index to list of focal rupture styles |
pal |
vector of colors |
Since the colors used by focal programs are arbitrary, this routines allows one to change the coloring scheme easily.
foc.icolor returns an index that is used to get the color associated with that style of faulting
Color for plotting, either a name or HEX RGB
Jonathan M. Lees <[email protected]>
foc.icolor
fcolors=c("DarkSeaGreen", "cyan1","SkyBlue1" , "RoyalBlue" ,"GreenYellow","orange","red") foc.color(3, fcolors)fcolors=c("DarkSeaGreen", "cyan1","SkyBlue1" , "RoyalBlue" ,"GreenYellow","orange","red") foc.color(3, fcolors)
Use Rake Angle to determine style of faulting
foc.icolor(rake)foc.icolor(rake)
rake |
degrees, rake angle of fault plane |
The styles are determined by the rake angle
strikeslip abs(rake) <= 15.0 or abs((180.0 - abs(rake))) <= 15.0
rev-obl strk-slp (rake >= 15.0 and rake < 45) or (rake >= 135 and rake < 165)
oblique reverse (rake >= 45.0 and rake < 75) or (rake >= 105 and rake < 135)
reverse rake >= 75.0 and rake < 105.0
norm-oblq strkslp (rake < -15.0 and rake >= -45) or (rake < -135 and rake >= -165)
oblq norm (rake < -45.0 and rake >= -75) or (rake < -105 and rake >= -135)
normal rake < -75.0 and rake >= -105
index (1-6)
Jonathan M. Lees <[email protected]>
foc.color
foc.icolor(25)foc.icolor(25)
Angles for focal planes
FOCangles(m)FOCangles(m)
m |
moment tensor |
Used in MapNonDouble and doNonDouble
vector of 6 angles, 3 for each plane
Lower Hemisphere.
Jonathan M. Lees<[email protected]>
MapNonDouble, doNonDouble, PTaxes, nodalLines
mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) di = dim(moments) number.of.events = di[1] moment_11 = moments[,2] moment_22 = moments[,3] moment_33 = moments[,4] moment_23 = moments[,5] moment_13 = moments[,6] moment_12 = moments[,7] i = 1 m=matrix( c(moment_11[i],moment_12[i],moment_13[i], moment_12[i],moment_22[i],moment_23[i], moment_13[i],moment_23[i],moment_33[i]), ncol=3, byrow=TRUE) angles.all = FOCangles(m) print(angles.all)mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) di = dim(moments) number.of.events = di[1] moment_11 = moments[,2] moment_22 = moments[,3] moment_33 = moments[,4] moment_23 = moments[,5] moment_13 = moments[,6] moment_12 = moments[,7] i = 1 m=matrix( c(moment_11[i],moment_12[i],moment_13[i], moment_12[i],moment_22[i],moment_23[i], moment_13[i],moment_23[i],moment_33[i]), ncol=3, byrow=TRUE) angles.all = FOCangles(m) print(angles.all)
Get character string describing type of fault from its style index
focleg(i)focleg(i)
i |
index to vector of focal styles |
character string used for setting text on plots
String of characters:
Strike slip fault
Reverse Oblique strike-slip fault
Reverse fault
Normal Oblique strike-slip fault
Oblique Normal fault
Formal fault
Jonathan M. Lees <[email protected]>
foc.icolor, foc.color
focleg(2)focleg(2)
Add points on equal-area focal plot
focpoint(az1, dip1, col = 2, pch = 5, lab = "", cex=1, UP = FALSE, PLOT = TRUE, ...)focpoint(az1, dip1, col = 2, pch = 5, lab = "", cex=1, UP = FALSE, PLOT = TRUE, ...)
az1 |
degrees, azimuth angle |
dip1 |
degrees, dip angle |
col |
color |
pch |
plot character for point |
lab |
text lable for point |
cex |
Character Size |
UP |
upper or lower hemisphere |
PLOT |
logical, PLOT=TRUE add points to current plot |
... |
graphical parameters |
List of x,y coordinates on the plot
Jonathan M. Lees <[email protected]>
Beachfoc, addmecpoints
### create focal mech ALIM=c(-1,-1, +1, +1) s=65 d=25 r=13 mc = CONVERTSDR(s,d,r ) MEC = MRake(mc$M) MEC$UP = FALSE MEC$icol = foc.icolor(MEC$rake1) MEC$ileg = focleg(MEC$icol) MEC$fcol = foc.color(MEC$icol) MEC$CNVRG = NA MEC$LIM = ALIM ### plot focal mech Beachfoc(MEC, fcol=MEC$fcol, fcolback="white") ### now add the F anf G axes focpoint(MEC$F$az, MEC$F$dip, pch=5, lab="F", UP=MEC$UP) focpoint(MEC$G$az, MEC$G$dip, pch=5, lab="G", UP=MEC$UP)### create focal mech ALIM=c(-1,-1, +1, +1) s=65 d=25 r=13 mc = CONVERTSDR(s,d,r ) MEC = MRake(mc$M) MEC$UP = FALSE MEC$icol = foc.icolor(MEC$rake1) MEC$ileg = focleg(MEC$icol) MEC$fcol = foc.color(MEC$icol) MEC$CNVRG = NA MEC$LIM = ALIM ### plot focal mech Beachfoc(MEC, fcol=MEC$fcol, fcolback="white") ### now add the F anf G axes focpoint(MEC$F$az, MEC$F$dip, pch=5, lab="F", UP=MEC$UP) focpoint(MEC$G$az, MEC$G$dip, pch=5, lab="G", UP=MEC$UP)
Force Right-Hand System
forcerighthand(U)forcerighthand(U)
U |
3 by 3 matrix |
Flip vectors so they form a right handed system
matrix
Jonathan M. Lees<[email protected]>
testrightHAND
Mtens = c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 = eigen(M1) testrightHAND(E1$vectors) E1$vectors = forcerighthand(E1$vectors) testrightHAND(E1$vectors)Mtens = c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 = eigen(M1) testrightHAND(E1$vectors) E1$vectors = forcerighthand(E1$vectors) testrightHAND(E1$vectors)
Read and reformat CMT solutions downloaded from the web.
getCMT(fn, skip=1)getCMT(fn, skip=1)
fn |
character file name |
skip |
number of lines to skip (e.g. header) |
Data can be extracted from web site: http://www.globalcmt.org/CMTsearch.html
The file must be cleaned prior to scanning - on download from the web site there are extra lines on top and bottom of file. Delete these. Leave one line on the top that describesthe columns. Data is separated by blanks. The files have a mixture of dates - some with 7 component dates (YYMMDD and others with 14 components YYYYMODDHHMM these are read in separately. Missing hours and minutes areset to zero.
list of CMT solution data:
lon |
lon of epicenter |
lat |
lat of epicenter |
str1 |
strike of fault plane |
dip1 |
dip of fault plane |
rake1 |
rake of fault plane |
str2 |
strike of auxilliary plane |
dip2 |
dip of auxilliary plane |
rake2 |
rake of auxilliary plane |
sc |
scale? |
iexp |
exponent? |
name |
name, includes the date |
Elat |
exploding latitude, set to lat initially |
Elon |
exploding longitude, set to lon initially |
jd |
julian day |
yr |
year |
mo |
month |
dom |
day of month |
Use ExplodeSymbols or explode to get new locations for expanding the plotting points.
Jonathan M. Lees<[email protected]>
http://www.globalcmt.org/CMTsearch.html
G. Ekstrom. Rapid earthquake analysis utilizes the internet. Computers in Physics, 8:632-638, 1994.
ExplodeSymbols, spherefocgeo, ternfocgeo
## Not run: g = getCMT("/home/lees/aleut.cmt") pg = prepFOCS(g) plot(range(pg$LONS), range(pg$LATS), type = "n", xlab = "LON", ylab = "LAT", asp = 1) for (i in 1:length(pg$LATS)) { mc = CONVERTSDR(g$str1[i], g$dip1[i], g$rake1[i]) MEC <- MRake(mc$M) MEC$UP = FALSE Fcol <- foc.color(foc.icolor(MEC$rake1), pal = 1) justfocXY(MEC, x = pg$LONS[i], y = pg$LATS[i], focsiz = 0.4, fcol = Fcol, xpd = FALSE) } ## End(Not run)## Not run: g = getCMT("/home/lees/aleut.cmt") pg = prepFOCS(g) plot(range(pg$LONS), range(pg$LATS), type = "n", xlab = "LON", ylab = "LAT", asp = 1) for (i in 1:length(pg$LATS)) { mc = CONVERTSDR(g$str1[i], g$dip1[i], g$rake1[i]) MEC <- MRake(mc$M) MEC$UP = FALSE Fcol <- foc.color(foc.icolor(MEC$rake1), pal = 1) justfocXY(MEC, x = pg$LONS[i], y = pg$LATS[i], focsiz = 0.4, fcol = Fcol, xpd = FALSE) } ## End(Not run)
Calculates rake angles for fault and auxilliary planes
GetRake(az1, dip1, az2, dip2, dir)GetRake(az1, dip1, az2, dip2, dir)
az1 |
azimuth in degrees of fault plane 1 |
dip1 |
dip in degrees of fault plane 1 |
az2 |
azimuth in degrees of auxilliary plane 2 |
dip2 |
dip in degrees of auxilliary plane 2 |
dir |
polarity |
uses output of CONVERTSDR or MEC structure
list of angles for fault plane and auxiallary plane
az1, dip1, rake1, dipaz1
|
strike, dip rake and downdip direction for plane 1 |
az2, dip2, rake2, dipaz2
|
strike, dip rake and downdip direction for plane 2 |
Jonathan M. Lees <[email protected]>
GetRakeSense, CONVERTSDR, Beachfoc, justfocXY
GetRake(345.000000, 25.000000, 122.000000, 71.000000, 1)GetRake(345.000000, 25.000000, 122.000000, 71.000000, 1)
Get the sense of the focal mechanism rake, from the U, V, P, T vectors
GetRakeSense(uaz, upl, vaz, vpl, paz, ppl, taz, tpl)GetRakeSense(uaz, upl, vaz, vpl, paz, ppl, taz, tpl)
uaz |
Azimuth of U vector |
upl |
dip of U vector |
vaz |
Azimuth of V vector |
vpl |
dip of V vector |
paz |
Azimuth of P vector |
ppl |
dip of P vector |
taz |
Azimuth of T vector |
tpl |
dip of T vector |
1, 0 to make sure the region of the T-axis is shaded and the P-axis is blank.
The convention is for the T-axis to be shaded, so this subroutine determines the order of the polygons to be plotted so that the appropriate regins are filled.
Jonathan M. Lees <[email protected]>
GetRake
mc =CONVERTSDR(65,25,13) angsense = GetRakeSense(mc$U$az, mc$U$dip, mc$V$az, mc$V$dip,mc$P$az, mc$P$dip,mc$T$az, mc$T$dip)mc =CONVERTSDR(65,25,13) angsense = GetRakeSense(mc$U$az, mc$U$dip, mc$V$az, mc$V$dip,mc$P$az, mc$P$dip,mc$T$az, mc$T$dip)
Get UW focal mechansims from a file. These are often called A and M cards
getUWfocs(amfile)getUWfocs(amfile)
amfile |
character, file name |
UW focal mechanisms are stored as A and M cards. The A card described the hypocenter the M card describes the focal mechanism.
List:
lon |
numeric, longitude |
lat |
numeric, latitude |
str1 |
numeric, strike of plane 1 |
dip1 |
numeric, dip of plane 1 |
rake1 |
numeric, rake of plane 1 |
str2 |
numeric, strike of plane 2 |
dip2 |
numeric, dip of plane 2 |
rake2 |
numeric, rake of plane 2 |
sc |
character, some GMT info for scale |
iexp |
character, some GMT info for scale |
name |
character, name |
yr |
numeric, year |
mo |
numeric, month |
dom |
numeric, day of month |
jd |
numeric, julian day |
hr |
numeric, hour |
mi |
numeric, minute |
se |
numeric, second |
z |
numeric, depth |
mag |
numeric, magnitude |
Uses UW2 format, so full 4 digit year is required
Jonathan M. Lees<[email protected]>
http://www.unc.edu/~leesj/XM_DOC/xm_hypo.doc.html
getCMT
## Not run: ##### uwpickfile is an ascii format file from University of Washington G1 = getUWfocs(uwpickfile) plot(G1$lon, G1$lat) MEKS = list(lon=G1$lon, lat=G1$lat, str1=G1$str1, dip1=G1$dip1, rake1=G1$rake1, dep=G1$z, name=G1$name) ## utm projection PROJ = GEOmap::setPROJ(type=2, LAT0=mean(G1$lat) , LON0=mean(G1$lon) ) XY = GEOmap::GLOB.XY(G1$lat, G1$lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.05) ## End(Not run)## Not run: ##### uwpickfile is an ascii format file from University of Washington G1 = getUWfocs(uwpickfile) plot(G1$lon, G1$lat) MEKS = list(lon=G1$lon, lat=G1$lat, str1=G1$str1, dip1=G1$dip1, rake1=G1$rake1, dep=G1$z, name=G1$name) ## utm projection PROJ = GEOmap::setPROJ(type=2, LAT0=mean(G1$lat) , LON0=mean(G1$lon) ) XY = GEOmap::GLOB.XY(G1$lat, G1$lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.05) ## End(Not run)
Hammer Equal Area projection
HAMMERprojXY(phi, lam)HAMMERprojXY(phi, lam)
phi |
Latitude, radians |
lam |
Longitude, radians |
list:
x |
coordinate for plotting |
y |
coordinate for plotting |
Jonathan M. Lees<[email protected]>
HAMMERprojXY(-25*pi/180, -16*pi/180)HAMMERprojXY(-25*pi/180, -16*pi/180)
Plot a Hudson plot as preparation for plotting T-k values for focal mechanisms.
hudson.net(add = FALSE, POINTS = TRUE, TEXT = TRUE, colint = "grey", colext = "black")hudson.net(add = FALSE, POINTS = TRUE, TEXT = TRUE, colint = "grey", colext = "black")
add |
logical, TRUE=add to existing plot |
POINTS |
logical, TRUE=add points |
TEXT |
logical, TRUE=add points |
colint |
color for interior lines |
colext |
color for exterior lines |
Draws a T-k plot for moment tensors
Graphical Side effects
Jonathan M. Lees<[email protected]>
Hudson, J.A., Pearce, R.G. and Rogers, R.M., 1989. Source time plot for inversion of the moment tensor, J. Geophys. Res., 94(B1), 765-774.
hudson.plot
hudson.net() Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) hudson.plot(E1$values)hudson.net() Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) hudson.plot(E1$values)
Hudson Source Type Plot
hudson.plot(m, col = "red", pch = 21, lwd = 2, cex = 1, bg="white")hudson.plot(m, col = "red", pch = 21, lwd = 2, cex = 1, bg="white")
m |
vector of eigen values, sorted |
col |
color |
pch |
plotting char |
lwd |
line width |
cex |
character expansion |
bg |
background color for filled symbols |
Add to existing Hudson net
Side effects
Jonathan M. Lees<[email protected]>
Hudson, J.A., Pearce, R.G. and Rogers, R.M., 1989. Source time plot for inversion of the moment tensor, J. Geophys. Res., 94(B1), 765-774.
hudson.net
hudson.net() Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) hudson.plot(E1$values)hudson.net() Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) hudson.plot(E1$values)
Amplitude of P-wave radiation pattern from Double-Couple earthquake
imageP(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)imageP(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)
phiS |
strike |
del |
dip |
lam |
lambda |
SCALE |
logical, TRUE=add scale on side of plot |
UP |
upper/lower hemisphere |
col |
color |
This program calls radP to calculate the radiation pattern and it plots the result using the standard image function
Used for the graphical side effect
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radP, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageP(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageP(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )
add scale to side of an image plot
imageSCALE(z, col, x, y = NULL, size = NULL, digits = 2, labels = c("breaks", "ranges"), nlab = 10)imageSCALE(z, col, x, y = NULL, size = NULL, digits = 2, labels = c("breaks", "ranges"), nlab = 10)
z |
elevation matrix |
col |
palette for plotting |
x |
x location on plot |
y |
y location on plot |
size |
length of scale |
digits |
digits on labels |
labels |
breaks to be plotted |
nlab |
number of breaks to be plotted |
Used for graphical side effect
Jonathan M. Lees <[email protected]>
data(volcano) image(volcano, col=rainbow(100) ) imageSCALE(volcano, rainbow(100), 1.015983, y = 0.874668, size = .01, digits = 2, labels = "breaks", nlab = 20)data(volcano) image(volcano, col=rainbow(100) ) imageSCALE(volcano, rainbow(100), 1.015983, y = 0.874668, size = .01, digits = 2, labels = "breaks", nlab = 20)
Amplitude of SH-wave radiation pattern from Double-Couple earthquake
imageSH(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)imageSH(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)
phiS |
strike |
del |
dip |
lam |
lambda |
SCALE |
logical, TRUE=add scale on side of plot |
UP |
upper/lower hemisphere |
col |
color |
This program calls radP to calculate the radiation pattern and it plots the result using the standard image function
Used for the graphical side effect
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radSH, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageSH(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageSH(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )
Amplitude of SV-wave radiation pattern from Double-Couple earthquake
imageSV(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)imageSV(phiS, del, lam, SCALE = FALSE, UP = FALSE, col = NULL)
phiS |
strike |
del |
dip |
lam |
lambda |
SCALE |
logical, TRUE=add scale on side of plot |
UP |
upper/lower hemisphere |
col |
color |
This program calls radP to calculate the radiation pattern and it plots the result using the standard image function
Used for the graphical side effect
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radSV, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageSV(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) imageSV(MEC$az1, MEC$dip1, MEC$rake1, SCALE=TRUE, UP=MEC$UP, col=rainbow(100) )
Inverse moment tensor from Tape angles.
inverseTAPE(GAMMA, BETA)inverseTAPE(GAMMA, BETA)
GAMMA |
Longitude, degrees |
BETA |
CoLatitude, degrees |
Uses Tape and Tape lune angles to estimate the moment tensor. This function is the inverse of the SourceType calculation. There are two solutions to the systems of equations.
Vectors are scaled by the maximum value.
Moment tensor list:
Va |
vector, First solution |
Vb |
vector, First solution |
The latitude is the CoLatitude.
Either vector can be used as a solution.
Orientation of moment tensor is not preserved int he lune plots.
Jonathan M. Lees<[email protected]>
Tape, W.,and C.Tape(2012), A geometric comparison of source-type plots for moment tensors, Geophys. J. Int., 190, 499-510.
SourceType
lats = seq(from = -80, to = 80, by=10) lons = seq(from=-30, to=30, by=10) i = 3 j = 3 u = inverseTAPE( lons[i], 90-lats[j] )lats = seq(from = -80, to = 80, by=10) lons = seq(from=-30, to=30, by=10) i = 3 j = 3 u = inverseTAPE( lons[i], 90-lats[j] )
Moment Tensors from the Harvard CMT
data(jimbo)data(jimbo)
A list of 9 moment tensors from the Kamchatka region.
http://www.globalcmt.org/CMTsearch.html
Ekstrom, G.; Nettles, M. & DziewoDski, A. The Global CMT Project 2004-2010: centroid-moment tensors for 13,017 earthquakes Physics of the Earth and Planetary Interiors, 2012.
Vertical Rotation matrix
JMAT(phi)JMAT(phi)
phi |
angle, degrees |
First rotate to plan, then within plane rotate to view angle.
3 by 3 matrix
Jonathan M. Lees<[email protected]>
ROTX, ROTZ, ROTY
phi = 18 MAT = JMAT(phi) v1 = c(1,1,0) v2 = MATphi = 18 MAT = JMAT(phi) v1 = c(1,1,0) v2 = MAT
Add simple focal mechanisms to plot
justfocXY(MEC, x = x, y = y, focsiz=1 , fcol = gray(0.9), fcolback = "white", xpd = TRUE)justfocXY(MEC, x = x, y = y, focsiz=1 , fcol = gray(0.9), fcolback = "white", xpd = TRUE)
MEC |
MEC structure |
x |
x-coordinate of center |
y |
y-coordinate of center |
focsiz |
size of focal sphere in inches |
fcol |
color of shaded region |
fcolback |
color of background region |
xpd |
logical, whether to extend the plot beyond, or to clip |
This routine can be used to add focal mechanisms on geographic map or other plot.
Used for graphical side effect
Jonathan M. Lees <[email protected]>
SDRfoc, foc.color
#### read in some data: Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) plot(MZ[,1], MZ[,2], type='n', xlab="LON", ylab="LAT", asp=1) for(i in 1:length(MZ[,1])) { paste(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = SDRfoc(MZ[i,3], MZ[i,4], MZ[i,5], u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) fcol = foc.color(foc.icolor(MEC$rake1), pal=1) justfocXY(MEC, x=MZ[i,1], y =MZ[i,2] , focsiz=.5, fcol =fcol , fcolback = "white", xpd = TRUE) }#### read in some data: Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) plot(MZ[,1], MZ[,2], type='n', xlab="LON", ylab="LAT", asp=1) for(i in 1:length(MZ[,1])) { paste(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = SDRfoc(MZ[i,3], MZ[i,4], MZ[i,5], u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) fcol = foc.color(foc.icolor(MEC$rake1), pal=1) justfocXY(MEC, x=MZ[i,1], y =MZ[i,2] , focsiz=.5, fcol =fcol , fcolback = "white", xpd = TRUE) }
Strike-Dip-Rake and Locations of Harvard CMT catalog for the intersection of the Kamchataka and Aleutian arcs
data(KAMCORN)data(KAMCORN)
The format is: chr "KAMCORN"
The data is selected fromt eh CMT catalog. Parameters are extracted from the normal distribution. Format of the list of data save in KAMCORN is: list(LAT=0 , LON =0 , DEPTH=0 , STRIKE=0 , DIP=0 , RAKE=0 )
http://www.globalcmt.org/CMTsearch.html
G. Ekstrom. Rapid earthquake analysis utilizes the internet. Computers in Physics, 8:632-638, 1994.
data(KAMCORN) plot(KAMCORN$LON, KAMCORN$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1) ###### Paz =vector() Pdip =vector() Taz =vector() Tdip =vector() h = vector() v = vector() IFcol = vector() Fcol = vector() for(i in 1:10) { Msdr = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] ) MEC = MRake(Msdr$M) MEC$UP = FALSE IFcol[i] = foc.icolor(MEC$rake1) Fcol[i] = foc.color(IFcol[i], 1) az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td ) Paz[i] = Msdr$M$paz Pdip[i] = Msdr$M$pd Taz[i] = Msdr$M$taz Tdip[i] = Msdr$M$td h[i] = V$h v[i] = V$v justfocXY( MEC, fcol = Fcol[i], KAMCORN$LON[i], KAMCORN$LAT[i] , focsiz = 0.4 ) }data(KAMCORN) plot(KAMCORN$LON, KAMCORN$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1) ###### Paz =vector() Pdip =vector() Taz =vector() Tdip =vector() h = vector() v = vector() IFcol = vector() Fcol = vector() for(i in 1:10) { Msdr = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] ) MEC = MRake(Msdr$M) MEC$UP = FALSE IFcol[i] = foc.icolor(MEC$rake1) Fcol[i] = foc.color(IFcol[i], 1) az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td ) Paz[i] = Msdr$M$paz Pdip[i] = Msdr$M$pd Taz[i] = Msdr$M$taz Tdip[i] = Msdr$M$td h[i] = V$h v[i] = V$v justfocXY( MEC, fcol = Fcol[i], KAMCORN$LON[i], KAMCORN$LAT[i] , focsiz = 0.4 ) }
takes azimuth and dip and projects the greaat circle on the focla sphere
lowplane(az, dip, col = par("col"), UP = FALSE, PLOT = TRUE)lowplane(az, dip, col = par("col"), UP = FALSE, PLOT = TRUE)
az |
degrees, azimuth of strike of plane |
dip |
degrees, dip |
col |
color of plane |
UP |
upper/lower hemisphere |
PLOT |
add to plot |
Here azimuth is measured from North, and represents the actual strike of the fault line.
list of x,y coordinates of plane
Jonathan M. Lees <[email protected]>
net
net() lowplane(65,23)net() lowplane(65,23)
Moment tensor to T-k
m2tk(m0)m2tk(m0)
m0 |
moment tensor eigenvalues, sorted decending |
Convert 3 eigen values of a moment tensor to T-k coordinates
list(t, k)
Keehoon Kim<[email protected]> Jonathan M. Lees<[email protected]>
Hudson
tk2uv, hudson.net, hudson.plot
v = c(2,-1,-1) m2tk(v)v = c(2,-1,-1) m2tk(v)
Given vertices of a 3D block, create a glyph structure (faces and normals)
makeblock3D(block1)makeblock3D(block1)
block1 |
matrix of vertices |
glyph structure list
aglyph |
list of faces (x,y,z) |
anorm |
Normals to faces |
Jonathan M. Lees<[email protected]>
ROTZ, ROTY, ROTX, BOXarrows3D, Z3Darrow, TRANmat
block1 = matrix(c(0,0,0, 1,0,0, 1,0.5,0, 0,0.5,0, 0,0,-2, 1,0,-2, 1,0.5,-2, 0,0.5,-2), byrow=TRUE, ncol=3) Bblock1 = makeblock3D(block1)block1 = matrix(c(0,0,0, 1,0,0, 1,0.5,0, 0,0.5,0, 0,0,-2, 1,0,-2, 1,0.5,-2, 0,0.5,-2), byrow=TRUE, ncol=3) Bblock1 = makeblock3D(block1)
Creates but does not plot an Equal-Angle (Schmidt) Stereonet
makenet()makenet()
list of x,y, values for drawing lines
x1 |
x-coordinate start of lines |
y1 |
y-coordinate start of lines |
x2 |
x-coordinate end of lines |
y2 |
y-coordinate end of lines |
Jonathan M. Lees <[email protected]>
Snyder, John P., 1987, Map Projections-a working manual, USGS-Professional Paper, 383p. pages 185-186
net, pnet
MN = makenet() pnet(MN)MN = makenet() pnet(MN)
Plot moment tensors on map
MapNonDouble(Locs, moments, sel = 1, siz = 0.2, col=rgb(1, .75, .75), PLANES = TRUE, add = FALSE, LEG=FALSE)MapNonDouble(Locs, moments, sel = 1, siz = 0.2, col=rgb(1, .75, .75), PLANES = TRUE, add = FALSE, LEG=FALSE)
Locs |
Locations, x,y |
moments |
list of moments: seven elements. See details. |
sel |
integer, index of which to plot |
siz |
size to plot, inches |
col |
color, either a single color, rgb, or a color palette. |
PLANES |
logical, whether to add nodal planes, default=TRUE |
add |
logical, whether to add to plot, default=FALSE |
LEG |
logical, whether to add focal mech legend based on color coding, default=FALSE |
Moment tensors are added to an existing plot. The first element of the list is the integer index of the event. The next six elements are the moments in the following order, c(Mxx, Myy, Mzz, Mzy, Mxz, Mxy) .
If the data is in spherical coordinates, one must switch the
sign of the Mrp and Mtp components, so:
Mrr = Mzz
Mtt = Mxx
Mpp = Myy
Mrt = Mxz
Mrp = -Myz
Mtp = -Mxy
A color palette can be provided for some details of the radiation patterns, e.g. col=rainbow(12). If col is NULL, the colors will be chosen according to focal.color from RFOC, based on rake of first nodal plane.
If col is NULL, then the colors are set by foc.color and it is appropriate to add a legend.
list:
FOC |
matrix, focal mechanism angles (strike, dip rake) |
LAB |
matrix, x-y location for labels |
If events are read in using spherical rather than cartesian
coordinates
need a conversion:
Mrr = Mzz
Mtt = Mxx
Mpp = Myy
Mrt = Mxz
Mrp = -Myz
Mtp = -Mxy
Jonathan M. Lees<[email protected]>
Ekstrom, G.; Nettles, M. & DziewoDski, A. The Global CMT Project 2004-2010: centroid-moment tensors for 13,017 earthquakes Physics of the Earth and Planetary Interiors, 2012.
doNonDouble, ShadowCLVD, angles, nodalLines, PTaxes, focal.color, foc.icolor
## Not run: library(maps) library(GEOmap) ########## load the data data(widdenMoments) ################# to read in the data from a file, ## GG = scan("widdenMoments.txt",sep=" ", ## what=list(ID=0,Event="",Lat=0,Long=0,Depth=0,Mw=0,ML=0,DC=0, ## CLVD=0,ISO=0,VR=0,nsta=0,Mxx=0,Mxy=0,Mxz=0, ## Myy=0,Myz=0,Mzz=0,Mo=0,Ftest=0) ) GG = widdenMoments Locs = list(y=GG$Lat,x=GG$Long) ef = 1e20 moments = cbind(GG$ID, ef*GG$Mxx, ef*GG$Myy, ef*GG$Mzz, ef*GG$Myz, ef*GG$Mxz,ef*GG$Mxy) UTAH = map('state', region = c('utah'), plot=FALSE ) mlon = mean(UTAH$x, na.rm=TRUE) mlat = mean(UTAH$y, na.rm=TRUE) Gutah = maps2GEOmap(UTAH) ############ for mercator projection PROJ = GEOmap::setPROJ(type = 1, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) ############ for UTM projection PROJ = GEOmap::setPROJ(type = 2, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) LIMlat = expandbound(Gutah$POINTS$lat) LIMlon = expandbound(Gutah$POINTS$lon) PLAT = pretty(LIMlat) PLON = pretty(LIMlon) ############### plot the map ######## Utah is a little rectangular dev.new(width=9, height=12) plotGEOmapXY(Gutah, LIM = c(min(PLON), min(PLAT) , max(PLON) , max(PLAT)) , PROJ=PROJ, axes=FALSE, xlab="", ylab="" ) ### add tic marks kbox = GEOmap::GLOB.XY(PLAT,PLON, PROJ) sqrTICXY(kbox , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ######## add focal mechs siz = 0.2 MapNonDouble(Glocs, moments,col=NULL, add=TRUE, LEG=TRUE) up = par("usr") ui = par("pin") ratx = (up[2] - up[1])/ui[1] raty = (up[4] - up[3])/ui[2] usizx = siz * ratx AXY = NoOverlap(Glocs$x,Glocs$y, usizx ) MapNonDouble(AXY, moments,col=NULL, add=TRUE, LEG=TRUE) #### MapNonDouble(NXY, moments,col=NULL, add=TRUE, LEG=TRUE) ## End(Not run)## Not run: library(maps) library(GEOmap) ########## load the data data(widdenMoments) ################# to read in the data from a file, ## GG = scan("widdenMoments.txt",sep=" ", ## what=list(ID=0,Event="",Lat=0,Long=0,Depth=0,Mw=0,ML=0,DC=0, ## CLVD=0,ISO=0,VR=0,nsta=0,Mxx=0,Mxy=0,Mxz=0, ## Myy=0,Myz=0,Mzz=0,Mo=0,Ftest=0) ) GG = widdenMoments Locs = list(y=GG$Lat,x=GG$Long) ef = 1e20 moments = cbind(GG$ID, ef*GG$Mxx, ef*GG$Myy, ef*GG$Mzz, ef*GG$Myz, ef*GG$Mxz,ef*GG$Mxy) UTAH = map('state', region = c('utah'), plot=FALSE ) mlon = mean(UTAH$x, na.rm=TRUE) mlat = mean(UTAH$y, na.rm=TRUE) Gutah = maps2GEOmap(UTAH) ############ for mercator projection PROJ = GEOmap::setPROJ(type = 1, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) ############ for UTM projection PROJ = GEOmap::setPROJ(type = 2, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) LIMlat = expandbound(Gutah$POINTS$lat) LIMlon = expandbound(Gutah$POINTS$lon) PLAT = pretty(LIMlat) PLON = pretty(LIMlon) ############### plot the map ######## Utah is a little rectangular dev.new(width=9, height=12) plotGEOmapXY(Gutah, LIM = c(min(PLON), min(PLAT) , max(PLON) , max(PLAT)) , PROJ=PROJ, axes=FALSE, xlab="", ylab="" ) ### add tic marks kbox = GEOmap::GLOB.XY(PLAT,PLON, PROJ) sqrTICXY(kbox , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ######## add focal mechs siz = 0.2 MapNonDouble(Glocs, moments,col=NULL, add=TRUE, LEG=TRUE) up = par("usr") ui = par("pin") ratx = (up[2] - up[1])/ui[1] raty = (up[4] - up[3])/ui[2] usizx = siz * ratx AXY = NoOverlap(Glocs$x,Glocs$y, usizx ) MapNonDouble(AXY, moments,col=NULL, add=TRUE, LEG=TRUE) #### MapNonDouble(NXY, moments,col=NULL, add=TRUE, LEG=TRUE) ## End(Not run)
takes the pole information from a steroplot and returns the cartesian coordinates
mc2cart(az, dip)mc2cart(az, dip)
az |
degrees, orientation angle, from North |
dip |
degrees, dip of pole |
list of x,y,z values
Jonathan M. Lees <[email protected]>
v1 = mc2cart(65,32) v2 = mc2cart(135,74)v1 = mc2cart(65,32) v2 = mc2cart(135,74)
Convert a normalized moment tensor from the CMT catalog to Strike-Dip-Rake.
mijsdr(mxx, myy, mzz, mxy, mxz, myz)mijsdr(mxx, myy, mzz, mxy, mxz, myz)
mxx |
moment tensor 1,1 |
myy |
moment tensor 2,2 |
mzz |
moment tensor 3,3 |
mxy |
moment tensor 1,2 |
mxz |
moment tensor 1,3 |
myz |
moment tensor 2,3 |
the coordinate system is modified to represent a system centered on the source.
Focal Mechanism list
This code will convert the output of the website, http://www.globalcmt.org/CMTsearch.html when dumped in the psmeca (GMT v>3.3) format.
Jonathan M. Lees<[email protected]>
http://www.globalcmt.org/CMTsearch.html
getCMT
mijsdr(-1.96, 1.07, 0.89, 0.51, 0.08, -0.68)mijsdr(-1.96, 1.07, 0.89, 0.51, 0.08, -0.68)
Calculate the distance between moment tensors based on quaternions.
MomentDist(E1, E2)MomentDist(E1, E2)
E1 |
Moment tensor |
E2 |
Moment tensor |
Moment tensors should be right handed.
angle in degrees
Jonathan M. Lees<[email protected]>
Tape and Tape, 2012
forcerighthand, testrightHAND
Mtens = c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) Mtens = c(5.054, -2.235, -2.819, -0.476, 5.420, 5.594) M2 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 = eigen(M1) ### make sure these are a right handed system, ### ie x1 cross x2 = x3 E2 = eigen(M2) ### make sure these are a right handed system, ### ie x1 cross x2 = x3 testrightHAND(E1$vectors) testrightHAND(E2$vectors) E1$vectors = forcerighthand(E1$vectors) E2$vectors = forcerighthand(E2$vectors) testrightHAND(E1$vectors) testrightHAND(E2$vectors) MomentDist(E1, E2)Mtens = c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) Mtens = c(5.054, -2.235, -2.819, -0.476, 5.420, 5.594) M2 = matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 = eigen(M1) ### make sure these are a right handed system, ### ie x1 cross x2 = x3 E2 = eigen(M2) ### make sure these are a right handed system, ### ie x1 cross x2 = x3 testrightHAND(E1$vectors) testrightHAND(E2$vectors) E1$vectors = forcerighthand(E1$vectors) E2$vectors = forcerighthand(E2$vectors) testrightHAND(E1$vectors) testrightHAND(E2$vectors) MomentDist(E1, E2)
Calculate various parameters associated with the Rake or Slip of an earthquake
MRake(M)MRake(M)
M |
list(uaz, ud, vaz, vd, paz, pd, taz, td) |
This routine takes the four poles U, V, P, T, and returns a MEC structure. (uaz, ud ) = U pole azimuth and dip ( vaz, vd)= V pole azimuth and dip (paz, pd)= P pole azimuth and dip (taz, td)= T pole azimuth and dip
returns a MEC structure
Jonathan M. Lees <[email protected]>
CONVERTSDR, GetRakeSense, GetRake
mc = CONVERTSDR(329, 8, 110 ) MEC = MRake(mc$M)mc = CONVERTSDR(329, 8, 110 ) MEC = MRake(mc$M)
Plot Equal Area Stereo-Net. Lambert azimuthal Equal-Area (Schmidt) from Snyder p. 185-186
net(add = FALSE, col = gray(0.7), border = "black", lwd = 1, LIM = c(-1, -1, +1, +1))net(add = FALSE, col = gray(0.7), border = "black", lwd = 1, LIM = c(-1, -1, +1, +1))
add |
logical, TRUE=add to existing plot |
col |
color of lines |
border |
color of outer rim of stereonet |
lwd |
linewidth of lines |
LIM |
bounding area for a new plot |
Used for graphical side effects
Jonathan M. Lees <[email protected]>
Snyder, John P., 1987, Map Projections-a working manual, USGS-Professional Paper, 383p. pages 185-186
pcirc
net(FALSE, col=rgb(.8,.7,.7) ,border='blue' )net(FALSE, col=rgb(.8,.7,.7) ,border='blue' )
Plots a fault plane and the slip vector. Used for geographic representation of numerous focal spheres.
nipXY(MEC, x = x, y = y, focsiz=1, fcol = gray(0.9), nipcol = "black", cex = 0.4)nipXY(MEC, x = x, y = y, focsiz=1, fcol = gray(0.9), nipcol = "black", cex = 0.4)
MEC |
MEC structure |
x |
coordinate on plot |
y |
coordinate on plot |
focsiz |
size in inches |
fcol |
color for plotting |
nipcol |
color of slip point |
cex |
character expansion for slip point |
Slip vector is the cross product of the poles to the fault plane and auxilliary planes.
LIST
Q |
output of qpoint |
N |
slip vector |
Jonathan M. Lees<[email protected]>
qpoint, CROSSL, lowplane, TOCART
set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1, xlab='km', ylab='km' ) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) nipXY(MEC, x = XY$x[i], y = XY$y[i], focsiz=0.5, fcol = jcol, nipcol = 'black' , cex = 1) }set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1, xlab='km', ylab='km' ) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) nipXY(MEC, x = XY$x[i], y = XY$y[i], focsiz=0.5, fcol = jcol, nipcol = 'black' , cex = 1) }
Add nodal planes to focal mechanism
nodalLines(strike, dip, rake, PLOT=TRUE)nodalLines(strike, dip, rake, PLOT=TRUE)
strike |
numeric, strike of fault |
dip |
numeric, dip of fault |
rake |
numeric, rake of fault |
PLOT |
logical, add lines to plot, default=TRUE |
Lower Hemisphere focal plane.
Side effects
Lower Hemisphere based on FOCangles.
Jonathan M. Lees<[email protected]>
doNonDouble, MapNonDouble, FOCangles
mo <- list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments <- cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)mo <- list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments <- cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)
Illustrate a normal fault using animation
normal.fault(ANG = (45), anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 4, Light = c(45, 45))normal.fault(ANG = (45), anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 4, Light = c(45, 45))
ANG |
Angle of dip |
anim |
animation vector |
KAPPA |
Phong parameter for lighting |
Light |
lighting point |
Program will animate a normal fault for educational purposes. Animation must be stopped by halting execution.
Graphical Side effects
Jonathan M. Lees<[email protected]>
strikeslip.fault, thrust.fault
normal.fault(45, anim=0, KAPPA=4, Light=c(-20, 80)) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) normal.fault(45, anim=anim, KAPPA=4, Light=c(-20, 80)) ## End(Not run)normal.fault(45, anim=0, KAPPA=4, Light=c(-20, 80)) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) normal.fault(45, anim=anim, KAPPA=4, Light=c(-20, 80)) ## End(Not run)
Add a circle to a plot, with cross-hairs
pcirc(gcol = "black", border = "black", ndiv = 36)pcirc(gcol = "black", border = "black", ndiv = 36)
gcol |
color of crosshairs |
border |
border color |
ndiv |
number of divisions for the circle |
no return values, used for side effects
Jonathan M. Lees <[email protected]>
net
net() pcirc(gcol = "green", border = "purple", ndiv = 36)net() pcirc(gcol = "green", border = "purple", ndiv = 36)
rotates a body in 3D and plots projection on existing plot
pglyph3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), anorms = list(), zee = c(0, 0, 1), col = "white", border = "black")pglyph3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), anorms = list(), zee = c(0, 0, 1), col = "white", border = "black")
aglyph |
glyph structure describing the vertices and normal vectors of a 3D body |
M |
rotation matrix 1 |
M2 |
rotation matrix 2 |
anorms |
up vector |
zee |
up vector |
col |
coor of body |
border |
color of border |
Hidden sides are removed and phong shading is introduced to create 3D effect.
The input consists of an object defined by a list structure, list(aglyph, anorm) where aglyph is list of 3D polygons (faces) and anorm are outward normals to these faces.
Used for side effect on plots
For unusual rotations or bizarre bodies, this routine may produce strange looking shapes.
Jonathan M. Lees <[email protected]>
Rogers and Adams, 1990, Mathematical Elements for Computer Graphics, McGraw-Hill, 611p.
Z3Darrow, ROTX, ROTY, ROTZ
### create the 3D object len = .7 basethick=.05 headlip=.02 headlen=.3 #### create a 3D glyph structure aglyph = Z3Darrow(len = len , basethick =basethick , headlen =headlen , headlip=headlip ) #### define the up vector myzee = matrix(c(0,0,1, 1), nrow=1, ncol=4) ##### set rotation angles: gamma =12 beta =39 alpha = 62 ######## set up rotation matrix R3 = ROTZ(gamma) R2 = ROTY(beta) R1 = ROTZ(alpha) ### create rotation matrix M = R1 M2 = R1 plot(c(-1,1), c(-1,1)) pglyph3D(aglyph$aglyph, anorms=aglyph$anorm , M=M, M2=M2, zee=myzee , col=rgb(.7, 0,0) )### create the 3D object len = .7 basethick=.05 headlip=.02 headlen=.3 #### create a 3D glyph structure aglyph = Z3Darrow(len = len , basethick =basethick , headlen =headlen , headlip=headlip ) #### define the up vector myzee = matrix(c(0,0,1, 1), nrow=1, ncol=4) ##### set rotation angles: gamma =12 beta =39 alpha = 62 ######## set up rotation matrix R3 = ROTZ(gamma) R2 = ROTY(beta) R1 = ROTZ(alpha) ### create rotation matrix M = R1 M2 = R1 plot(c(-1,1), c(-1,1)) pglyph3D(aglyph$aglyph, anorms=aglyph$anorm , M=M, M2=M2, zee=myzee , col=rgb(.7, 0,0) )
Create phong shading for faces showing on the 3D block
phong3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), Light = c(45, 45), anorms = list(), zee = c(0, 0, 1), col = "white", border = "black")phong3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), Light = c(45, 45), anorms = list(), zee = c(0, 0, 1), col = "white", border = "black")
aglyph |
3-D body list of faces and normals |
M |
Rotation Matrix |
M2 |
Viewing Matrix |
Light |
light source direction |
anorms |
normals to faces |
zee |
Up vector for Body |
col |
color for faces |
border |
border color for sides |
Uses a standard phong shading model based ont eh dot product of the face normal vector and direction of incoming light.
Graphical Side effect
Jonathan M. Lees<[email protected]>
Watt, Alan. Fundamentals of Three-dimensional Computer Graphics, Addison-Wesley, 1989, 430p.
makeblock3D, BOXarrows3D, PROJ3D, Z3Darrow, pglyph3D
########### create a block and rotation matrix, then color it ANG=(45) DEGRAD = pi/180 y1 = 1.5 y2 = y1 - 1/tan((ANG)*DEGRAD) z1 = 1 x1 = 1 Ablock1 = matrix(c(0,0,0, 1,0,0, 1,y1,0, 0,y1,0, 0,0,-1, 1,0,-1, 1,y2,-1, 0,y2,-1), byrow=TRUE, ncol=3) Nblock1 = makeblock3D(Ablock1) Light=c(45,45) angz = -45 angx = -45 R1 = ROTZ(angz) R2 = ROTX(angx) M = R1 Z2 = PROJ3D(Nblock1$aglyph, M=M, anorms=Nblock1$anorm , zee=c(0,0,1)) RangesX = range(attr(Z2, "RangesX")) RangesY = range(attr(Z2, "RangesY")) plot( RangesX, RangesY, type='n', asp=1, ann=FALSE, axes=FALSE) phong3D(Nblock1$aglyph, M=M, anorms=Nblock1$anorm , Light = Light, zee=c(0,0,1), col=rgb(.7,.5, .5) , border="black")########### create a block and rotation matrix, then color it ANG=(45) DEGRAD = pi/180 y1 = 1.5 y2 = y1 - 1/tan((ANG)*DEGRAD) z1 = 1 x1 = 1 Ablock1 = matrix(c(0,0,0, 1,0,0, 1,y1,0, 0,y1,0, 0,0,-1, 1,0,-1, 1,y2,-1, 0,y2,-1), byrow=TRUE, ncol=3) Nblock1 = makeblock3D(Ablock1) Light=c(45,45) angz = -45 angx = -45 R1 = ROTZ(angz) R2 = ROTX(angx) M = R1 Z2 = PROJ3D(Nblock1$aglyph, M=M, anorms=Nblock1$anorm , zee=c(0,0,1)) RangesX = range(attr(Z2, "RangesX")) RangesY = range(attr(Z2, "RangesY")) plot( RangesX, RangesY, type='n', asp=1, ann=FALSE, axes=FALSE) phong3D(Nblock1$aglyph, M=M, anorms=Nblock1$anorm , Light = Light, zee=c(0,0,1), col=rgb(.7,.5, .5) , border="black")
P and T-axes and Locations of Harvard CMT catalog for the intersection of the Kamchataka and Aleutian arcs
data(PKAM)data(PKAM)
The format is: chr "PKAM"
The data is selected from the CMT catalog. Parameters are extracted from the standard web distribution. Format of the list of data save in PKAM is:
itemPazP-axis azimuth angle itemPdipP-axis dip angle itemTazT-axis azimuth angle itemTdipT-axis dip angle itemhhorizontal point to plot on ternary plot itemvvertical point to plot on ternary plot itemfcolscolors, not used itemLATSLatitude itemLONSLongitude itemIFcolinteger pointer to internal color itemyryear, not used itemJDHMJulian Day, hour, minute, not used itemJDHMSJulian Day, hour, minute, seconds
http://www.globalcmt.org/CMTsearch.html
G. Ekstrom. Rapid earthquake analysis utilizes the internet. Computers in Physics, 8:632-638, 1994.
data(PKAM) ## ###### plot the locations: plot( RPMG::fmod(PKAM$LONS, 360), PKAM$LATS) ###### PlotTernfoc(PKAM$h,PKAM$v,x=0, y=0, siz=1, fcols='black', add=FALSE, LAB=TRUE) ###### change the colors for the plot acols = rainbow(7) fcols = acols[PKAM$IFcol] ###### PlotTernfoc(PKAM$h,PKAM$v,x=0, y=0, siz=1, fcols=fcols, add=FALSE, LAB=TRUE)data(PKAM) ## ###### plot the locations: plot( RPMG::fmod(PKAM$LONS, 360), PKAM$LATS) ###### PlotTernfoc(PKAM$h,PKAM$v,x=0, y=0, siz=1, fcols='black', add=FALSE, LAB=TRUE) ###### change the colors for the plot acols = rainbow(7) fcols = acols[PKAM$IFcol] ###### PlotTernfoc(PKAM$h,PKAM$v,x=0, y=0, siz=1, fcols=fcols, add=FALSE, LAB=TRUE)
Takes a MEC structure and plots all three radiation patterns.
plotfoc(MEC)plotfoc(MEC)
MEC |
MEC list |
Plot makes three figures after calling par(mfrow=c(3,1)).
Graphical Side Effects.
Basic MEC List Structure
| az1 | azimuth angle plane 1, degrees |
| dip1 | dip angle plane 1, degrees |
| az2 | azimuth angle plane 2, degrees |
| dip2 | dip angle plane 2, degrees |
| dir | 0,1 to determine which section of focal sphere is shaded |
| rake1 | rake angle plane 1, degrees |
| dipaz1 | dip azimuth angle plane 1, degrees |
| rake2 | rake angle plane 2, degrees |
| dipaz2 | dip azimuth angle plane 2, degrees |
| P | pole list(az, dip) P-axis |
| T | pole list(az, dip) T-axis |
| U | pole list(az, dip) U-axis |
| V | pole list(az, dip) V-axis |
| F | pole list(az, dip) F-axis |
| G | pole list(az, dip) G-axis |
| sense | 0,1 to determine which section of focal sphere is shaded |
| M | list of focal parameters used in some calculations |
| UP | logical, TRUE=upper hemisphere |
| icol | index to suite of colors for focal mechanism |
| ileg | Kind of fault |
| fcol | color of focal mechanism |
| CNVRG | Character, note on convergence of solution |
| LIM | vector plotting region (x1, y1, x2, y2) |
Jonathan M. Lees<[email protected]>
SDRfoc, Mrake, Pradfoc, radiateSH, radP, radSV, SVradfoc, radiateP, radiateSV, radSH, SHradfoc, imageP, imageSH, imageSV
M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=FALSE) plotfoc(M)M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=FALSE) plotfoc(M)
Plot a long list of focal mechanisms
plotmanyfoc(MEK, PROJ, focsiz = 0.5, foccol = NULL, UP=TRUE, focstyle=1, PMAT = NULL, LEG = FALSE, DOBAR = FALSE)plotmanyfoc(MEK, PROJ, focsiz = 0.5, foccol = NULL, UP=TRUE, focstyle=1, PMAT = NULL, LEG = FALSE, DOBAR = FALSE)
MEK |
List of Focal Mechanisms, see details |
PROJ |
Projection |
focsiz |
focal size, inches |
foccol |
focal color |
UP |
logical, UP=TRUE means plot upper hemisphere (DEFAULT=TRUE) |
focstyle |
integer, 1=beach ball, 2=nipplot, 3=strike-slip, 4=P-T, 5=P, 6=T |
PMAT |
Projection Matrix from persp |
LEG |
logical, TRUE= add focal legend for color codes |
DOBAR |
add strike dip bar at epicenter |
Input MEK list contains
MEKS = list(lon=0, lat=0, str1=0, dip1=0, rake1=0, dep=0, name="", Elat=0, Elon=0)
Graphical Side Effects
Jonathan M. Lees<[email protected]>
Lees, J. M., Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences, Computers & Geosciences, 26, 7, 751-761, 2000.
justfocXY
set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.5)set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.5)
Plot a Focal Mechanism
plotMEC(x, detail = 0, up = FALSE)plotMEC(x, detail = 0, up = FALSE)
x |
Mechanism list |
detail |
level of detail |
up |
logical, Upper or lower hemisphere |
Side Effects
Jonathan M. Lees<[email protected]>
mc = CONVERTSDR(65, 32, -34 ) plotMEC(mc, detail=2, up=FALSE)mc = CONVERTSDR(65, 32, -34 ) plotMEC(mc, detail=2, up=FALSE)
Plot both fault and auxilliary planes
PlotPlanes(MEC, col1 = 1, col2 = 3)PlotPlanes(MEC, col1 = 1, col2 = 3)
MEC |
MEC structure |
col1 |
color for plane 1 |
col2 |
color for plane 2 |
Given MEC structure and focal mechanism plot both planes. This code adds to existing plot, so net() should be called.
Graphical Side Effects
Jonathan M. Lees <[email protected]>
net, lowplane
net() MFOC1 = SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) PlotPlanes(MFOC1, 'green', 'red' )net() MFOC1 = SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) PlotPlanes(MFOC1, 'green', 'red' )
Project PT axes on the sphere and smooth the image. This function requires function kde2d, from the MASS library.
PlotPTsmooth(paz, pdip, x = 0, y = 0, siz = 1, bcol = "white", border ="black", IMAGE = TRUE, CONT = TRUE, cont.col = "black", pal = terrain.colors(100), LABS = FALSE, add = FALSE, NCP=50, NIP=200)PlotPTsmooth(paz, pdip, x = 0, y = 0, siz = 1, bcol = "white", border ="black", IMAGE = TRUE, CONT = TRUE, cont.col = "black", pal = terrain.colors(100), LABS = FALSE, add = FALSE, NCP=50, NIP=200)
paz |
vector of Axis azimuths, degrees |
pdip |
vector of dip angles, degrees |
x |
x-location of plot center in user coordinates |
y |
y-location of plot center in user coordinates |
siz |
siz of plot in user coordinates |
bcol |
color |
border |
border color |
IMAGE |
logical, TRUE=create an image plot |
CONT |
logical, TRUE=add contour lines |
cont.col |
color of contour lines |
pal |
pallete for image plot |
LABS |
text Label for image |
add |
logical, TRUE=add to plot |
NCP |
integer, Number of points to use for calculating smoothed contours, default=50 |
NIP |
integer, Number of points to use for calculating smoothed image, default=200 |
Program requires MASS libary for 2D smoothing routine kde2d.
For calculating contours the kde2d program creates a smoothed 2D image using NCP points per side. For the images, NIP points are used. To reduce the size of plots, or, if the subplots are very small, reduce NIP to a smaller value for faster plotting.
Graphical Side Effect
Points that fall on the opposite hemisphere are reflected through the origin.
Jonathan M. Lees<[email protected]>
kde2d
plot(c(-1,1), c(-1,1), asp=1, type='n') paz = rnorm(100, mean=297, sd=10) pdip = rnorm(100, mean=52, sd=8) PlotPTsmooth(paz, pdip, x=0.5, y=.5, siz=.3, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=FALSE) taz = rnorm(100, mean=138, sd=10) tdip = rnorm(100, mean=12, sd=8) PlotPTsmooth(taz, tdip, x=-.5, y=.4, siz=.3, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=TRUE) ########### put them together plot(c(-1,1), c(-1,1), asp=1, type='n') PlotPTsmooth(paz, pdip, x=0, y=, siz=1, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=FALSE) PlotPTsmooth(taz, tdip, x=0, y=, siz=1, border=NA, bcol='white' , LABS=FALSE, add=TRUE, IMAGE=FALSE, CONT=TRUE)plot(c(-1,1), c(-1,1), asp=1, type='n') paz = rnorm(100, mean=297, sd=10) pdip = rnorm(100, mean=52, sd=8) PlotPTsmooth(paz, pdip, x=0.5, y=.5, siz=.3, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=FALSE) taz = rnorm(100, mean=138, sd=10) tdip = rnorm(100, mean=12, sd=8) PlotPTsmooth(taz, tdip, x=-.5, y=.4, siz=.3, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=TRUE) ########### put them together plot(c(-1,1), c(-1,1), asp=1, type='n') PlotPTsmooth(paz, pdip, x=0, y=, siz=1, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=FALSE) PlotPTsmooth(taz, tdip, x=0, y=, siz=1, border=NA, bcol='white' , LABS=FALSE, add=TRUE, IMAGE=FALSE, CONT=TRUE)
Create and plot a ternary diagram using rake angle to distribute focal mechanisms on a ternary diagram.
PlotTernfoc(h, v, x = 0, y = 0, siz = 1, fcols = "black", LABS = FALSE, add = FALSE)PlotTernfoc(h, v, x = 0, y = 0, siz = 1, fcols = "black", LABS = FALSE, add = FALSE)
h |
x-coordinate on ternary plot |
v |
y-coordinate of ternary plot |
x |
x Location of center of Ternary plot |
y |
y Location of center of Ternary plot |
siz |
size of plot in user coordinates |
fcols |
vector of colors associated with each focal mechanism |
LABS |
logical, TRUE=add labels at vertices of Ternary plot |
add |
logical, add to plot=TRUE |
Used for graphical side effect.
Jonathan M. Lees <[email protected]>
J. M. Lees. Geotouch: Software for three and four dimensional gis in the earth sciences. Computers & Geosciences, 26(7):751–761, 2000
ternfoc.point, Bfocvec
Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) h = vector() v = vector() Fcol = vector() for(i in 1:length(MZ[,3])) { Msdr = CONVERTSDR(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td ) h[i] = V$h v[i] = V$v Fcol[i] = foc.color(foc.icolor(MEC$rake1), pal=1) } PlotTernfoc(h,v,x=0, y=0, siz=1, fcols=Fcol, add=FALSE, LAB=TRUE) MFOC1 = SDRfoc(65,90,1, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol1 = foc.color(foc.icolor(MFOC1$rake1), pal=1) MFOC2 = SDRfoc(135,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol2 = foc.color(foc.icolor(MFOC2$rake1), pal=1) MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol3 = foc.color(foc.icolor(MFOC3$rake1), pal=1) justfocXY( MFOC3, fcol = Fcol3, 1.2, -0.9, focsiz = 0.4 ) justfocXY( MFOC2, fcol = Fcol2, -1.2, -0.9, focsiz = 0.4 ) justfocXY( MFOC1, fcol = Fcol1, 0, 1.414443+.2, focsiz = 0.4 )Z1 = c(159.33,51.6,206,18,78, 161.89,54.5,257,27,133, 170.03,53.57,-44,13,171, 154.99,50.16,-83,19,-40, 151.09,47.15,123,23,-170, 176.31,51.41,-81,22,122, 153.71,46.63,205,28,59, 178.39,51.21,-77,16,126, 178.27,51.1,-86,15,115, 177.95,51.14,-83,25,126, 178.25,51.18,215,16,27 ) MZ = matrix(Z1, ncol=5, byrow=TRUE) h = vector() v = vector() Fcol = vector() for(i in 1:length(MZ[,3])) { Msdr = CONVERTSDR(MZ[i,3], MZ[i,4], MZ[i,5]) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td ) h[i] = V$h v[i] = V$v Fcol[i] = foc.color(foc.icolor(MEC$rake1), pal=1) } PlotTernfoc(h,v,x=0, y=0, siz=1, fcols=Fcol, add=FALSE, LAB=TRUE) MFOC1 = SDRfoc(65,90,1, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol1 = foc.color(foc.icolor(MFOC1$rake1), pal=1) MFOC2 = SDRfoc(135,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol2 = foc.color(foc.icolor(MFOC2$rake1), pal=1) MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Fcol3 = foc.color(foc.icolor(MFOC3$rake1), pal=1) justfocXY( MFOC3, fcol = Fcol3, 1.2, -0.9, focsiz = 0.4 ) justfocXY( MFOC2, fcol = Fcol2, -1.2, -0.9, focsiz = 0.4 ) justfocXY( MFOC1, fcol = Fcol1, 0, 1.414443+.2, focsiz = 0.4 )
Plot an arc of a circle with cross-hairs.
PLTcirc(gcol = "black", border = "black", ndiv = 36, angs = c(-pi, pi), PLOT = TRUE, add = FALSE)PLTcirc(gcol = "black", border = "black", ndiv = 36, angs = c(-pi, pi), PLOT = TRUE, add = FALSE)
gcol |
cross hairs color |
border |
border color |
ndiv |
number of divisions |
angs |
vector from angs[1] to angs[2] in radians |
PLOT |
logical, if TRUE plot |
add |
logical, if TRUE add to existing plot |
list used for plotting:
x |
x coordinates |
y |
y coordinates |
phi |
angles, radians |
Jonathan M. Lees <[email protected]>
PLTcirc(gcol = "purple", border = "black", ndiv = 36, angs = c(-pi, pi), PLOT = TRUE, add = FALSE) PLTcirc(gcol = NULL, border = "green" , ndiv = 36, angs = c(-pi/4, pi/4), PLOT = TRUE, add = TRUE)PLTcirc(gcol = "purple", border = "black", ndiv = 36, angs = c(-pi, pi), PLOT = TRUE, add = FALSE) PLTcirc(gcol = NULL, border = "green" , ndiv = 36, angs = c(-pi/4, pi/4), PLOT = TRUE, add = TRUE)
Plots stereonet created by makenet
pnet(MN, add = FALSE, col = gray(0.7), border = "black", lwd = 1)pnet(MN, add = FALSE, col = gray(0.7), border = "black", lwd = 1)
MN |
Net strucutre created by makenet |
add |
TRUE= add to existing plot |
col |
color of lines |
border |
color for outside border |
lwd |
line width |
Used Graphical Side Effects.
Jonathan M. Lees <[email protected]>
Snyder, John P., 1987, Map Projections-a working manual, USGS-Professional Paper, 383p. pages 185-186
net, pnet
MN = makenet() pnet(MN)MN = makenet() pnet(MN)
Calculate the projection of the focal mechanism polygon
polyfoc(strike1, dip1, strike2, dip2, PLOT = FALSE, UP = TRUE)polyfoc(strike1, dip1, strike2, dip2, PLOT = FALSE, UP = TRUE)
strike1 |
strike of plane 1, degrees |
dip1 |
dip of plane 1, degrees |
strike2 |
strike of plane 1, degrees |
dip2 |
dip of plane 2, degrees |
PLOT |
logical, TRUE = add to plot |
UP |
upper/lower hemisphere |
List of coordinates of polygon
Px |
x-coordinates of polygon |
Py |
y-coordinates of polygon |
Jonathan M. Lees <[email protected]>
faultplane
MEC = SDRfoc(13,59,125, PLOT=FALSE) net() ply = polyfoc(MEC$az1, MEC$dip1, MEC$az2, MEC$dip2, PLOT = TRUE, UP = TRUE)MEC = SDRfoc(13,59,125, PLOT=FALSE) net() ply = polyfoc(MEC$az1, MEC$dip1, MEC$az2, MEC$dip2, PLOT = TRUE, UP = TRUE)
Plot P-wave radiation with information from the pickfile and waveform data
Pradfoc(A, MEC, GU, pscale, col)Pradfoc(A, MEC, GU, pscale, col)
A |
Pickfile structure |
MEC |
MEC structure |
GU |
Waveform Event Structure |
pscale |
logical (not used) |
col |
color palette |
Image plot of the P radiation pattern
Graphical Side effects
Jonathan M. Lees<[email protected]>
imageP
MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Pradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) Pradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )
Takes a vector to a pole and reflects it to the lower hemisphere
Preflect(az, dip)Preflect(az, dip)
az |
azimuth angle, degrees |
dip |
dip in degrees |
list
az |
azimuth angle, degrees |
dip |
dip in degrees |
...
Jonathan M. Lees <[email protected]>
REFLECT
z = Preflect(65, -23) z = Preflect(265, -23)z = Preflect(65, -23) z = Preflect(265, -23)
Prepare Focals for plotting. Program cycles through data and prepares a relevant data for further plotting and analysis.
prepFOCS(CMTSOL)prepFOCS(CMTSOL)
CMTSOL |
see getCMT for the format for the input here. |
Used internally in spherefocgeo and ternfocgeo.
List:
Paz |
P-axis azimuth |
Pdip |
P-axis dip |
Taz |
T-axis azimuth |
Tdip |
T-axis dip |
h |
horizontal distance on ternary plot |
v |
vertical distance on ternary plot |
fcols |
focal color |
LATS |
latitudes |
LONS |
longitudes |
IFcol |
index of color |
yr |
year |
JDHM |
character identification |
JDHMS |
character identification |
Jonathan M. Lees<[email protected]>
getCMT, spherefocgeo, ternfocgeo
Print focal mechanism
printMEC(x, digits = max(3, getOption("digits") - 3), ...)printMEC(x, digits = max(3, getOption("digits") - 3), ...)
x |
Mechanism list |
digits |
digits for numeric information |
... |
standard printing parameters |
Side Effects
Jonathan M. Lees<[email protected]>
mc = CONVERTSDR(65, 32, -34 ) printMEC(mc)mc = CONVERTSDR(65, 32, -34 ) printMEC(mc)
Project a 3D body after rotation and translation
PROJ3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), anorms = list(), zee = c(0, 0, 1))PROJ3D(aglyph, M = diag(1, nrow = 4), M2 = diag(1, nrow = 4), anorms = list(), zee = c(0, 0, 1))
aglyph |
glyph structure |
M |
rotation matrix |
M2 |
rotation matrix |
anorms |
normals to structure |
zee |
Up direction of body |
This function takes a 3D body, rotates it and projects it for plotting. An example glyph is found in Z3Darrow.
Glyph structure
x, y, z
|
coordinates of rotated body faces |
xp |
rotated normal vectors |
zd |
depth mean value of each face |
Jonathan M. Lees<[email protected]>
makeblock3D, ROTZ, ROTY, ROTX, BOXarrows3D, Z3Darrow, TRANmat
block1 = matrix(c(0,0,0, 1,0,0, 1,0.5,0, 0,0.5,0, 0,0,-2, 1,0,-2, 1,0.5,-2, 0,0.5,-2), byrow=TRUE, ncol=3) Bblock1 = makeblock3D(block1) R3 = ROTX(-40) R2 = ROTY(0) R1 = ROTZ(20) T = TRANmat(.1, 0, 0 ) M = R1 %*% R2 %*% R3 %*% T T2 = TRANmat(1, 0.5, 0 ) MT = T2 %*% R1 %*% R2 %*% R3 %*% T Z1 = PROJ3D(Bblock1$aglyph, M=MT, anorms=Bblock1$anorm , zee=c(0,0,1))block1 = matrix(c(0,0,0, 1,0,0, 1,0.5,0, 0,0.5,0, 0,0,-2, 1,0,-2, 1,0.5,-2, 0,0.5,-2), byrow=TRUE, ncol=3) Bblock1 = makeblock3D(block1) R3 = ROTX(-40) R2 = ROTY(0) R1 = ROTZ(20) T = TRANmat(.1, 0, 0 ) M = R1 %*% R2 %*% R3 %*% T T2 = TRANmat(1, 0.5, 0 ) MT = T2 %*% R1 %*% R2 %*% R3 %*% T Z1 = PROJ3D(Bblock1$aglyph, M=MT, anorms=Bblock1$anorm , zee=c(0,0,1))
Plot P-T axis on CLVD
PTaxes(strike, dip, rake)PTaxes(strike, dip, rake)
strike |
strike |
dip |
dip |
rake |
rake |
Lower Hemisphere. Add PT axes on a moment tensor plot
Side effects
Jonathan M. Lees<[email protected]>
doNonDouble, MapNonDouble
mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)mo = list(n=1, m1=1.035675e+017, m2=-1.985852e+016, m3=-6.198052e+014, m4=1.177936e+017, m5=-7.600627e+016, m6=-3.461405e+017) moments = cbind(mo$n, mo$m1, mo$m2, mo$m3, mo$m4, mo$m5, mo$m6) doNonDouble(moments)
given a focal mechanism, add P-T lines to a plot
PTXY2(x = x, y = y, MEC, focsiz, pt = 0, ...)PTXY2(x = x, y = y, MEC, focsiz, pt = 0, ...)
x |
x-location on plot |
y |
y-location on plot |
MEC |
Focal Mechanism list from SDRFOC |
focsiz |
size of mechanism, inches |
pt |
pt = 0(plot both), 1=only P axes, 2=only T axes, default=0 |
... |
graphical parameters |
This is a summary plot to be used instead of Beach Balls.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
Lees, J. M., Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences, Computers & Geosciences, 26, 7, 751-761, 2000.
nipXY, justfocXY
### HAiti Earthquake Jan, 2010 MEC <- SDRfoc(71, 64, 25 , u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) plot(c(0, 1), c(0,1), type='n', asp=1) u <- par("usr") justfocXY(MEC, x=.5, y= .5, focsiz=0.5, fcol ='brown' , fcolback = "white", xpd = TRUE) PTXY2(1.0, .5 , MEC ,0.5, col="purple", lwd=3 ) nipXY(MEC, x = 0.25, y = .5, focsiz=0.5, fcol ='purple', nipcol = "black", cex = 0.4) ##### or set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) PTXY2(XY$x[i], XY$y[i] , MEC ,focsiz=0.5, col=jcol, lwd=3) }### HAiti Earthquake Jan, 2010 MEC <- SDRfoc(71, 64, 25 , u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) plot(c(0, 1), c(0,1), type='n', asp=1) u <- par("usr") justfocXY(MEC, x=.5, y= .5, focsiz=0.5, fcol ='brown' , fcolback = "white", xpd = TRUE) PTXY2(1.0, .5 , MEC ,0.5, col="purple", lwd=3 ) nipXY(MEC, x = 0.25, y = .5, focsiz=0.5, fcol ='purple', nipcol = "black", cex = 0.4) ##### or set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) PTXY2(XY$x[i], XY$y[i] , MEC ,focsiz=0.5, col=jcol, lwd=3) }
Plot a set of (azimuths, takeoff) angles on a stereonet.
qpoint(az, iang, col = 2, pch = 5, lab = "", POS = 4, UP = FALSE, PLOT = FALSE, cex = 1)qpoint(az, iang, col = 2, pch = 5, lab = "", POS = 4, UP = FALSE, PLOT = FALSE, cex = 1)
az |
vector of azimuths, degrees |
iang |
vector of incident angles, degrees |
col |
color |
pch |
plotting character |
lab |
text labels |
POS |
position for labels |
UP |
logical, TRUE=upper |
PLOT |
logical, add to existing plot |
cex |
character expansion of labels |
The iang argument represents the takeoff angle, and is measured from the nadir (z-axis pointing down).
List
x |
coordinate on plot |
y |
coordinate on plot |
Jonathan M. Lees<[email protected]>
FixDip, focpoint
d = runif(10, 0, 90) a = runif(10, 0,360) net() qpoint(a, d)d = runif(10, 0, 90) a = runif(10, 0,360) net() qpoint(a, d)
Plots focal mechanism and makes radiation plot with mark up
radiateP(MEC, SCALE = FALSE, col = col, TIT = FALSE)radiateP(MEC, SCALE = FALSE, col = col, TIT = FALSE)
MEC |
focal mechanism structure |
SCALE |
logical, TRUE=add scale |
col |
color palette |
TIT |
title for plot |
Used for side graphical effect
Jonathan M. Lees <[email protected]>
radP, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateP(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateP(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)
Plots focal mechanism and makes radiation plot with mark up
radiateSH(MEC, SCALE = FALSE, col = col, TIT = FALSE)radiateSH(MEC, SCALE = FALSE, col = col, TIT = FALSE)
MEC |
focal mechanism structure |
SCALE |
logical, TRUE=add scale |
col |
color palette |
TIT |
title for plot |
Used for side graphical effect
Jonathan M. Lees <[email protected]>
radSH, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateSH(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateSH(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)
Plots focal mechanism and makes radiation plot with mark up
radiateSV(MEC, SCALE = FALSE, col = col, TIT = FALSE)radiateSV(MEC, SCALE = FALSE, col = col, TIT = FALSE)
MEC |
focal mechanism structure |
SCALE |
logical, TRUE=add scale |
col |
color palette |
TIT |
title for plot |
Used for side graphical effect
Jonathan M. Lees <[email protected]>
radSV, SDRfoc
MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateSV(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)MEC =SDRfoc(65,25,13, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) radiateSV(MEC, SCALE = FALSE, col = rainbow(100) , TIT = FALSE)
calculate the radiation patterns for P waves
radP(del, phiS, lam, ichi, phi)radP(del, phiS, lam, ichi, phi)
del |
degrees, angle |
phiS |
degrees,angle |
lam |
degrees, angle |
ichi |
degrees, take off angle |
phi |
degrees, take off azimuth |
Given a focal mechanism strike-dip-rake and a given incident angle (take-off angle) and azimuth, return the P amplitude
Amplitude of the P wave
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radP, radSV, imageP
phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radP(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radP(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)
calculate the radiation patterns for SH waves
radSH(del, phiS, lam, ichi, phi)radSH(del, phiS, lam, ichi, phi)
del |
degrees, angle |
phiS |
degrees,angle |
lam |
degrees, angle |
ichi |
degrees, take off angle |
phi |
degrees, take off azimuth |
Given a focal mechanism strike-dip-rake and a given incident angle (take-off angle) and azimuth, return the SH amplitude
Amplitude of the SH wave
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radP, radSV, imageSH
phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radSH(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radSH(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)
calculate the radiation patterns for SV waves
radSV(del, phiS, lam, ichi, phi)radSV(del, phiS, lam, ichi, phi)
del |
degrees, angle |
phiS |
degrees,angle |
lam |
degrees, angle |
ichi |
degrees, take off angle |
phi |
degrees, take off azimuth |
Given a focal mechanism strike-dip-rake and a given incident angle (take-off angle) and azimuth, return the SV amplitude
Amplitude of the SV wave
Jonathan M. Lees <[email protected]>
K.~Aki and P.~G. Richards.Quantitative seismology. University Science Books, Sausalito, Calif., 2nd edition, 2002.
radP, radSH, imageSV
phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radSV(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)phiS=65 del=25 lam=13 x = seq(-1, 1, 0.01) y = x X = matrix(rep(x, length(y)), nrow= length(x)) Y = t(X) RAD2DEG = 180/pi p = RAD2DEG*(pi/2 -atan2(Y, X)) p[p<0] = p[p<0] + 360 R = sqrt(X^2+Y^2) R[R>1] = NaN dip =RAD2DEG*2*asin(R/sqrt(2)) ### Calculate the radiation pattern G = radSV(del, phiS, lam, dip, p) ### plot values image(x,y,G, asp=1)
Focal Legend based on rake
rakelegend(corn="topright", pal=1)rakelegend(corn="topright", pal=1)
corn |
position of legend, default="topright" |
pal |
palette number, default=1 |
Colors are based on earlier publication of Geotouch program.
For pal = 1, colors are , DarkSeaGreen, cyan1, SkyBlue1, RoyalBlue, GreenYellow, orange, red.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
Lees, J. M., (1999) Geotouch: Software for Three and Four-Dimensional GIS in the Earth Sciences, Computers and Geosciences, 26(7) 751-761.
foc.color,focleg
plot(c(0,1), c(0,1), type='n') rakelegend(corn="topleft", pal=1)plot(c(0,1), c(0,1), type='n') rakelegend(corn="topleft", pal=1)
Read and plot a CMT solution copied from the Harvard CMT website.
readCMT(filename, PLOT=TRUE)readCMT(filename, PLOT=TRUE)
filename |
character, file name |
PLOT |
Logical, TRUE=plot mechanisms sequentially |
Uses the standard output format.
List of mechanisms and graphical Side effects. Each element in the list consists of a list including: FIRST,yr,mo,dom,hr,mi,sec,name,tshift,half,lat,lon,z,Mrr,Mtt,Mpp,Mrt,Mrp,Mtp. The FIRST element is simply a duplicate of the PDE solution card.
Other formats are available.
Jonathan M. Lees<[email protected]>
Ekstrom, G.; Nettles, M. and DziewoDski, A. The Global CMT Project 2004-2010: centroid-moment tensors for 13,017 earthquakes Physics of the Earth and Planetary Interiors, 2012.
doNonDouble, MapNonDouble
## Not run: Hcmt = readCMT("CMT_FULL_FORMAT.txt") ######## or, Hcmt = readCMT("CMT_FULL_FORMAT.txt", PLOT=FALSE) moments = matrix(ncol=7, nrow=length(Hcmt)) Locs = list(y=vector(length=length(Hcmt)) ,x=vector(length=length(Hcmt))) for(i in 1:length(Hcmt)) { P1 = Hcmt[[i]] ######### Note the change of sign for cartesian coordinates moments[i,] = cbind(i, P1$Mtt, P1$Mpp, P1$Mrr, -P1$Mrp, P1$Mrt ,-P1$Mtp) Locs$y[i] = P1$lat Locs$x[i] = P1$lon } mlon = mean(Locs$x, na.rm=TRUE) mlat = mean(Locs$y, na.rm=TRUE) PROJ = GEOmap::setPROJ(type = 1, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) LIMlat = expandbound(Locs$y) LIMlon = expandbound(Locs$x) PLAT = pretty(LIMlat) PLON = pretty(LIMlon) data(worldmap) par(xpd=FALSE) plotGEOmapXY(worldmap, LIM = c(LIMlon[1],LIMlat[1] ,LIMlon[2],LIMlat[2]) , PROJ=PROJ, axes=FALSE, xlab="", ylab="" ) ### add tic marks kbox = GEOmap::GLOB.XY(PLAT,PLON, PROJ) sqrTICXY(kbox , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ######## add focal mechs MapNonDouble(Glocs, moments, col=NULL, add=TRUE) ## End(Not run)## Not run: Hcmt = readCMT("CMT_FULL_FORMAT.txt") ######## or, Hcmt = readCMT("CMT_FULL_FORMAT.txt", PLOT=FALSE) moments = matrix(ncol=7, nrow=length(Hcmt)) Locs = list(y=vector(length=length(Hcmt)) ,x=vector(length=length(Hcmt))) for(i in 1:length(Hcmt)) { P1 = Hcmt[[i]] ######### Note the change of sign for cartesian coordinates moments[i,] = cbind(i, P1$Mtt, P1$Mpp, P1$Mrr, -P1$Mrp, P1$Mrt ,-P1$Mtp) Locs$y[i] = P1$lat Locs$x[i] = P1$lon } mlon = mean(Locs$x, na.rm=TRUE) mlat = mean(Locs$y, na.rm=TRUE) PROJ = GEOmap::setPROJ(type = 1, LAT0 = mlat , LON0 = mlon) Glocs = GEOmap::GLOB.XY(Locs$y, Locs$x, PROJ ) LIMlat = expandbound(Locs$y) LIMlon = expandbound(Locs$x) PLAT = pretty(LIMlat) PLON = pretty(LIMlon) data(worldmap) par(xpd=FALSE) plotGEOmapXY(worldmap, LIM = c(LIMlon[1],LIMlat[1] ,LIMlon[2],LIMlat[2]) , PROJ=PROJ, axes=FALSE, xlab="", ylab="" ) ### add tic marks kbox = GEOmap::GLOB.XY(PLAT,PLON, PROJ) sqrTICXY(kbox , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ######## add focal mechs MapNonDouble(Glocs, moments, col=NULL, add=TRUE) ## End(Not run)
Given a set of (x,y) points, partition the field into rectangles each containing a minimum number of points
RectDense(INx, INy, icut = 1, u = par("usr"), ndivs = 10)RectDense(INx, INy, icut = 1, u = par("usr"), ndivs = 10)
INx |
x-coordinates |
INy |
y-coordinates |
icut |
cut off for number of points |
u |
user coordinates |
ndivs |
number of divisions in x-coordinate |
Based on the user coordinates as returned from par('usr'). Each rectangular region is tested for the number of points that fall within icut or greater.
List:
icorns |
matrix of corners that passed test |
ilens |
vector,number of points in each icorns box |
ipass |
vector, index of the corners that passed icut |
corners |
matrix of all corners |
lens |
vector,number of points for each box |
Jonathan M. Lees<[email protected]>
x = rnorm(100) y = rnorm(100) plot(x,y) u = par('usr') RI = RectDense(x, y, icut=3, u=u, ndivs=10) rect(RI$icorns[,1],RI$icorns[,2],RI$icorns[,3],RI$icorns[,4], col=NA, border='blue')x = rnorm(100) y = rnorm(100) plot(x,y) u = par('usr') RI = RectDense(x, y, icut=3, u=u, ndivs=10) rect(RI$icorns[,1],RI$icorns[,2],RI$icorns[,3],RI$icorns[,4], col=NA, border='blue')
Reflect pole to lower hemisphere
REFLECT(A)REFLECT(A)
A |
structure of azimuth and Dips in degrees |
list of:cartesian coordinates of reflected pole
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
az |
azimuth, degrees |
dip |
dip, degrees |
Jonathan M. Lees <[email protected]>
Preflect
A = list(az=231, dip = -65) REFLECT(A)A = list(az=231, dip = -65) REFLECT(A)
Rotate mechanism to vertical plan at specified angle
rotateFoc(MEX, phi)rotateFoc(MEX, phi)
MEX |
Focal Mechanism list |
phi |
angle in degrees |
Assumed vertical plane, outer hemisphere
Focal Mechanism
Jonathan M. Lees<[email protected]>
plotfoc, SDRfoc,Beachfoc, TEACHFOC, plotmanyfoc, getUWfocs
a1 = SDRfoc(90, 90, 90, u = TRUE , PLOT = TRUE) par(mfrow=c(2,2)) SDRfoc(a1$az1, a1$dip1, a1$rake1, u = TRUE, PLOT = TRUE) ra1 = rotateFoc(a1, -90) SDRfoc(ra1$az1, ra1$dip1, ra1$rake1, u = TRUE , PLOT = TRUE) ra1 = rotateFoc(a1, 0) SDRfoc(a1$az1, a1$dip1, a1$rake1, u = TRUE, PLOT = TRUE) SDRfoc(ra1$az1, ra1$dip1, ra1$rake1, u = TRUE , PLOT = TRUE)a1 = SDRfoc(90, 90, 90, u = TRUE , PLOT = TRUE) par(mfrow=c(2,2)) SDRfoc(a1$az1, a1$dip1, a1$rake1, u = TRUE, PLOT = TRUE) ra1 = rotateFoc(a1, -90) SDRfoc(ra1$az1, ra1$dip1, ra1$rake1, u = TRUE , PLOT = TRUE) ra1 = rotateFoc(a1, 0) SDRfoc(a1$az1, a1$dip1, a1$rake1, u = TRUE, PLOT = TRUE) SDRfoc(ra1$az1, ra1$dip1, ra1$rake1, u = TRUE , PLOT = TRUE)
Rotate Focal Mechanism into the vertical plane by a certain number of degrees
Rotfocphi(phi, urot, udip, vrot, vdip, az1, d1, az2, d2, prot, pdip, trot, tdip)Rotfocphi(phi, urot, udip, vrot, vdip, az1, d1, az2, d2, prot, pdip, trot, tdip)
phi |
degrees in plane to rotate |
urot |
U-vector azimuth |
udip |
U-vector dip |
vrot |
V-vector azimuth |
vdip |
V-vector dip |
az1 |
First plane - azimuth |
d1 |
First plane - dip |
az2 |
Second plane - azimuth |
d2 |
Second plane - dip |
prot |
P-axis azimuth |
pdip |
P-axis dip |
trot |
T-axis azimuth |
tdip |
T-axis dip |
Rotate the focal mech by phi degrees
list:
Jonathan M. Lees<[email protected]>
xsecmanyfoc, rotateFoc
Rotate T-P axes
RotTP(rotmat, strk1, dip1)RotTP(rotmat, strk1, dip1)
rotmat |
rotation matrix, 3 by 3 |
strk1 |
strike angle |
dip1 |
dip angle |
These are used as functions auxiallry to rotateFoc.
list:
strk |
strike angle |
dip |
dip angle |
Jonathan M. Lees<[email protected]>
Rotfocphi, TP2XYZ
phi = 18 MAT = JMAT(phi) RotTP(MAT, 30, 40)phi = 18 MAT = JMAT(phi) RotTP(MAT, 30, 40)
Matrix rotation about the X-axis
ROTX(deg)ROTX(deg)
deg |
Angle in degrees |
A 4 by 4 matrix for rotation and translation for 3-D transformation
Jonathan M. Lees <[email protected]>
Rogers and Adams, 1990, Mathematical Elements for Computer Graphics, McGraw-Hill, 611p.
ROTY, ROTZ
v = c(1,4,5) A = ROTX(23) vp = c(v, 1)v = c(1,4,5) A = ROTX(23) vp = c(v, 1)
3x3 Rotation about the x axis
rotx3(deg)rotx3(deg)
deg |
angle, degrees |
returns a 3 by 3 rotation matrix
matrix, 3 by 3
Jonathan M. Lees <[email protected]>
roty3, rotz3, ROTX, ROTZ, ROTY
a = 45 rotx3(a)a = 45 rotx3(a)
Matrix rotation about the Y-axis
ROTY(deg)ROTY(deg)
deg |
Angle in degrees |
A 4 by 4 matrix for rotation and translation for 3-D transformation
Jonathan M. Lees <[email protected]>
Rogers and Adams, 1990, Mathematical Elements for Computer Graphics, McGraw-Hill, 611p.
ROTX, ROTZ
v = c(1,4,5) A = ROTY(23) vp = c(v, 1)v = c(1,4,5) A = ROTY(23) vp = c(v, 1)
3x3 Rotation about the y axis
roty3(deg)roty3(deg)
deg |
angle, degrees |
returns a 3 by 3 rotation matrix
matrix, 3 by 3
Jonathan M. Lees <[email protected]>
rotz3, rotx3, ROTX, ROTZ, ROTY
a = 45 roty3(a)a = 45 roty3(a)
Matrix rotation about the Z-axis
ROTZ(deg)ROTZ(deg)
deg |
Angle in degrees |
A 4 by 4 matrix for rotation and translation for 3-D transformation
Jonathan M. Lees <[email protected]>
Rogers and Adams, 1990, Mathematical Elements for Computer Graphics, McGraw-Hill, 611p.
ROTX, ROTY
v = c(1,4,5) A = ROTZ(23) vp = c(v, 1)v = c(1,4,5) A = ROTZ(23) vp = c(v, 1)
3x3 Rotation about the z axis
rotz3(deg)rotz3(deg)
deg |
angle, degrees |
returns a 3 by 3 rotation matrix
matrix, 3 by 3
Jonathan M. Lees <[email protected]>
roty3, rotx3, ROTX, ROTZ, ROTY
a = 45 rotz3(a)a = 45 rotz3(a)
Given Strike-Dip-Rake plot a focal mechanism
SDRfoc(s, d, r, u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT = TRUE)SDRfoc(s, d, r, u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT = TRUE)
s |
strike, degrees |
d |
dip, degrees |
r |
rake, degrees |
u |
logical, TRUE=upper hemisphere |
ALIM |
bounding box on plot |
PLOT |
logical, TRUE=add to plot |
The ALIM vector allows one to zoom into portions of the focal mechanism for details when points are tightly clustered.
MEC structure
Basic MEC List Structure
| az1 | azimuth angle plane 1, degrees |
| dip1 | dip angle plane 1, degrees |
| az2 | azimuth angle plane 2, degrees |
| dip2 | dip angle plane 2, degrees |
| dir | 0,1 to determine which section of focal sphere is shaded |
| rake1 | rake angle plane 1, degrees |
| dipaz1 | dip azimuth angle plane 1, degrees |
| rake2 | rake angle plane 2, degrees |
| dipaz2 | dip azimuth angle plane 2, degrees |
| P | pole list(az, dip) P-axis |
| T | pole list(az, dip) T-axis |
| U | pole list(az, dip) U-axis |
| V | pole list(az, dip) V-axis |
| F | pole list(az, dip) F-axis |
| G | pole list(az, dip) G-axis |
| sense | 0,1 to determine which section of focal sphere is shaded |
| M | list of focal parameters used in some calculations |
| UP | logical, TRUE=upper hemisphere |
| icol | index to suite of colors for focal mechanism |
| ileg | Kind of fault |
| fcol | color of focal mechanism |
| CNVRG | Character, note on convergence of solution |
| LIM | vector plotting region (x1, y1, x2, y2) |
Jonathan M. Lees <[email protected]>
CONVERTSDR
M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=TRUE)M = SDRfoc(-25, 34, 16,u = FALSE, ALIM = c(-1, -1, +1, +1), PLOT=TRUE)
Plot non-double couple part of the focal mechanism provided in the moment tensor.
ShadowCLVD(m, PLOT = TRUE, col=rgb(1, .75, .75))ShadowCLVD(m, PLOT = TRUE, col=rgb(1, .75, .75))
m |
moment tensor |
PLOT |
logical, TRUE means plot |
col |
color, either a single color, rgb, or a color palette |
This code is meant to be used with doNonDouble or MapNonDouble functions for plotting the non-double couple components of the moment tensor. A color palette can be provided for some details of the radiation patterns, e.g. col=rainbow(12).
Side effects and image list
Lower Hemisphere.
Jonathan M. Lees<[email protected]>
doNonDouble, MapNonDouble
############ moment tensor from Harvard CMT catalog sponent = 26 ef = 1*10^(sponent) Mrr = 2.375*ef Mtt = -2.777*ef Mpp = 0.403*ef Mrt = 2.800*ef Mrp = 1.190*ef Mtp = -0.539*ef ############ convert to cartesian coordinates Mzz=Mrr Mxx= Mtt Myy= Mpp Mxz= Mrt Myz= -Mrp Mxy= -Mtp m=matrix( c(Mxx,Mxy,Mxz, Mxy,Myy,Myz, Mxz,Myz,Mzz), ncol=3, byrow=TRUE) Fi=seq(from=0, by=0.1, to=361) ### dev.new() plot(cos(Fi*pi/180.0),sin(Fi*pi/180.0),type='l', asp=1 , ann=FALSE, axes=FALSE) ShadowCLVD(m, col='red')############ moment tensor from Harvard CMT catalog sponent = 26 ef = 1*10^(sponent) Mrr = 2.375*ef Mtt = -2.777*ef Mpp = 0.403*ef Mrt = 2.800*ef Mrp = 1.190*ef Mtp = -0.539*ef ############ convert to cartesian coordinates Mzz=Mrr Mxx= Mtt Myy= Mpp Mxz= Mrt Myz= -Mrp Mxy= -Mtp m=matrix( c(Mxx,Mxy,Mxz, Mxy,Myy,Myz, Mxz,Myz,Mzz), ncol=3, byrow=TRUE) Fi=seq(from=0, by=0.1, to=361) ### dev.new() plot(cos(Fi*pi/180.0),sin(Fi*pi/180.0),type='l', asp=1 , ann=FALSE, axes=FALSE) ShadowCLVD(m, col='red')
Plot SH-wave radiation with information from the pickfile and waveform data
SHradfoc(A, MEC, GU, pscale, col)SHradfoc(A, MEC, GU, pscale, col)
A |
Pickfile structure |
MEC |
MEC structure |
GU |
Waveform Event Structure |
pscale |
logical (not used) |
col |
color palette |
Image plot of the SH radiation pattern
Graphical Side effects
Jonathan M. Lees<[email protected]>
imageSH
MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) SHradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) SHradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )
Given a vector of EigenValues, extract the source type.
SourceType(v)SourceType(v)
v |
vector of decreasing eigenvalues |
plotting for -30 to 30 degree quadrant.
list:
phi |
latitude angle in degrees |
lam |
longitude angle in degrees |
Jonathan M. Lees<[email protected]>
Tape, W.,and C.Tape(2012), A geometric comparison of source-type plots for moment tensors, Geophys. J. Int., 190, 499-510.
HAMMERprojXY, TapeBase, TapePlot
SourceType(c(1,-1,1) ) T1 = TapeBase() m1 = list(Mxx=1.543, Mxy=0.786, Myy=0.336, Mxz=-2.441, Myz=0.353, Mzz=0.961) i = 1 M1=matrix( c(m1$Mxx[i],m1$Mxy[i],m1$Mxz[i], m1$Mxy[i],m1$Myy[i],m1$Myz[i], m1$Mxz[i],m1$Myz[i],m1$Mzz[i]), ncol=3, byrow=TRUE) E1 = eigen(M1) h = SourceType( sort(E1$values, decreasing=TRUE) ) h$dip = 90-h$phi ## cat(paste(h$dip, h$lam, sep=" "), sep="\n") h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) TapePlot(T1) points(h1$x, h1$y, pch=21, bg="red" )SourceType(c(1,-1,1) ) T1 = TapeBase() m1 = list(Mxx=1.543, Mxy=0.786, Myy=0.336, Mxz=-2.441, Myz=0.353, Mzz=0.961) i = 1 M1=matrix( c(m1$Mxx[i],m1$Mxy[i],m1$Mxz[i], m1$Mxy[i],m1$Myy[i],m1$Myz[i], m1$Mxz[i],m1$Myz[i],m1$Mzz[i]), ncol=3, byrow=TRUE) E1 = eigen(M1) h = SourceType( sort(E1$values, decreasing=TRUE) ) h$dip = 90-h$phi ## cat(paste(h$dip, h$lam, sep=" "), sep="\n") h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) TapePlot(T1) points(h1$x, h1$y, pch=21, bg="red" )
Spherical Projections of PT axes distributed geographically.
spherefocgeo(CMTSOL, PROJ = NULL, icut = 5, ndivs = 10, bbox=c(0,1, 0, 1), PLOT = TRUE, add = FALSE, RECT = FALSE, pal = terrain.colors(100))spherefocgeo(CMTSOL, PROJ = NULL, icut = 5, ndivs = 10, bbox=c(0,1, 0, 1), PLOT = TRUE, add = FALSE, RECT = FALSE, pal = terrain.colors(100))
CMTSOL |
see output of getCMT for list input |
PROJ |
Map projection |
icut |
cut off for number of points in box, default=5 |
ndivs |
divisions of map area, default=10 |
bbox |
bounding box for dividing the area, given as minX, maxX, minY, maxY; default=usr coordinates from par() |
PLOT |
logical, default=TRUE |
add |
logical, add to existing plot |
RECT |
logical, TRUE=plot rectangles |
pal |
palette fo rimages in each box |
Program divides the area into blocks, tests each one for minimum number per block and projects the P and T axes onto an equal area stereonet.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
PlotPTsmooth, ternfocgeo, prepFOCS, RectDense
N = 100 LATS = c(7.593004, 25.926045) LONS = c(268.1563 , 305) lon=rnorm(N, mean=mean(LONS), sd=diff(LONS)/2 ) lat=rnorm(N, mean=mean(LATS), sd=diff(LATS)/2) str1=runif(N,50,100) dip1=runif(N,10, 80) rake1=runif(N,5, 180) dep=runif(N,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) points(XY$x, XY$y) spherefocgeo(MEKS, PROJ, PLOT=TRUE, icut = 3, ndivs = 4, add=TRUE, pal=terrain.colors(100), RECT=TRUE ) ## Not run: plot(x=range(IZ$x), y=range(IZ$y), type='n', asp=1, axes=FALSE, ann=FALSE) image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE) image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE) plotGEOmapXY(haiti.map, LIM = c(Lon.range[1],Lat.range[1] , Lon.range[2] ,Lat.range[2]), PROJ =PROJ, MAPstyle = 2, MAPcol = 'black' , add=TRUE ) H = rectPERIM(JMAT$xo, JMAT$yo) antipolygon(H$x ,H$y, col=grey(.85) , corner=1, pct=.4) sqrTICXY(H , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) spherefocgeo(OLDCMT, PROJ, PLOT=TRUE, add=TRUE, pal=topo.colors(100) ) ## End(Not run)N = 100 LATS = c(7.593004, 25.926045) LONS = c(268.1563 , 305) lon=rnorm(N, mean=mean(LONS), sd=diff(LONS)/2 ) lat=rnorm(N, mean=mean(LATS), sd=diff(LATS)/2) str1=runif(N,50,100) dip1=runif(N,10, 80) rake1=runif(N,5, 180) dep=runif(N,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) points(XY$x, XY$y) spherefocgeo(MEKS, PROJ, PLOT=TRUE, icut = 3, ndivs = 4, add=TRUE, pal=terrain.colors(100), RECT=TRUE ) ## Not run: plot(x=range(IZ$x), y=range(IZ$y), type='n', asp=1, axes=FALSE, ann=FALSE) image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE) image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE) plotGEOmapXY(haiti.map, LIM = c(Lon.range[1],Lat.range[1] , Lon.range[2] ,Lat.range[2]), PROJ =PROJ, MAPstyle = 2, MAPcol = 'black' , add=TRUE ) H = rectPERIM(JMAT$xo, JMAT$yo) antipolygon(H$x ,H$y, col=grey(.85) , corner=1, pct=.4) sqrTICXY(H , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) spherefocgeo(OLDCMT, PROJ, PLOT=TRUE, add=TRUE, pal=topo.colors(100) ) ## End(Not run)
Given a set of points, draw a spline and affix an arrow at the end.
spline.arrow(x, y = 0, kdiv = 20, arrow = 1, length = 0.2, col = "black", thick = 0.01, headlength = 0.2, headthick = 0.1, code = 2, ...)spline.arrow(x, y = 0, kdiv = 20, arrow = 1, length = 0.2, col = "black", thick = 0.01, headlength = 0.2, headthick = 0.1, code = 2, ...)
x |
vector, x-coordinates |
y |
vector, y-coordinates |
kdiv |
Number of divisions |
arrow |
style of arrow, 1=simple arrow, 2=fancy arrow |
length |
length of head |
col |
color of arrow |
thick |
thickness of arrow stem |
headlength |
length of arrow head |
headthick |
thickness of arrow head |
code |
code, 1=arrow on end of spline, 3=arrow on beginning. |
... |
graphical parameters for the line |
Can use either simple arrows or fancy arrows.
list of x,y coordinates of the spline and Graphical Side effects
Jonathan M. Lees<[email protected]>
fancyarrows
plot(c(0,1), c(0,1), type='n') G=list() G$x=c(0.1644,0.1227,0.0659,0.0893,0.2346, 0.3514,0.5518,0.7104,0.6887,0.6903,0.8422) G$y=c(0.8816,0.8305,0.7209,0.6086,0.5372, 0.6061,0.6545,0.6367,0.4352,0.3025,0.0475) spline.arrow(G$x, G$y)plot(c(0,1), c(0,1), type='n') G=list() G$x=c(0.1644,0.1227,0.0659,0.0893,0.2346, 0.3514,0.5518,0.7104,0.6887,0.6903,0.8422) G$y=c(0.8816,0.8305,0.7209,0.6086,0.5372, 0.6061,0.6545,0.6367,0.4352,0.3025,0.0475) spline.arrow(G$x, G$y)
Given a focal mechanism, add Strike Dip lines to a plot.
StrikeDip(x = x, y = y, MEC, focsiz, addDIP = TRUE, ...)StrikeDip(x = x, y = y, MEC, focsiz, addDIP = TRUE, ...)
x |
x-location on plot |
y |
y-location on plot |
MEC |
Focal Mechanism list from SDRFOC |
focsiz |
size of mechanism, inches |
addDIP |
Logical, TRUE = add dip line perpendicular to strike |
... |
graphical parameters |
This is a summary plot to be used instead of Beach Balls.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
Lees, J. M., Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences, Computers & Geosciences, 26, 7, 751-761, 2000.
nipXY, justfocXY, plotmanyfoc
### HAiti Earthquake Jan, 2010 MEC <- SDRfoc(71, 64, 25 , u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) plot(c(0, 1), c(0,1), type='n', asp=1) u <- par("usr") focsiz <- 0.5 justfocXY(MEC, x=.5, y= .5, focsiz=0.5, fcol ='brown' , fcolback = "white", xpd = TRUE) StrikeDip(1.0, .5 , MEC ,focsiz, col="purple", lwd=3 ) nipXY(MEC, x = 0.25, y = .5, focsiz=0.5, fcol ='purple', nipcol = "black", cex = 1) ##### or set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) StrikeDip(XY$x[i], XY$y[i] , MEC ,focsiz, col=jcol, lwd=3 ) }### HAiti Earthquake Jan, 2010 MEC <- SDRfoc(71, 64, 25 , u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) plot(c(0, 1), c(0,1), type='n', asp=1) u <- par("usr") focsiz <- 0.5 justfocXY(MEC, x=.5, y= .5, focsiz=0.5, fcol ='brown' , fcolback = "white", xpd = TRUE) StrikeDip(1.0, .5 , MEC ,focsiz, col="purple", lwd=3 ) nipXY(MEC, x = 0.25, y = .5, focsiz=0.5, fcol ='purple', nipcol = "black", cex = 1) ##### or set.seed(2015) N = 20 lon=runif(20, 268.1563 , 305) lat=runif(20, 7.593004, 25.926045) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) for(i in 1:length(XY$x)) { Msdr = CONVERTSDR(MEKS$str1[i], MEKS$dip1[i],MEKS$rake1[i]) MEC = MRake(Msdr$M) MEC$UP = FALSE jcol = foc.color(foc.icolor(MEC$rake1), pal=1) StrikeDip(XY$x[i], XY$y[i] , MEC ,focsiz, col=jcol, lwd=3 ) }
Illustrate a strikeslip fault using animation
strikeslip.fault(anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 2, Light = c(45, 45))strikeslip.fault(anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 2, Light = c(45, 45))
anim |
animation vector |
KAPPA |
Phong parameter for lighting |
Light |
lighting point |
Program will animate a strikeslip fault for educational purposes. Animation must be stopped by halting execution.
Graphical Side effects
Jonathan M. Lees<[email protected]>
normal.fault, thrust.fault
strikeslip.fault(anim=0, Light=c(45,90) ) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) strikeslip.fault(anim=anim, Light=c(45,90) ) ## End(Not run)strikeslip.fault(anim=0, Light=c(45,90) ) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) strikeslip.fault(anim=anim, Light=c(45,90) ) ## End(Not run)
Plot SV-wave radiation with information from the pickfile and waveform data
SVradfoc(A, MEC, GU, pscale, col)SVradfoc(A, MEC, GU, pscale, col)
A |
Pickfile structure |
MEC |
MEC structure |
GU |
Waveform Event Structure |
pscale |
logical (not used) |
col |
color palette |
Image plot of the SV radiation pattern
Graphical Side effects
Jonathan M. Lees<[email protected]>
imageSV
MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) SVradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )MEC = SDRfoc(65, 32, -34, u=TRUE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE) SVradfoc(NULL, MEC , NULL, TRUE, rainbow(100) )
Create a structure of Tape Base lines
TapeBase()TapeBase()
Program returns the lines and points for plotting a Tape plot. Based on the Hammer projection.
List
The list includes points and other information
Jonathan M. Lees<[email protected]>
Tape, W., and C. Tape (2012), A geometric comparison of source-type plots for moment tensors, Geophys. J. Int., 190, 499-510.
TapePlot, HAMMERprojXY
T1 =TapeBase() TapePlot(T1)T1 =TapeBase() TapePlot(T1)
Tape style Lune Plot using Hammer projection
TapePlot(TapeList = list(), add = FALSE, ann = TRUE, pcol = c(grey(0), grey(0.85), grey(0.95)))TapePlot(TapeList = list(), add = FALSE, ann = TRUE, pcol = c(grey(0), grey(0.85), grey(0.95)))
TapeList |
List of strokes from TapeBase |
add |
logical, TRUE=add to existing plot |
ann |
logical, TRUE=annotape |
pcol |
3-vector of colors: inner lines, upper polygon, lower polygon |
Plot an Tape net from the TapeBase function.
Side effects
Jonathan M. Lees<[email protected]>
Tape, W., and C. Tape (2012), A geometric comparison of source-type plots for moment tensors, Geophys. J. Int., 190, 499-510. https://doi.org/10.1111/j.1365-246X.2012.05490.x
TapeBase, HAMMERprojXY
T1 = TapeBase() TapePlot(T1) data(widdenMoments) WM = widdenMoments par(mfrow=c(1,1), mai=c(0,0,0,0)) T1 = TapeBase() TapePlot(T1) for(i in 1:length(WM$Mxx)) { M1=matrix( c(WM$Mxx[i],WM$Mxy[i],WM$Mxz[i], WM$Mxy[i],WM$Myy[i],WM$Myz[i], WM$Mxz[i],WM$Myz[i],WM$Mzz[i]), ncol=3, byrow=TRUE) E1 = eigen(M1) h = SourceType( sort(E1$values, decreasing=TRUE) ) h$dip = 90-h$phi ## cat(paste(h$dip, h$lam, sep=" "), sep="\n") h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) points(h1$x, h1$y, pch=21, bg="orange" ) }T1 = TapeBase() TapePlot(T1) data(widdenMoments) WM = widdenMoments par(mfrow=c(1,1), mai=c(0,0,0,0)) T1 = TapeBase() TapePlot(T1) for(i in 1:length(WM$Mxx)) { M1=matrix( c(WM$Mxx[i],WM$Mxy[i],WM$Mxz[i], WM$Mxy[i],WM$Myy[i],WM$Myz[i], WM$Mxz[i],WM$Myz[i],WM$Mzz[i]), ncol=3, byrow=TRUE) E1 = eigen(M1) h = SourceType( sort(E1$values, decreasing=TRUE) ) h$dip = 90-h$phi ## cat(paste(h$dip, h$lam, sep=" "), sep="\n") h1 = HAMMERprojXY(h$dip*pi/180, h$lam*pi/180) points(h1$x, h1$y, pch=21, bg="orange" ) }
Plots Beachball figure with numerous vectors and points added and labeled. Useful for teaching about focal mechanisms.
TEACHFOC(s, d, r, up = FALSE)TEACHFOC(s, d, r, up = FALSE)
s |
strike |
d |
dip |
r |
rake |
up |
logical, TRUE = upper |
Graphical side effects
Jonathan M. Lees<[email protected]>
CONVERTSDR, MRake,foc.icolor,focleg, foc.color, focpoint, PlotPlanes, nipXY , fancyarrows
TEACHFOC(65, 32, -34, up=TRUE)TEACHFOC(65, 32, -34, up=TRUE)
Add a point to a ternary plot
ternfoc.point(deltaB, deltaP, deltaT)ternfoc.point(deltaB, deltaP, deltaT)
deltaB |
angle, degrees |
deltaP |
angle, degrees |
deltaT |
angle, degrees |
Plot point on a Ternary diagram using Froelich's algorithm.
List
h |
vector of x coordinates |
v |
vector of y coordinates |
Use Bfocvec(az1, dip1, az2, dip2) to get the deltaB angle.
Jonathan M. Lees <[email protected]>
C. Frohlich. Triangle diagrams: ternary graphs to display similarity and diversity of earthquake focal mechanisms. Physics of the Earth and Planetary Interiors, 75:193-198, 1992.
Bfocvec
Msdr = CONVERTSDR(55.01, 165.65, 29.2 ) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td )Msdr = CONVERTSDR(55.01, 165.65, 29.2 ) MEC = MRake(Msdr$M) MEC$UP = FALSE az1 = Msdr$M$az1 dip1 = Msdr$M$d1 az2 = Msdr$M$az2 dip2 = Msdr$M$d2 BBB = Bfocvec(az1, dip1, az2, dip2) V = ternfoc.point(BBB$Bdip, Msdr$M$pd, Msdr$M$td )
Ternary plots of rake categories (strike-slip, normal, thrust) distributed geographically.
ternfocgeo(CMTSOL, PROJ = NULL, icut = 5, ndivs = 10, bbox=c(0,1, 0, 1), PLOT = TRUE, add = FALSE, RECT = FALSE)ternfocgeo(CMTSOL, PROJ = NULL, icut = 5, ndivs = 10, bbox=c(0,1, 0, 1), PLOT = TRUE, add = FALSE, RECT = FALSE)
CMTSOL |
see output of getCMT for list input |
PROJ |
Map projection |
icut |
cut off for number of points in box, default=5 |
ndivs |
divisions of map area, default=10 |
bbox |
bounding box for dividing the area, given as minX, maxX, minY, maxY; default=usr coordinates from par() |
PLOT |
logical, default=TRUE |
add |
logical, add to existing plot |
RECT |
logical, TRUE=plot rectangles |
Program divides the area into blocks, tests each one for minimum number per block and plots a ternary plot for each block.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
PlotTernfoc, spherefocgeo, prepFOCS, RectDense
N = 100 LATS = c(7.593004, 25.926045) LONS = c(268.1563 , 305) lon=rnorm(N, mean=mean(LONS), sd=diff(LONS)/2 ) lat=rnorm(N, mean=mean(LATS), sd=diff(LATS)/2) str1=runif(N,50,100) dip1=runif(N,10, 80) rake1=runif(N,5, 180) dep=runif(N,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) ## points(XY$x, XY$y) ternfocgeo(MEKS , PROJ, PLOT=TRUE, icut = 3, ndivs = 4, add=TRUE, RECT=TRUE) points(XY$x, XY$y, pch=8, col="purple" ) ################# next restrict the boxes to a specific region plot(range(XY$x), range(XY$y), type='n', asp=1) points(XY$x, XY$y) ternfocgeo(MEKS , PROJ, PLOT=TRUE, icut = 3, ndivs = 5, bbox=c(-2000,2000,-2000,2000) , add=TRUE, RECT=TRUE) ## Not run: ##### this example shows a real application with a map plot(x=range(IZ$x), y=range(IZ$y), type='n', asp=1, axes=FALSE, ann=FALSE) image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE) image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE) plotGEOmapXY(haiti.map, LIM = c(Lon.range[1],Lat.range[1] , Lon.range[2] ,Lat.range[2]), PROJ =PROJ, MAPstyle = 2, MAPcol = 'black' , add=TRUE ) H = rectPERIM(JMAT$xo, JMAT$yo) antipolygon(H$x ,H$y, col=grey(.85) , corner=1, pct=.4) sqrTICXY(H , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ternfocgeo(OLDCMT, PROJ, PLOT=TRUE, add=TRUE) ## End(Not run)N = 100 LATS = c(7.593004, 25.926045) LONS = c(268.1563 , 305) lon=rnorm(N, mean=mean(LONS), sd=diff(LONS)/2 ) lat=rnorm(N, mean=mean(LATS), sd=diff(LATS)/2) str1=runif(N,50,100) dip1=runif(N,10, 80) rake1=runif(N,5, 180) dep=runif(N,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) ## points(XY$x, XY$y) ternfocgeo(MEKS , PROJ, PLOT=TRUE, icut = 3, ndivs = 4, add=TRUE, RECT=TRUE) points(XY$x, XY$y, pch=8, col="purple" ) ################# next restrict the boxes to a specific region plot(range(XY$x), range(XY$y), type='n', asp=1) points(XY$x, XY$y) ternfocgeo(MEKS , PROJ, PLOT=TRUE, icut = 3, ndivs = 5, bbox=c(-2000,2000,-2000,2000) , add=TRUE, RECT=TRUE) ## Not run: ##### this example shows a real application with a map plot(x=range(IZ$x), y=range(IZ$y), type='n', asp=1, axes=FALSE, ann=FALSE) image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE) image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE) plotGEOmapXY(haiti.map, LIM = c(Lon.range[1],Lat.range[1] , Lon.range[2] ,Lat.range[2]), PROJ =PROJ, MAPstyle = 2, MAPcol = 'black' , add=TRUE ) H = rectPERIM(JMAT$xo, JMAT$yo) antipolygon(H$x ,H$y, col=grey(.85) , corner=1, pct=.4) sqrTICXY(H , PROJ, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) ) ternfocgeo(OLDCMT, PROJ, PLOT=TRUE, add=TRUE) ## End(Not run)
Test Right Hand of tensor
testrightHAND(U)testrightHAND(U)
U |
3 by 3 matrix |
The fuction eigen does not always produce a right-handed eigenvector matrix. The code tests each cross product to see if it creates a right-hand system.
logical vector
Jonathan M. Lees<[email protected]>
forcerighthand
Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) testrightHAND(E1$vectors)Mtens <- c(-0.412, 0.084, 0.328 ,0.398, -1.239, 1.058) M1 <- matrix(c(Mtens[1], Mtens[4], Mtens[5], Mtens[4], Mtens[2], Mtens[6], Mtens[5],Mtens[6], Mtens[3]), ncol=3, nrow=3, byrow=TRUE) E1 <- eigen(M1) testrightHAND(E1$vectors)
Illustrate a thrust fault using animation
thrust.fault(anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 2, Light = c(45, 45))thrust.fault(anim = seq(from = 0, to = 1, by = 0.1), KAPPA = 2, Light = c(45, 45))
anim |
animation vector |
KAPPA |
Phong parameter for lighting |
Light |
lighting point |
Program will animate a thrust fault for educational purposes. Animation must be stopped by halting execution.
Graphical Side effects
Jonathan M. Lees<[email protected]>
strikeslip.fault, thrust.fault
thrust.fault(anim=0, KAPPA=4, Light=c(-20, 80)) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) thrust.fault(anim=anim, KAPPA=4, Light=c(-20, 80)) ## End(Not run)thrust.fault(anim=0, KAPPA=4, Light=c(-20, 80)) ## Not run: #### execute a stop command to stop this animation anim= seq(from=0, to=1, by=.1) thrust.fault(anim=anim, KAPPA=4, Light=c(-20, 80)) ## End(Not run)
Tk plot to u-v coordinate transformation
tk2uv(T, k)tk2uv(T, k)
T |
T-value |
k |
k-value |
T and k come from moment tensor analysis.
List: u and v
Keehoon Kim<[email protected]> Jonathan M. Lees<[email protected]>
Hudson
m2tk, hudson.net, hudson.plot
v = c(2,-1,-1) m = m2tk(v) tk2uv(m$T, m$k)v = c(2,-1,-1) m = m2tk(v) tk2uv(m$T, m$k)
Convert cartesian coordinates to strike and dip
to.spherical(x, y, z)to.spherical(x, y, z)
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
LIST
az |
angle, degrees |
dip |
angle, degrees |
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
Jonathan M. Lees <[email protected]>
SDRfoc
to.spherical(3, 4, 5)to.spherical(3, 4, 5)
Convert azimuth and dip to cartesian coordinates
TOCART.DIP(az, dip)TOCART.DIP(az, dip)
az |
azimuth, degrees |
dip |
dip, degrees |
LIST
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
az |
azimuth, degrees |
dip |
dip, degrees |
Jonathan M. Lees <[email protected]>
to.spherical
TOCART.DIP(134, 32)TOCART.DIP(134, 32)
Convert azimuth-dip to cartesian coordinates with list as argument
tocartL(A)tocartL(A)
A |
|
List
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
x positive north, y positive east, z positive downward
Jonathan M. Lees <[email protected]>
TOCART.DIP, RSEIS::TOCART, tosphereL, to.spherical
A = list(az=23, dip=84) tocartL(A)A = list(az=23, dip=84) tocartL(A)
Get Azimuth and Dip from Cartesian vector on a sphere.
TOSPHERE(x, y, z)TOSPHERE(x, y, z)
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
az |
azimuth angle, degrees |
dip |
dip, degrees |
Jonathan M. Lees<[email protected]>
TOSPHERE.DIP, tosphereL, to.spherical
TOSPHERE(3, 4, 5)TOSPHERE(3, 4, 5)
convert to spherical coordinates
TOSPHERE.DIP(x, y, z)TOSPHERE.DIP(x, y, z)
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
takes three components and returns azimuth and dip
List
az |
azimuth, degrees |
dip |
Dip, degrees |
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
Jonathan M. Lees<[email protected]>
to.spherical
TOSPHERE.DIP(3, 4, 5)TOSPHERE.DIP(3, 4, 5)
convert to spherical coordinates
tosphereL(A)tosphereL(A)
A |
list (x,y,z) |
takes list of three components and returns azimuth and dip
List
az |
azimuth, degrees |
dip |
Dip, degrees |
x |
x-coordinate |
y |
y-coordinate |
z |
z-coordinate |
Jonathan M. Lees<[email protected]>
TOSPHERE
A = list(x=12 ,y=2, z=-3 ) tosphereL(A)A = list(x=12 ,y=2, z=-3 ) tosphereL(A)
Convert trend and dip to cartesian coordinates.
TP2XYZ(trend, dip)TP2XYZ(trend, dip)
trend |
trend angle, degrees |
dip |
dip angle, degrees |
These are used as functions auxiallry to rotateFoc.
vector: x, y, z
Jonathan M. Lees<[email protected]>
RotTP
TP2XYZ(34, 40)TP2XYZ(34, 40)
Create a 4 by 4 translation matrix
TRANmat(x, y, z)TRANmat(x, y, z)
x |
x-translation |
y |
y-translation |
z |
z-translation |
Matrix suitaqble for translating a 3D body.
Jonathan M. Lees <[email protected]>
Rogers and Adams, 1990, Mathematical Elements for Computer Graphics, McGraw-Hill, 611p.
ROTX, ROTZ, ROTY
zT = TRANmat(5, 4, 2)zT = TRANmat(5, 4, 2)
Cartesian Moment Tensors from Varvryuk
data(Vmoments)data(Vmoments)
A list of 9 moment tensors from Vaclav Varvryuk
http://www.ig.cas.cz/en/research-&-teaching/software-download/
http://www.ig.cas.cz/en/research-&-teaching/software-download/
Cartesian Moment Tensors from Widden Paper in Utah
data(widdenMoments)data(widdenMoments)
A list of 48 moment tensors from Utah
SRL paper
Seismological Research Letters
plot a Wulff Stereonet (Equal-Angle)
Wnet(add = FALSE, col = gray(0.7), border = "black", lwd = 1)Wnet(add = FALSE, col = gray(0.7), border = "black", lwd = 1)
add |
Logical, TRUE=add to existing plot |
col |
color |
border |
border color |
lwd |
line width |
Plots equal-angle stereonet as opposed to equal-area.
graphical side effects
Jonathan M. Lees <[email protected]>
net, pnet
Wnet()Wnet()
Adds points to Wulff Equal-Angle Stereonet
Wpoint(az1, dip1, col = 2, pch = 5, lab = "", UP = FALSE)Wpoint(az1, dip1, col = 2, pch = 5, lab = "", UP = FALSE)
az1 |
azimuth angle, degrees |
dip1 |
dip angle, degrees |
col |
color |
pch |
plotting character |
lab |
label for point |
UP |
logical, TRUE=Upperhemisphere |
Wulff net point is added to existing plot.
graphical side effects
Jonathan M. Lees <[email protected]>
Wnet
Wnet() Wpoint(23, 34)Wnet() Wpoint(23, 34)
Plot Focal Mechs at X-Y positions on cross sectionsor other plots that do not have geographic coordinates and projection.
xsecmanyfoc(MEK, theta=NULL, focsiz = 0.5, foccol = NULL, UP=TRUE, focstyle=1, LEG = FALSE, DOBAR = FALSE)xsecmanyfoc(MEK, theta=NULL, focsiz = 0.5, foccol = NULL, UP=TRUE, focstyle=1, LEG = FALSE, DOBAR = FALSE)
MEK |
List of Focal Mechanisms, see details |
focsiz |
focal size, inches |
theta |
degrees, angle from north for projecting the focal mechs |
foccol |
focal color, default is to calculate based on rake |
UP |
logical, UP=TRUE means plot upper hemisphere (DEFAULT=TRUE) |
focstyle |
integer, 1=beach ball, 2=nipplot |
LEG |
logical, TRUE= add focal legend for color codes |
DOBAR |
add strike dip bar at epicenter |
Input MEK list contains
MEKS = list(lon=0, lat=0, str1=0, dip1=0, rake1=0, dep=0, name="", Elat=0, Elon=0, x=0, y=0)
The x, y coordinates of the input list are location where the focals will be plotted. For cross sections x=distance along the section and y would be depth. The focal mechs are added to the current plot.
Graphical Side Effects
If theta is NULL focals are plotted as if they were on a plan view. If theta is provided, however, the mechs are plotted with view from the vertical cross section. The cross section is taken at two points. Theta should be determined by viewing the cross section with the first point on the left and the second on the right. The view angle is through the section measured in degrees from north.
Jonathan M. Lees<[email protected]>
Lees, J. M., Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences, Computers & Geosciences, 26, 7, 751-761, 2000.
justfocXY, plotmanyfoc
############# create and plot the mechs in plan view: N = 20 lon=runif(20, 235, 243) lat=runif(20, 45.4, 49) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.5) ex = range(XY$x) why = range(XY$y) JJ = list(x=ex, y=why) SWA = GEOmap::eqswath(XY$x, XY$y, MEKS$dep, JJ, width = diff(why) , PROJ = PROJ) MEKS$x = rep(NA, length(XY$x)) MEKS$y = rep(NA, length(XY$y)) MEKS$x[SWA$flag] = SWA$r MEKS$y[SWA$flag] = -SWA$depth bigR = sqrt( (JJ$x[2]-JJ$x[1])^2 + (JJ$y[2]-JJ$y[1])^2) plot(c(0,bigR) , c(0, min(-SWA$depth)) , type='n', xlab="Distance, KM", ylab="Depth") points(SWA$r, -SWA$depth) xsecmanyfoc(MEKS, focsiz=0.5, LEG = TRUE, DOBAR=FALSE) title("cross section: focals are plotted as if in plan view") ang1 = atan2( JJ$y[2]-JJ$y[1] , JJ$x[2]-JJ$x[1]) degang = ang1*180/pi xsecmanyfoc(MEKS, focsiz=0.5, theta=degang, LEG = TRUE, DOBAR=FALSE) title("cross section: focals are view from the side projection (outer hemisphere)")############# create and plot the mechs in plan view: N = 20 lon=runif(20, 235, 243) lat=runif(20, 45.4, 49) str1=runif(20,50,100) dip1=runif(20,10, 80) rake1=runif(20,5, 180) dep=runif(20,1,15) name=seq(from=1, to=length(lon), by=1) Elat=NULL Elon=NULL yr = rep(2017, times=N) jd = runif(N, min=1, max=365) MEKS = list(lon=lon, lat=lat, str1=str1, dip1=dip1, rake1=rake1, dep=dep, name=name, yr=yr, jd = jd) PROJ = GEOmap::setPROJ(type=2, LAT0=mean(lat) , LON0=mean(lon) ) ## utm XY = GEOmap::GLOB.XY(lat, lon, PROJ) plot(range(XY$x), range(XY$y), type='n', asp=1) plotmanyfoc(MEKS, PROJ, focsiz=0.5) ex = range(XY$x) why = range(XY$y) JJ = list(x=ex, y=why) SWA = GEOmap::eqswath(XY$x, XY$y, MEKS$dep, JJ, width = diff(why) , PROJ = PROJ) MEKS$x = rep(NA, length(XY$x)) MEKS$y = rep(NA, length(XY$y)) MEKS$x[SWA$flag] = SWA$r MEKS$y[SWA$flag] = -SWA$depth bigR = sqrt( (JJ$x[2]-JJ$x[1])^2 + (JJ$y[2]-JJ$y[1])^2) plot(c(0,bigR) , c(0, min(-SWA$depth)) , type='n', xlab="Distance, KM", ylab="Depth") points(SWA$r, -SWA$depth) xsecmanyfoc(MEKS, focsiz=0.5, LEG = TRUE, DOBAR=FALSE) title("cross section: focals are plotted as if in plan view") ang1 = atan2( JJ$y[2]-JJ$y[1] , JJ$x[2]-JJ$x[1]) degang = ang1*180/pi xsecmanyfoc(MEKS, focsiz=0.5, theta=degang, LEG = TRUE, DOBAR=FALSE) title("cross section: focals are view from the side projection (outer hemisphere)")
Create the list structure for a 3D arrow.
Z3Darrow(len = 1, basethick = 0.1, headlen = 0.6, headlip = 0.1)Z3Darrow(len = 1, basethick = 0.1, headlen = 0.6, headlip = 0.1)
len |
Length in user coordinates |
basethick |
Thickness of the base |
headlen |
Length of the head |
headlip |
Width of the overhang lip |
Creates a strucutre suitable for plotting rotated and translated 3D arrows.
List
aglyph |
List of vertices of the faces |
anorm |
Outward facing normal vectors to faces |
Jonathan M. Lees <[email protected]>
PROJ3D, pglyph3D, phong3D
ZA = Z3Darrow(len = 1, basethick = 0.1, headlen = 0.6, headlip = 0.1)ZA = Z3Darrow(len = 1, basethick = 0.1, headlen = 0.6, headlip = 0.1)