MapNonDouble: Map moment tensors

Description

Plot moment tensors on map

Usage

MapNonDouble(Locs, moments, sel = 1, siz = 0.2,
col=rgb(1, .75, .75), PLANES = TRUE, add = FALSE, LEG=FALSE)

Value

list:

FOC: matrix, focal mechanism angles (strike, dip rake)
LAB: matrix, x-y location for labels

Arguments

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

Author

Jonathan M. Lees<jonathan.lees@unc.edu>

Details

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.

References

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.

Examples

Run this code


if (FALSE) {

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)



}

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning