Read Harvard CMT moment
Read and plot a CMT solution copied from the Harvard CMT website.
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<jonathan.lees@unc.edu>
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.