| Title: | Seismic Hypocenter Determination |
|---|---|
| Description: | Non-linear inversion for hypocenter estimation and analysis of seismic data collected continuously, or in trigger mode. The functions organize other functions from 'RSEIS' and 'GEOmap' to help researchers pick, locate, and store hypocenters for detailed seismic investigation. Error ellipsoids and station influence are estimated via jackknife analysis. References include Iversen, E. S., and J. M. Lees (1996)<doi:10.1785/BSSA0860061853>. |
| Authors: | Jonathan M. Lees [aut, cre], Baptiste Auguie [ctb] |
| Maintainer: | Jonathan M. Lees <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.5-1 |
| Built: | 2025-03-01 03:04:46 UTC |
| Source: | https://github.com/cran/Rquake |
Non-linear earthquake locations are estimated by sequential convergence to hypocenter solutions, along with error ellipsoids and 3D-plotting, using a coordination of functions from 'RSEIS', 'GEOmap', 'RFOC' and others for a complete seismic analysis from field campaign data or data extracted from online websites. Interactive codes for seismic phase picking can be combined with event location to go from raw seismic time series to earthquake analysis and spatial statistics.
Rquake is a package for analaysis of seismic data collected continuously, or in trigger mode. The functions organize other functions from 'RSEIS' and 'GEOmap' to help researchers pick, locate, and store hypocenters for detailed seismic investigation.
CONTPF EQXYresid INITpickfile NLSlocate PFoutput RQ SavePF UPdateEQLOC XYSETUP Y2Pphase chak contPFarrivals doAmap gMAP getregionals prepPDE viewCHAC
Jonathan M. Lees<jonathan.lees.edu> Maintainer:Jonathan M. Lees<jonathan.lees.edu>
Lee, W.H.K., and S.W. Stewart, Principles and Applications of Microearthquake Networks, Academic Press, New York, 1981.
library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE)library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE)
1D Velocity Ecuador
data(ASW.vel)data(ASW.vel)
a list of velocities for hypocenter relocation
Mario Ruiz
data(ASW.vel) data(wu_coso.vel) data(fuj1.vel) data(LITHOS.vel) RSEIS::Comp1Dvels(c("ASW.vel","wu_coso.vel", "fuj1.vel" , "LITHOS.vel" ))data(ASW.vel) data(wu_coso.vel) data(fuj1.vel) data(LITHOS.vel) RSEIS::Comp1Dvels(c("ASW.vel","wu_coso.vel", "fuj1.vel" , "LITHOS.vel" ))
Perform jackknife on earthquake location by eliminating stations.
BLACKJACK(Ldat, vel)BLACKJACK(Ldat, vel)
Ldat |
event list |
vel |
Velocity model |
Stations are eliminated, not rows.
event list with pseudo values
events are located with P and S-wave arrivals, but code here should eliminate just stations.
Jonathan M. Lees<[email protected]>
Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.
Vlocate, plotJACKLLZ
###### lps=list of files names to be read data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) B = BLACKJACK(Ldat, vel) ## the code HiJACK ### runs BLACKJACK on many pickfiles stored in files ### COSOjack = HiJACK(lps, sta) ### plotJACKLLZ(COSOjack, sta, proj)###### lps=list of files names to be read data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) B = BLACKJACK(Ldat, vel) ## the code HiJACK ### runs BLACKJACK on many pickfiles stored in files ### COSOjack = HiJACK(lps, sta) ### plotJACKLLZ(COSOjack, sta, proj)
Check to see if location data has the minimally correct list components.
checkLOCATEinput(Ldat, EQ, vel = NULL)checkLOCATEinput(Ldat, EQ, vel = NULL)
Ldat |
list, must inlude: x,y,err, sec, cor (see details) |
EQ |
list, must inlude: x,y,z, t |
vel |
list, 1D velocity structure |
Input pick list must have at x,y,z, sec, cor, err elements for each station.
logical: FALSE mean problem with data
Jonathan M. Lees<[email protected]>
XYlocate
library(RSEIS) library(GEOmap) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) #### get station X-Y values in km XY = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon, proj) ### add to Ldat list Ldat$x = XY$x Ldat$y = XY$y wstart = which.min(Ldat$sec) EQ = list(x=XY$x[wstart], y=XY$y[wstart], z=6, t=Ldat$sec[wstart] ) checkLOCATEinput(Ldat, EQ)library(RSEIS) library(GEOmap) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) #### get station X-Y values in km XY = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon, proj) ### add to Ldat list Ldat$x = XY$x Ldat$y = XY$y wstart = which.min(Ldat$sec) EQ = list(x=XY$x[wstart], y=XY$y[wstart], z=6, t=Ldat$sec[wstart] ) checkLOCATEinput(Ldat, EQ)
Given a pick file in WPX format, break the picks apart clustered accoring to single link cluster analysis.
clusterWPX(twpx, tol = 200, PLOT = FALSE)clusterWPX(twpx, tol = 200, PLOT = FALSE)
twpx |
WPX list |
tol |
tolerance in seconds - all pick distances less than tol will be set to zero to force these to be associated. |
PLOT |
logical, if TRUE, add verbose plotting |
If there is not significant separation of picks, only one cluster is returned. To avoid spurious clusters, increase the tolerance.
list of WPX lists
Cluster depends on what one considers a cluster.
Jonathan M. Lees<[email protected]>
RSEIS::addWPX, RSEIS::catWPX, RSEIS::checkWPX,RSEIS::cleanWPX, PCsaveWPX, RSEIS::setWPX, RSEIS::repairWPX
s1 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s2 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=5, mi=2, sec = runif(5)) s3 = RSEIS::catWPX(s1,s2) twpx = data.frame(s3) L3 = clusterWPX(twpx)s1 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s2 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=5, mi=2, sec = runif(5)) s3 = RSEIS::catWPX(s1,s2) twpx = data.frame(s3) L3 = clusterWPX(twpx)
Button to Contour Pickfile Arrivals, used internally in swig.
CONTPF(nh, g, idev = 3)CONTPF(nh, g, idev = 3)
nh |
RSEIS list |
g |
swig parameters |
idev |
device for plotting |
Driver for contPFarrivals
Side effects
Jonathan M. Lees<[email protected]>
contPFarrivals
if(interactive()){ ###### interactive: addition of button in swig data(GH, package='RSEIS') buts = "CONTPF" RSEIS::swig(GH, PADDLAB=buts, SHOWONLY=FALSE ) }if(interactive()){ ###### interactive: addition of button in swig data(GH, package='RSEIS') buts = "CONTPF" RSEIS::swig(GH, PADDLAB=buts, SHOWONLY=FALSE ) }
Contour plot of arrival times recorded in a pickfile list.
contPFarrivals(PF, stas, proj=NULL, cont=TRUE, POINTS=TRUE, image=FALSE , col=RSEIS::tomo.colors(50), gcol="black", phase="P", add=TRUE)contPFarrivals(PF, stas, proj=NULL, cont=TRUE, POINTS=TRUE, image=FALSE , col=RSEIS::tomo.colors(50), gcol="black", phase="P", add=TRUE)
PF |
Pickfile list in RSEIS format |
stas |
station list |
proj |
projection from GEOmap |
cont |
logical, add contour to plot |
POINTS |
logical, add mark up (stations) to plot |
image |
logical, add image to plot |
col |
color palette for image |
gcol |
color for contour lines |
phase |
character, phase to contour |
add |
logical, TRUE=add to existing plot |
Contours the arrival time. The earliest arrival is subtracted from each time pick. Uses only the phase indicated and there can be only one phase per station - default is earliest at each station.
Graphical Side Effects
Jonathan M. Lees<[email protected]>
doAmap
library(RSEIS) data(GH, package='RSEIS') sta = GH$stafile g1 = GH$pickfile proj = GEOmap::setPROJ(type=2, LAT0 =median(sta$lat) , LON0 = median(sta$lon)) grcol = grey(seq(from=0.3, to=0.95, length=50)) contPFarrivals(g1, sta, proj=proj,cont=TRUE, POINTS=TRUE, image=TRUE , col=grcol, phase="P", add=FALSE )library(RSEIS) data(GH, package='RSEIS') sta = GH$stafile g1 = GH$pickfile proj = GEOmap::setPROJ(type=2, LAT0 =median(sta$lat) , LON0 = median(sta$lon)) grcol = grey(seq(from=0.3, to=0.95, length=50)) contPFarrivals(g1, sta, proj=proj,cont=TRUE, POINTS=TRUE, image=TRUE , col=grcol, phase="P", add=FALSE )
Coso Station Location file, 1989-1999
data(coso_sta_LLZ)data(coso_sta_LLZ)
Name, Lat, Lon, Z
Personal Files
Wu, H. and J. M. Lees (1996). Attenuation Structure of Coso Geothermal Area, California, from P Wave Pulse Widths, Bull. Seismol. Soc. Am., 86, 1574-1590.
Lees, J. M. (1998), Multiplet analysis at Coso Geothermal,Bull. Seismol. Soc. Am. 88(5) 1127-1143.
Set of selected seismic arrival files with hypocenter locations.
data("cosopix")data("cosopix")
List consisting of:
PF: original text version of file, as read from disk
AC: Acard: hypocenter information
LOC: location
MC: Fault Mechanizm card
STAS: Station information
LIP: Error Ellipse
E: E-card
F: F-card
filename: original file location
UWFILEID: UW file identification
comments: Comments on event location
OSTAS: Station names
H: High resolution location numbers
N: Stations Not used in location
Each element of this list is an individual earthquake record.
data(cosopix) A = sapply(cosopix, '[[', 'LOC') ### gather stations ST.name = vector(mode='character') ST.lat = vector(mode='numeric') ST.lon = vector(mode='numeric') ST.z = vector(mode='numeric') for(i in 1:length(cosopix)) { g = cosopix[[i]] g = data.frame(g$STAS ) w = which(!is.na(g$lat) ) ST.name = c(ST.name, g$name[w]) ST.lat = c(ST.lat, g$lat[w]) ST.lon = c(ST.lon, g$lon[w]) ST.z = c(ST.z, g$z[w]) } notdup = !duplicated(ST.name) name = ST.name[notdup ] lat = ST.lat[notdup ] lon =ST.lon[notdup ] z = ST.z[notdup ] plot(range(c(A[9, ], lon)) , range(c(A[8, ], lat)) , type='n', xlab='Lon', ylab='Lat') points(lon, lat, pch=6) text(lon, lat, labels=name, pos=3) points(A[9, ], A[8, ])data(cosopix) A = sapply(cosopix, '[[', 'LOC') ### gather stations ST.name = vector(mode='character') ST.lat = vector(mode='numeric') ST.lon = vector(mode='numeric') ST.z = vector(mode='numeric') for(i in 1:length(cosopix)) { g = cosopix[[i]] g = data.frame(g$STAS ) w = which(!is.na(g$lat) ) ST.name = c(ST.name, g$name[w]) ST.lat = c(ST.lat, g$lat[w]) ST.lon = c(ST.lon, g$lon[w]) ST.z = c(ST.z, g$z[w]) } notdup = !duplicated(ST.name) name = ST.name[notdup ] lat = ST.lat[notdup ] lon =ST.lon[notdup ] z = ST.z[notdup ] plot(range(c(A[9, ], lon)) , range(c(A[8, ], lat)) , type='n', xlab='Lon', ylab='Lat') points(lon, lat, pch=6) text(lon, lat, labels=name, pos=3) points(A[9, ], A[8, ])
Default Velocity Function is returned in the event no velocity function is available.
defaultVEL(kind = 1)defaultVEL(kind = 1)
kind |
integer, 1=fuj1, 2=LITHOS |
A set of default velocity functions are available.
velocity list, P and S waves
Jonathan M. Lees<[email protected]>
fuj1.vel
v = defaultVEL(1)v = defaultVEL(1)
Distance weighting for non-linear earthquake location.
DistWeight(dist, err, distwt) DistWeightLL(lat, lon, elat, elon, err, distwt) DistWeightXY(x, y, ex, ey, err, distwt)DistWeight(dist, err, distwt) DistWeightLL(lat, lon, elat, elon, err, distwt) DistWeightXY(x, y, ex, ey, err, distwt)
dist |
distance in km |
err |
sigma error in seconds |
distwt |
distance weighting parameter |
lat |
Latitude |
lon |
Longitude |
elat |
Event Latitude |
elon |
Event Longitude |
x |
station X(km) |
y |
station Y(km) |
ex |
event X (km) |
ey |
event Y (km) |
Based on Lquake scheme from University of Washington. If you need to reduce the effect of distance weighting, increase distwt.
Since the hypocenter moves between each iteration, the distance weighting is updated.
vector of weights
Jonathan M. Lees<[email protected]>
DistWeight(1:10, .4, 20)DistWeight(1:10, .4, 20)
Plot a map of station locations
doAmap(stas, doproj = TRUE)doAmap(stas, doproj = TRUE)
stas |
station list |
doproj |
logical, if TRUE, project (UTM) the data so plot is in units of km with the median lat-lon as the center. If FALSE, use the lat-lon coordinates. |
The range of the plot is expanded by 10 percent prior to plotting.
list, GEOmap projection
Jonathan M. Lees<[email protected]>
gMAP,expandbound,GLOB.XY
data(coso_sta_LLZ) ### or read in from file: ## fsta = "staLLZ.txt" ## stas = scan(file=fsta,what=list(name="", lat=0, lon=0, z=0)) ## stas$z = stas$z/1000 stas = coso_sta_LLZ STA = doAmap(stas, doproj = TRUE)data(coso_sta_LLZ) ### or read in from file: ## fsta = "staLLZ.txt" ## stas = scan(file=fsta,what=list(name="", lat=0, lon=0, z=0)) ## stas$z = stas$z/1000 stas = coso_sta_LLZ STA = doAmap(stas, doproj = TRUE)
Error Elipse for Hypocenter Location
eqlipse(x, y, cov, wcols = c(1, 2), dof = 2, pct=0.05, ...)eqlipse(x, y, cov, wcols = c(1, 2), dof = 2, pct=0.05, ...)
x |
X-location for drawing |
y |
Y-location for drawing |
cov |
matrix, 3 by 3 Covariance matrix |
wcols |
vector, which columns to extract from cov, see details. |
dof |
Degrees of Freedom for 95 percent confidence |
pct |
Percent used for 2-sided confidence bounds, default=0.05 |
... |
graphical parameters, par |
The 3 by 3 matrix is supplied and a 2 by 2 matrix is subtracted depending on which components are being drawn. For X-Y projections, use wcols=c(1,2). For vertical cross sections, rotate the cov matrix and then extract the columns.
Side effects, graphical
Jonathan M. Lees<[email protected]>
eqwrapup
library(RSEIS) data(GH, package='RSEIS') data(VELMOD1D, package='RSEIS' ) vel = VELMOD1D gpf = GH$pickfile w1 = which(gpf$STAS$phase=="P" | gpf$STAS$phase=="S" ) N = length(w1) Ldat = list( name = gpf$STAS$name[w1], sec = gpf$STAS$sec[w1], phase = gpf$STAS$phase[w1], lat=gpf$STAS$lat[w1], lon = gpf$STAS$lon[w1], z = gpf$STAS$z[w1], err= gpf$STAS$err[w1], yr = rep(gpf$LOC$yr , times=N), jd = rep(gpf$LOC$jd, times=N), mo = rep(gpf$LOC$mo, times=N), dom = rep(gpf$LOC$dom, times=N), hr =rep( gpf$LOC$hr, times=N), mi = rep(gpf$LOC$mi, times=N) ) EQ = GH$pickfile$LOC EQ$t = EQ$sec kuality = eqwrapup(Ldat, EQ, vel, distwt = 20, verbose = TRUE ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) XYSTAS = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon , proj) eqxy = GEOmap::GLOB.XY(EQ$lat, EQ$lon, proj) plot(range(c(XYSTAS$x, eqxy$x)), range(c(XYSTAS$y, eqxy$y)), type='n', asp=1, xlab="km", ylab="km" ) points(XYSTAS$x, XYSTAS$y, pch=6) points(eqxy$x, eqxy$y, pch=8, col='red') #### covariance matrix KOV = kuality$cov[2:4, 2:4] #### add uncertainty eqlipse(eqxy$x, eqxy$y , KOV, wcols = c(1,2) , dof=kuality$ndf, border="blue" )library(RSEIS) data(GH, package='RSEIS') data(VELMOD1D, package='RSEIS' ) vel = VELMOD1D gpf = GH$pickfile w1 = which(gpf$STAS$phase=="P" | gpf$STAS$phase=="S" ) N = length(w1) Ldat = list( name = gpf$STAS$name[w1], sec = gpf$STAS$sec[w1], phase = gpf$STAS$phase[w1], lat=gpf$STAS$lat[w1], lon = gpf$STAS$lon[w1], z = gpf$STAS$z[w1], err= gpf$STAS$err[w1], yr = rep(gpf$LOC$yr , times=N), jd = rep(gpf$LOC$jd, times=N), mo = rep(gpf$LOC$mo, times=N), dom = rep(gpf$LOC$dom, times=N), hr =rep( gpf$LOC$hr, times=N), mi = rep(gpf$LOC$mi, times=N) ) EQ = GH$pickfile$LOC EQ$t = EQ$sec kuality = eqwrapup(Ldat, EQ, vel, distwt = 20, verbose = TRUE ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) XYSTAS = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon , proj) eqxy = GEOmap::GLOB.XY(EQ$lat, EQ$lon, proj) plot(range(c(XYSTAS$x, eqxy$x)), range(c(XYSTAS$y, eqxy$y)), type='n', asp=1, xlab="km", ylab="km" ) points(XYSTAS$x, XYSTAS$y, pch=6) points(eqxy$x, eqxy$y, pch=8, col='red') #### covariance matrix KOV = kuality$cov[2:4, 2:4] #### add uncertainty eqlipse(eqxy$x, eqxy$y , KOV, wcols = c(1,2) , dof=kuality$ndf, border="blue" )
Calculate error and summary information on earthquake location.
eqwrapup(Ldat, EQ, vel, distwt=20, lambdareg = 0.0, verbose=FALSE)eqwrapup(Ldat, EQ, vel, distwt=20, lambdareg = 0.0, verbose=FALSE)
Ldat |
List of station arrival times, lat-lon, and uncertainty |
EQ |
List of earthquake location: Lat-Lon-z-t |
vel |
velocity model |
distwt |
distance weight, default=20 |
lambdareg |
numeric, regularization parameter (default=0) |
verbose |
logical, TRUE=print information to screen |
Earthquakes are located with a generalized inverse (SVD). covariance matrix is extracted and 95% confidence bounds are calculated. Quality factors Q1 and Q1 estimate the quality iof the location based on the gap, minimum distance and rms.
List
rms |
Root Mean Square Residual |
meanres |
Mean Residual |
sdres |
Standard Dev of residuals |
sdmean |
Standard error of mean residual |
sswres |
Sum squared weighted residuals |
ndf |
Number of Degrees of Freedom |
sterrx |
km, error in X (East-West) |
sterry |
km, error in Y (North-South) |
sterrz |
km, error in Z, (depth) |
sterrt |
s, Delta-time |
cov |
covariance matrix (used for error ellipsoids) |
lam |
lambda |
gap |
Spatial gap (max subtended angle) |
herr |
Horizontal error |
distmin |
Minimum distance to epicenter |
Q1 |
Quality Factor based on Gap and RMS |
Q2 |
Quality factor based on RMS, depth and min-Distance |
The Damping parameter (lambda) is set to zero. In the UW lquake program, lambda is set to 0.02.
Jonathan M. Lees<[email protected]>
Klocate, Glocate, getGAP
library(RSEIS, package='RSEIS') data(GH, package='RSEIS') data(wu_coso.vel, package='Rquake' ) vel = wu_coso.vel gpf = GH$pickfile w1 = which(gpf$STAS$phase=="P" | gpf$STAS$phase=="S" ) N = length(w1) Ldat = list( name = gpf$STAS$name[w1], sec = gpf$STAS$sec[w1], phase = gpf$STAS$phase[w1], lat=gpf$STAS$lat[w1], lon = gpf$STAS$lon[w1], z = gpf$STAS$z[w1], err= gpf$STAS$err[w1], yr = rep(gpf$LOC$yr , times=N), jd = rep(gpf$LOC$jd, times=N), mo = rep(gpf$LOC$mo, times=N), dom = rep(gpf$LOC$dom, times=N), hr =rep( gpf$LOC$hr, times=N), mi = rep(gpf$LOC$mi, times=N) ) EQ = GH$pickfile$LOC EQ$t = EQ$sec kuality = eqwrapup(Ldat, EQ, vel, distwt = 20, verbose = TRUE ) names(kuality)library(RSEIS, package='RSEIS') data(GH, package='RSEIS') data(wu_coso.vel, package='Rquake' ) vel = wu_coso.vel gpf = GH$pickfile w1 = which(gpf$STAS$phase=="P" | gpf$STAS$phase=="S" ) N = length(w1) Ldat = list( name = gpf$STAS$name[w1], sec = gpf$STAS$sec[w1], phase = gpf$STAS$phase[w1], lat=gpf$STAS$lat[w1], lon = gpf$STAS$lon[w1], z = gpf$STAS$z[w1], err= gpf$STAS$err[w1], yr = rep(gpf$LOC$yr , times=N), jd = rep(gpf$LOC$jd, times=N), mo = rep(gpf$LOC$mo, times=N), dom = rep(gpf$LOC$dom, times=N), hr =rep( gpf$LOC$hr, times=N), mi = rep(gpf$LOC$mi, times=N) ) EQ = GH$pickfile$LOC EQ$t = EQ$sec kuality = eqwrapup(Ldat, EQ, vel, distwt = 20, verbose = TRUE ) names(kuality)
Given an earthquake hypocenter and a list of station information, retrieve the station residuals.
EQXYresid(XY, vel = list(), h1 = c(0, 0, 0, 0), PLOT = FALSE)EQXYresid(XY, vel = list(), h1 = c(0, 0, 0, 0), PLOT = FALSE)
XY |
matrix of station location and arrival times. |
vel |
list, RSEIS velocity model |
h1 |
hypocenter location, c(x,y,z,t) |
PLOT |
logical, TRUE=plot the residuals |
The XY mtrix is in cartesian coordinates, i.e. it has been projected into units of km. Only 1D velocity models are used at this time. Only residuals of P and S wave arrivals are estimated.
vector, right hand side of the least squares problem.
Jonathan M. Lees<[email protected]>
travel.time1D,UPdateEQLOC
#### get sample data data(GH, package='RSEIS') pstas = GH$pickfile ###### get velocity file v = GH$velfile #### project to flatten proj = GEOmap::setPROJ(type = 2, LAT0 = mean(pstas$STAS$lat), LON0 = mean(pstas$STAS$lon) ) XY = GEOmap::GLOB.XY(pstas$STAS$lat, pstas$STAS$lon, proj) ####### elevation corrections elcor = rep(0, length(pstas$STAS$lat)) DZ = pstas$STAS$z - mean(pstas$STAS$z) elcor[pstas$STAS$phase=="P"] = DZ[pstas$STAS$phase=="P"]/v$vp[1] elcor[pstas$STAS$phase=="S"] = DZ[pstas$STAS$phase=="S"]/v$vs[1] ###### set up requisite vectors XY$cor = elcor XY$phase = pstas$STAS$phase XY$sec = pstas$STAS$sec sol = c(GH$pickfile$LOC$lat, GH$pickfile$LOC$lon, GH$pickfile$LOC$z, GH$pickfile$LOC$sec) eqXY = GEOmap::GLOB.XY(sol[1], sol[2], proj) ####### get residuals res = EQXYresid(XY, vel=v , h1=c(eqXY$x, eqXY$y, sol[3], sol[4] ) , PLOT=FALSE)#### get sample data data(GH, package='RSEIS') pstas = GH$pickfile ###### get velocity file v = GH$velfile #### project to flatten proj = GEOmap::setPROJ(type = 2, LAT0 = mean(pstas$STAS$lat), LON0 = mean(pstas$STAS$lon) ) XY = GEOmap::GLOB.XY(pstas$STAS$lat, pstas$STAS$lon, proj) ####### elevation corrections elcor = rep(0, length(pstas$STAS$lat)) DZ = pstas$STAS$z - mean(pstas$STAS$z) elcor[pstas$STAS$phase=="P"] = DZ[pstas$STAS$phase=="P"]/v$vp[1] elcor[pstas$STAS$phase=="S"] = DZ[pstas$STAS$phase=="S"]/v$vs[1] ###### set up requisite vectors XY$cor = elcor XY$phase = pstas$STAS$phase XY$sec = pstas$STAS$sec sol = c(GH$pickfile$LOC$lat, GH$pickfile$LOC$lon, GH$pickfile$LOC$z, GH$pickfile$LOC$sec) eqXY = GEOmap::GLOB.XY(sol[1], sol[2], proj) ####### get residuals res = EQXYresid(XY, vel=v , h1=c(eqXY$x, eqXY$y, sol[3], sol[4] ) , PLOT=FALSE)
Given three angles return rotation matrix.
euler_passive(phi, theta, psi)euler_passive(phi, theta, psi)
phi |
angle with x-axis |
theta |
angle with y-axis |
psi |
angle with z-axis |
Code borrowed from cpp code in package cda. used in rgl.ellipsoid.
3 by 3 rotation matrix.
Jonathan M. Lees<[email protected]>, Baptiste Auguie<[email protected]>
rgl.ellipsoid
options(rgl.useNULL = TRUE) phi=30*pi/180 ; theta= 20*pi/180; psi = 6*pi/180 rr = euler_passive(phi,theta,psi)options(rgl.useNULL = TRUE) phi=30*pi/180 ; theta= 20*pi/180; psi = 6*pi/180 rr = euler_passive(phi,theta,psi)
Given a covariance matrix calculated with Vlocate, extract euler's angles for plotting in rgl
getEulers(R)getEulers(R)
R |
covarince matrix |
Extract the euler angles for plotting an ellipsoid. psi about X-axis, theta about Y axis, phi about Z-axis.
vector, phi theta psi
Used in conjunction with ROTcovQUAKE
Jonathan M. Lees<[email protected]>
ROTcovQUAKE
options(rgl.useNULL = TRUE) R = matrix( runif(9), ncol=3) getEulers(R)options(rgl.useNULL = TRUE) R = matrix( runif(9), ncol=3) getEulers(R)
Given an earthquake and a set of stations, return the maximum angle subtended between adjacent stations relative to the epicenter.
getGAP(EQ, Ldat, PLOT = FALSE)getGAP(EQ, Ldat, PLOT = FALSE)
EQ |
List, Earthequake location, elements (lat, lon) must be present |
Ldat |
List, station information, (lat, lon) must be present |
PLOT |
logical, plot the stations and show the gap |
Theangles are calculated in cartesian coordinates with the epicenter at the origin using a UTM projection.
numeric, gap in degrees
Jonathan M. Lees<[email protected]>
eqwrapup
set.seed(0) N = 10 snames = paste(sep="", "A", as.character(1:N)) stas = list(name=snames, lat=runif(N, 35.9823, 36.1414), lon=runif(N, -118.0031, -117.6213)) NEQ = 3 WEQ = list(lat=runif(NEQ, 35.9823, 36.1414), lon=runif(NEQ, -118.0031, -117.6213)) MLAT = median(stas$lat) MLON = median(stas$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) XYSTAS = GEOmap::GLOB.XY(stas$lat, stas$lon , proj) eqxy = GEOmap::GLOB.XY(WEQ$lat, WEQ$lon, proj) plot(range(c(XYSTAS$x, eqxy$x)), range(c(XYSTAS$y, eqxy$y)), type='n', asp=1, xlab="km", ylab="km" ) points(XYSTAS$x, XYSTAS$y, pch=6) for(i in 1:NEQ) { EQ = list(lat=WEQ$lat[i], lon=WEQ$lon[i]) g = getGAP(EQ, stas, PLOT=FALSE) points(eqxy$x[i], eqxy$y[i], pch=8, col='red') text(eqxy$x[i], eqxy$y[i], labels=paste("gap=", format(g)), pos=3) }set.seed(0) N = 10 snames = paste(sep="", "A", as.character(1:N)) stas = list(name=snames, lat=runif(N, 35.9823, 36.1414), lon=runif(N, -118.0031, -117.6213)) NEQ = 3 WEQ = list(lat=runif(NEQ, 35.9823, 36.1414), lon=runif(NEQ, -118.0031, -117.6213)) MLAT = median(stas$lat) MLON = median(stas$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) XYSTAS = GEOmap::GLOB.XY(stas$lat, stas$lon , proj) eqxy = GEOmap::GLOB.XY(WEQ$lat, WEQ$lon, proj) plot(range(c(XYSTAS$x, eqxy$x)), range(c(XYSTAS$y, eqxy$y)), type='n', asp=1, xlab="km", ylab="km" ) points(XYSTAS$x, XYSTAS$y, pch=6) for(i in 1:NEQ) { EQ = list(lat=WEQ$lat[i], lon=WEQ$lon[i]) g = getGAP(EQ, stas, PLOT=FALSE) points(eqxy$x[i], eqxy$y[i], pch=8, col='red') text(eqxy$x[i], eqxy$y[i], labels=paste("gap=", format(g)), pos=3) }
Get Pand S travel times and derivatives
GETpsTT(phase, eqz = 6, staz = 0, delx = 1, dely = 1, deltadis = 6, vel)GETpsTT(phase, eqz = 6, staz = 0, delx = 1, dely = 1, deltadis = 6, vel)
phase |
character vector, phase |
eqz |
event depth |
staz |
station elevation |
delx |
km, delta X |
dely |
km, delta Y |
deltadis |
km, distance |
vel |
velocity models (P and S) |
Creates a vector of travel times, and a matrix and derivatives used for inversion.
list:
TT |
travel time vector |
Derivs |
matrix of derivatives, dtdx, dtdy, dtdz |
Jonathan M. Lees<[email protected]>
many.time1D
library(RSEIS) library(GEOmap) data(GH, package='RSEIS') data(VELMOD1D, package='RSEIS') vel = VELMOD1D p1 = GH$pickfile$STAS loc = GH$pickfile$LOC proj = GEOmap::setPROJ(type = 2, LAT0 =loc$lat, LON0 = loc$lon) XYsta = GEOmap::GLOB.XY(p1$lat, p1$lon, proj) XYq = GEOmap::GLOB.XY(loc$lat, loc$lon, proj) delx = XYq$x-XYsta$x dely = XYq$y-XYsta$y dists = sqrt(delx^2+dely^2) G1 = GETpsTT(p1$phase, eqz=loc$z, staz=0, delx=delx, dely=dely, deltadis=dists , vel)library(RSEIS) library(GEOmap) data(GH, package='RSEIS') data(VELMOD1D, package='RSEIS') vel = VELMOD1D p1 = GH$pickfile$STAS loc = GH$pickfile$LOC proj = GEOmap::setPROJ(type = 2, LAT0 =loc$lat, LON0 = loc$lon) XYsta = GEOmap::GLOB.XY(p1$lat, p1$lon, proj) XYq = GEOmap::GLOB.XY(loc$lat, loc$lon, proj) delx = XYq$x-XYsta$x dely = XYq$y-XYsta$y dists = sqrt(delx^2+dely^2) G1 = GETpsTT(p1$phase, eqz=loc$z, staz=0, delx=delx, dely=dely, deltadis=dists , vel)
Extract regional events from a hypocenter list (catalog)
getregionals(KAT, Mlat, Mlon, rad = 1000, t1 = 1, t2 = 2)getregionals(KAT, Mlat, Mlon, rad = 1000, t1 = 1, t2 = 2)
KAT |
catalog list: must include lat, lon and jsec. |
Mlat |
central latitude |
Mlon |
central longitude |
rad |
radius (km) |
t1 |
start time (julian days) |
t2 |
end time (julian days) |
Given an earthquake catalog from PDEs, for example, extract the events that are close to a network in a given time frame. The limited data set may be used to help predict arrival times for known hypocenter locations.
The time jsec is in julian days, i.e. jsec=jd+hr/24+mi/(24*60)+sec/(24*3600) so that they can be compared to t1 and t2.
Catalog
Jonathan M. Lees<[email protected]>
RSEIS::Mine.seis, RSEIS::swig
set.seed(1) Mlat = 36.00833 Mlon = -117.8048 N = 100 degz = 5 KAT = list(lat=runif(N, Mlat-degz,Mlat+degz) , lon=runif(N,Mlon-degz,Mlon+degz) ) ###### ranfdom times in January KAT$jsec = runif(N, 1, 30) + runif(N, 0, 24)/(24) + runif(N, 0, 59)/(24*60) ###### extract regional events localeqs = getregionals(KAT, Mlat, Mlon, rad=200 , t1=NULL, t2=NULL) plot(KAT$lon, KAT$lat, pch=8, col=grey(0.75) ) points(KAT$lon[localeqs], KAT$lat[localeqs], pch=1, col='red', cex=1.5 )set.seed(1) Mlat = 36.00833 Mlon = -117.8048 N = 100 degz = 5 KAT = list(lat=runif(N, Mlat-degz,Mlat+degz) , lon=runif(N,Mlon-degz,Mlon+degz) ) ###### ranfdom times in January KAT$jsec = runif(N, 1, 30) + runif(N, 0, 24)/(24) + runif(N, 0, 59)/(24*60) ###### extract regional events localeqs = getregionals(KAT, Mlat, Mlon, rad=200 , t1=NULL, t2=NULL) plot(KAT$lon, KAT$lat, pch=8, col=grey(0.75) ) points(KAT$lon[localeqs], KAT$lat[localeqs], pch=1, col='red', cex=1.5 )
Given an earthquake location and a set of arrival times, return a vector of residuals.
getresidTT(Ldat, EQ, stas, vel)getresidTT(Ldat, EQ, stas, vel)
Ldat |
List of arrival times |
EQ |
List of event location, (lat, lon, z, and time) |
stas |
station location list |
vel |
list, velocity structure |
1D travel time calculation.
vector of residuals
Jonathan M. Lees<[email protected]>
travel.time1D
######### LF is a vector of arrival time files ##### KAM is a set of locations data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) zip = LeftjustTime(twpx) w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) resids = getresidTT(Ldat, g1$LOC, g1$STAS , vel)######### LF is a vector of arrival time files ##### KAM is a set of locations data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) zip = LeftjustTime(twpx) w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) resids = getresidTT(Ldat, g1$LOC, g1$STAS , vel)
Extract the lat lon from the pick file.
Gfirstguess(Ldat, type = "first")Gfirstguess(Ldat, type = "first")
Ldat |
Matrix of data information |
type |
one of "first", "mean", or "median" |
Either the earliest arrival or the average station is returned. Used internally in the earthquake location program to provide a first guess.
vector, lat, lon, z and tee
Jonathan M. Lees<[email protected]>
Klocate
data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) g1 = Gfirstguess(twpx, type = "first")data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) g1 = Gfirstguess(twpx, type = "first")
Generic Map Button
gMAP(nh, g, idev = 3)gMAP(nh, g, idev = 3)
nh |
RSEIS structure |
g |
parameters used in swig |
idev |
device for plotting (not used) |
This is a button used internally in swig
Graphical Side Effects
Jonathan M. Lees<[email protected]>
swig
if(interactive()){ #### this is interactive ### adds button to swig menu data(GH, package='RSEIS') buts = "gMAP" RSEIS::swig(GH, PADDLAB = buts ) }if(interactive()){ #### this is interactive ### adds button to swig menu data(GH, package='RSEIS') buts = "gMAP" RSEIS::swig(GH, PADDLAB = buts ) }
defining functions for swig
GPIX(nh, g)GPIX(nh, g)
nh |
waveform list for RSEIS |
g |
plotting parameter list for interactive program |
Buttons can be defined on the fly.
Multiple picks on a panel
The return value depends on the nature of the function as it is returned to the main code swig. Choices for returning to swig are: break, replot, revert, replace, donothing, exit.
Jonathan M. Lees<[email protected]>
swig, XTR
if(interactive()){ ###### interactive addition of buttons in swig STDLAB=c("DONE", "QUIT", "SELBUT" , "GPIX" ) data(GH, package='RSEIS') JJ = RSEIS::swig(GH, sel=1:10, STDLAB=STDLAB) }if(interactive()){ ###### interactive addition of buttons in swig STDLAB=c("DONE", "QUIT", "SELBUT" , "GPIX" ) data(GH, package='RSEIS') JJ = RSEIS::swig(GH, sel=1:10, STDLAB=STDLAB) }
Jackknife a list of events
HiJACK(lps, sta, vel)HiJACK(lps, sta, vel)
lps |
list of earthquake event pickfiles, each element is an individual pickfile list, with STAS: relative timing of phase arrivals |
sta |
staiton list |
vel |
velocity list |
Driver for BLACKJACK
jackknife pseudovalues for each event
Jonathan M. Lees<[email protected]>
Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.
BLACKJACK
##### uses external files, runs Vlocate on each one #### lps = list of file names to be read data(cosopix) data(wu_coso.vel) data(coso_sta_LLZ) COSOjack = HiJACK(cosopix, coso_sta_LLZ, wu_coso.vel) proj = GEOmap::setPROJ(2, mean(coso_sta_LLZ$lat), mean(coso_sta_LLZ$lon)) #### show stats plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=1 ) #### show maps plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=2 )##### uses external files, runs Vlocate on each one #### lps = list of file names to be read data(cosopix) data(wu_coso.vel) data(coso_sta_LLZ) COSOjack = HiJACK(cosopix, coso_sta_LLZ, wu_coso.vel) proj = GEOmap::setPROJ(2, mean(coso_sta_LLZ$lat), mean(coso_sta_LLZ$lon)) #### show stats plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=1 ) #### show maps plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=2 )
Plot contours/image of Influence scores.
imageINFLUENCE(B, sta, proj)imageINFLUENCE(B, sta, proj)
B |
Pseudovalue list |
sta |
station location list |
proj |
projection list |
Following jackknife - plot results. this function is called by plotJACKLLZ.
side effects
Jonathan M. Lees<[email protected]>
Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.
plotJACKLLZ
Initialize a pickfile
INITpickfile(stas = NULL, src = NULL, WPX = NULL)INITpickfile(stas = NULL, src = NULL, WPX = NULL)
stas |
station list |
src |
hypocenter location |
WPX |
GPIX or PPIX picks from swig |
Initialize a pickfile with a set of picks extracted from swig.
list, pickfile
Jonathan M. Lees<[email protected]>
EmptyPickfile
data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) PF = INITpickfile(stas=GH$stafile, src=NULL, WPX=WW )data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) PF = INITpickfile(stas=GH$stafile, src=NULL, WPX=WW )
Earthquake Hypocenter Location
Klocate(Ldat, sol = c(0, 0, 0, 0), vel=defaultVEL(6), distwt = 20, errtol = c(0.01, 0.01, 0.01), maxit = 20, Lambda = 1, guessdepth = 6, APLOT = FALSE, stas = list(name = "", lat = NA, lon = NA, z = NA))Klocate(Ldat, sol = c(0, 0, 0, 0), vel=defaultVEL(6), distwt = 20, errtol = c(0.01, 0.01, 0.01), maxit = 20, Lambda = 1, guessdepth = 6, APLOT = FALSE, stas = list(name = "", lat = NA, lon = NA, z = NA))
Ldat |
swig pick list |
sol |
vector, initial solution |
vel |
velocity list |
distwt |
distance weight parameter |
errtol |
error tolerance |
maxit |
Maximum number of iterations |
Lambda |
damping parameter |
guessdepth |
initial depth for guess |
APLOT |
logical, plot intermediate solutions |
stas |
station list |
Inversion is done with SVD.
Event location in Lat-Lon-Z-T.
Damped least squares.
Jonathan M. Lees<[email protected]>
swig, defaultVEL
###### could read from a list of files on disk ### LF = list.files(path=pdir, pattern="p$", full.names=TRUE ) data(GH, package='RSEIS') g1 = GH$pickfile ## points(g1$H$lon, g1$H$lat, pch=8, col='red') w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) ###### let the code determine the initial guess NEW = Klocate(Ldat )###### could read from a list of files on disk ### LF = list.files(path=pdir, pattern="p$", full.names=TRUE ) data(GH, package='RSEIS') g1 = GH$pickfile ## points(g1$H$lon, g1$H$lat, pch=8, col='red') w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) ###### let the code determine the initial guess NEW = Klocate(Ldat )
'RSEIS' Button: Restore Last WPX file from memory. Function is used internally in swig.
lastPIX(nh, g) editPIX(nh, g)lastPIX(nh, g) editPIX(nh, g)
nh |
GH list from RSEIS |
g |
parameters from swig |
New WPX list attached to g
Jonathan M. Lees<[email protected]>
Given an existing list of seismic picks, add Latitude, Longitude and Elevation associated with the indicated station.
latlonz2wpx(twpx, stas)latlonz2wpx(twpx, stas)
twpx |
List of picks from swig |
stas |
station list |
The names of the stations are matched to the station names int he station file.
Pick file with LLZ added as list members.
Jonathan M. Lees<[email protected]>
Klocate
data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS )data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS )
List location data
LDATlist(g1, w1)LDATlist(g1, w1)
g1 |
loc list |
w1 |
index |
side effects
Jonathan M. Lees<[email protected]>
Adjust times relative to least minute.
LeftjustTime(g1)LeftjustTime(g1)
g1 |
list with times, yr, jd, hr, mi, sec |
Reutrns the list with the times adjusted to the least minimum (left adjusted)
list is returned.
Jonathan M. Lees<[email protected]>
recdate
set.seed(0) d1 = list(yr=rep(2005, 4), jd=rep(5, 4), hr=rep(6, 4), mi=c(1,1,2,3), sec=runif(4, 0, 60)) LeftjustTime(d1)set.seed(0) d1 = list(yr=rep(2005, 4), jd=rep(5, 4), hr=rep(6, 4), mi=c(1,1,2,3), sec=runif(4, 0, 60)) LeftjustTime(d1)
Check WPX list for legitimate picks
legitWPX(twpx, quiet=TRUE)legitWPX(twpx, quiet=TRUE)
twpx |
pick information list (WPX) |
quiet |
logical, default=TRUE, FALSE generates an error message |
Used internall to test if a WPX list has legitimate picks. Initially a list is generated with NA and 0 values in the place holders. If no legitimate picks are added, the list still exists, but the picks are bogus, so this routine will return 0.
integer: 0=not legitimate, 1=legitimate
Currently only the name is tested for all(NA), but this might be changed int he future for a more sophisticated test.
Jonathan M. Lees<[email protected]>
PCsaveWPX
### test fails library(RSEIS) jk = RSEIS::cleanWPX() legitWPX(jk) #### test passes: data(GH, package='RSEIS') gwpx = RSEIS::uwpfile2ypx(GH$pickfile) legitWPX(gwpx)### test fails library(RSEIS) jk = RSEIS::cleanWPX() legitWPX(jk) #### test passes: data(GH, package='RSEIS') gwpx = RSEIS::uwpfile2ypx(GH$pickfile) legitWPX(gwpx)
calculate the mean km distance of a set of Lat-lon pairs
MeanStaDist(Ldat)MeanStaDist(Ldat)
Ldat |
station list with elements of Lat-Lon |
Given a list with elements named lat and lon, find the mean station distance.
scalar
Jonathan M. Lees<[email protected]>
setPROJ, GLOB.XY, dist
data(GH, package='RSEIS') MeanStaDist(GH$pickfile$STAS)data(GH, package='RSEIS') MeanStaDist(GH$pickfile$STAS)
Nonlinear Least Squares Location using Gieger's method
NLSlocate(GH, vel = list(), init = c(0, 0, 0, 0), PLOT = FALSE)NLSlocate(GH, vel = list(), init = c(0, 0, 0, 0), PLOT = FALSE)
GH |
List, RSEIS |
vel |
velocity model |
init |
initial guess for event location |
PLOT |
logical, TRUE=plot |
This is an adaptation of non-linear least squares inversion for earthquake location. A residual function is supplied, and iterations are performed until the location is determined.
vector, new location
At this stage there are no weighting mechanisms or code to eliminate data that has residuals that are too large.
Jonathan M. Lees<[email protected]>
Lee, W.H.K., and S.W. Stewart, Principles and Applications of Microearthquake Networks, Academic Press, New York, 1981.
swig
data(GH, package='RSEIS') ### location is: eqsol = NLSlocate(GH, vel=GH$velfile, PLOT=TRUE )data(GH, package='RSEIS') ### location is: eqsol = NLSlocate(GH, vel=GH$velfile, PLOT=TRUE )
Require only one pick per station of a specified phase.
OnePerSta(twpx, phase = "Y")OnePerSta(twpx, phase = "Y")
twpx |
WPX list |
phase |
character, specific phase |
This is used to reduce the number of picks for specific station and phase. The purpose is avoid multiple P-wave phases for each station in the earthquake location routines.
WPX list
For S-waves there may be multiple S-wave arrivals, as in the case for shear wave splitting. In that case it is probably best to name the phases differently, as in S1, S2, for example.
Jonathan M. Lees<[email protected]>
cleanWPX, repairWPX
s1 = RSEIS::setWPX(name="HI", phase="P", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s2 = RSEIS::setWPX(name="BYE", phase="P", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s3 = RSEIS::catWPX(s1, s2) s4 = OnePerSta(s3, phase = "P")s1 = RSEIS::setWPX(name="HI", phase="P", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s2 = RSEIS::setWPX(name="BYE", phase="P", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) s3 = RSEIS::catWPX(s1, s2) s4 = OnePerSta(s3, phase = "P")
Create a character string from a date for naming unique output files.
PCfiledatetime(orgtim, tims)PCfiledatetime(orgtim, tims)
orgtim |
time vector of length 5: c(yr, jd, hr, mi, sec) |
tims |
seconds to add to orgtim |
filename |
character string |
Jonathan M. Lees<[email protected]>
library(RSEIS) data(GH, package='RSEIS') g1 = getGHtime(GH) g2 = unlist(g1) PCfiledatetime(g2, 1)library(RSEIS) data(GH, package='RSEIS') g1 = getGHtime(GH) g2 = unlist(g1) PCfiledatetime(g2, 1)
Save a WPX list to a file on the local file system.
PCsaveWPX(twpx, destdir = NULL)PCsaveWPX(twpx, destdir = NULL)
twpx |
WPX list |
destdir |
character, destination directory, default=NULL |
Creates a file with the list as in native binary format. This file can be loaded with the standard load function in R. The name of the file is created by using the minimum time extracted from the WPX list. The suffix on the file name is RDATA. When reading in, the object created is named "twpx" for further processing.
destdir must be set, otherwise the destination directory will be temporary. Typically this is set to a local directory where the user has write access.
Side effects on file system. The name of the output file is returned.
User must have write access to the destination directory.
Jonathan M. Lees<[email protected]>
RSEIS::addWPX, RSEIS::catWPX, RSEIS::checkWPX, RSEIS::cleanWPX, RSEIS::clusterWPX, RSEIS::repairWPX, RSEIS::setWPX
##### save files as RDS to users disk s1 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) hh = PCsaveWPX(s1, destdir= tempdir() ) ### read in the data twpx = readRDS(hh) data.frame(twpx)##### save files as RDS to users disk s1 = RSEIS::setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) hh = PCsaveWPX(s1, destdir= tempdir() ) ### read in the data twpx = readRDS(hh) data.frame(twpx)
Write a pickfile to disk, after updating the earthquake location, in a variety of formats.
PFoutput(PF, stas = NULL, sol = NULL, format = 0, destdir=NULL)PFoutput(PF, stas = NULL, sol = NULL, format = 0, destdir=NULL)
PF |
Pickfile list from RSEIS |
stas |
station list |
sol |
solution vector, (lat, lon, z, t0) |
format |
integer, 0=all formats, 1=native R, 2=UW, 3=csv) |
destdir |
character, full path to destination directory, default=NULL ) |
Writes files to disk in local directory.
Side effects: writes files to user's disk
The destdir (destination directory) must be provided for the file to be save properly.
Creates a file name and writes to disk in a variety of formats.
A destdir that is NULL will result in writing to a temporary file.
Jonathan M. Lees<[email protected]>
SavePF, RSEIS
data(GH, package='RSEIS') g1 = GH$pickfile #### saves pick files to disk PFoutput(g1, stas = NULL, sol = NULL, format = 1, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 2, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 3, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 0, destdir=tempdir() )data(GH, package='RSEIS') g1 = GH$pickfile #### saves pick files to disk PFoutput(g1, stas = NULL, sol = NULL, format = 1, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 2, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 3, destdir=tempdir() ) PFoutput(g1, stas = NULL, sol = NULL, format = 0, destdir=tempdir() )
Picking functions for swig
Pick3(nh, g)Pick3(nh, g)
nh |
waveform list for RSEIS |
g |
plotting parameter list for interactive program |
Buttons can be defined on the fly.
Multiple picks on a panel
The return value depends on the nature of the function as it is returned to the main code swig. Choices for returning to swig are: break, replot, revert, replace, donothing, exit.
Jonathan M. Lees<[email protected]>
swig, PickWin
if(interactive()){ ###### interactive addition of button in swig library(RSEIS) MYFUNC<-function(nh, g) { print("pressed MYFUNC") d = data.frame(list(stations=nh$STNS, components=nh$COMPS)) print(d) g$action = "replot" invisible(list(global.vars=g)) } STDLAB=c("DONE", "QUIT", "SELBUT" , "MYFUNC" ) data(GH, package='RSEIS') JJ = RSEIS::swig(GH, sel=1:10, STDLAB=STDLAB) }if(interactive()){ ###### interactive addition of button in swig library(RSEIS) MYFUNC<-function(nh, g) { print("pressed MYFUNC") d = data.frame(list(stations=nh$STNS, components=nh$COMPS)) print(d) g$action = "replot" invisible(list(global.vars=g)) } STDLAB=c("DONE", "QUIT", "SELBUT" , "MYFUNC" ) data(GH, package='RSEIS') JJ = RSEIS::swig(GH, sel=1:10, STDLAB=STDLAB) }
Plot Earthquake location
plotEQ(Ldat, AQ, add = FALSE, prep = FALSE, TIT = "UTM Projected Stations", proj = NULL, xlim = NULL, ylim = NULL)plotEQ(Ldat, AQ, add = FALSE, prep = FALSE, TIT = "UTM Projected Stations", proj = NULL, xlim = NULL, ylim = NULL)
Ldat |
Data list |
AQ |
Earthquake solution (location) |
add |
logical, TRUE=add to plot |
prep |
preparation |
TIT |
title |
proj |
projection list |
xlim |
2-vector, x limits (km) |
ylim |
2-vector, y limits (km) |
used internally in RElocateEQ
graphical side effects
Jonathan M. Lees<[email protected]>
RElocateEQ
BoxPlot Jackknife of station locations
plotJACKLLZ(hjack, sta, proj = NULL, PLOT=1)plotJACKLLZ(hjack, sta, proj = NULL, PLOT=1)
hjack |
Output of hijack |
sta |
station location list |
proj |
projection list |
PLOT |
plotting flag, 1,2. if plot =1 plot only boxplot, if plot=2 plot only map. Default=0 |
takes the output of the HiJack function and extracts the pseudovalues and influence information for boxplots.
Graphical side effects and
X |
influence of lon |
Y |
influence of lat |
Z |
influence of depth |
Jonathan M. Lees<[email protected]>
Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.
HiJACK, BLACKJACK,imageINFLUENCE, Vlocate
data(cosopix) data(wu_coso.vel) data(coso_sta_LLZ) COSOjack = HiJACK(cosopix, coso_sta_LLZ, wu_coso.vel) proj = GEOmap::setPROJ(2, mean(coso_sta_LLZ$lat), mean(coso_sta_LLZ$lon)) #### show stats plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=1 ) #### show maps plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=2 )data(cosopix) data(wu_coso.vel) data(coso_sta_LLZ) COSOjack = HiJACK(cosopix, coso_sta_LLZ, wu_coso.vel) proj = GEOmap::setPROJ(2, mean(coso_sta_LLZ$lat), mean(coso_sta_LLZ$lon)) #### show stats plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=1 ) #### show maps plotJACKLLZ(COSOjack, coso_sta_LLZ, proj, PLOT=2 )
Post Processing on EQrquake
PostREQquake(XQ, proj)PostREQquake(XQ, proj)
XQ |
List of Earthquakes |
proj |
projection list |
Following event locations, plot.
graphical side effects
Jonathan M. Lees<[email protected]>
Plotting error ellipsoids of many events
PostVquake(MANYeq, GX, GY, XY, proj, add=FALSE, ...)PostVquake(MANYeq, GX, GY, XY, proj, add=FALSE, ...)
MANYeq |
List of earthquakes following Vlocate |
GX |
X-bounds for plot |
GY |
Y-bounds for plot |
XY |
station locations in km |
proj |
projection list |
add |
logical; if TRUE, add to existing plot (DEFAULT=FALSE) |
... |
graphical parameters for plotting (see par) |
Plots the event and the error ellipsoids
Graphical side effects
This is used to plot many event locations and their error ellipsoids
Jonathan M. Lees<[email protected]>
eqlipse
Return the range of dates and times for any list with a date/time list
Qrangedatetime(D)Qrangedatetime(D)
D |
info list from RSEIS seismic data list |
min |
date time list |
max |
date time list |
Jonathan M. Lees<[email protected]>
library(RSEIS) data(GH, package='RSEIS') v = Qrangedatetime(GH$info)library(RSEIS) data(GH, package='RSEIS') v = Qrangedatetime(GH$info)
Button to reset the choices of station and component in swig and Mine.seis
ReSet(nh, g)ReSet(nh, g)
nh |
RSEIS list |
g |
swig parameters |
Driver for SELstaDB
Side effects
Jonathan M. Lees<[email protected]>
SELstaDB, Mine.seis
if(interactive()){ data(GH, package='RSEIS') buts = "ReSet" RSEIS::swig(GH, PADDLAB=buts) }if(interactive()){ data(GH, package='RSEIS') buts = "ReSet" RSEIS::swig(GH, PADDLAB=buts) }
Extract Event location information following Vlocate
ripper(AQ)ripper(AQ)
AQ |
event location list |
Extract lat-lon from event locations to track intermediate solutions and convergence
2 by N matrix, lat-lon
Jonathan M. Lees<[email protected]>
plotEQ
library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE) qtip = ripper(AQ)library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE) qtip = ripper(AQ)
Selects which rows in the hypocenter determination to keep during non-linear itaerations based on robust rsidual elimination.
Rowz2Keep(Ldat, EQ, G1, RESMAX)Rowz2Keep(Ldat, EQ, G1, RESMAX)
Ldat |
List of station arrivals |
EQ |
Earthquake location |
G1 |
derivative and travel time estimates |
RESMAX |
2-vector for P and S-wave residual maxima |
This is a utility used internally.
Residuals greater than the respective maxima provided are eliminated in the svd inversion. If fewer than 4 remain, the smallest 4 rows are returned.
Index of good rows
Jonathan M. Lees<[email protected]>
XYlocate
Driver for NLSlocate
RQ(nh, g, idev = 3)RQ(nh, g, idev = 3)
nh |
RSEIS list |
g |
parameters from swig |
idev |
device for plotting |
Button to be called from within swig after picking.
new hypocenter
Jonathan M. Lees<[email protected]>
NLSlocate, EQXYresid, XYSETUP, swig,chak
if(interactive()){ ##### interactive data(GH, package='RSEIS') buts = c("GPIX","PPIX", "PickWin", "fspread", "gMAP", "RQ" , "CONTPF") RSEIS::swig(GH, PADDLAB=buts) }if(interactive()){ ##### interactive data(GH, package='RSEIS') buts = c("GPIX","PPIX", "PickWin", "fspread", "gMAP", "RQ" , "CONTPF") RSEIS::swig(GH, PADDLAB=buts) }
Save a pick file from within swig
SavePF(nh, g)SavePF(nh, g)
nh |
RSEIS data list |
g |
list of parameters internal to swig |
Uses PFoutput to save a pickfile to disk.
Side Effects
Pickfile is saved as a native R file with wpx extension
Jonathan M. Lees<[email protected]>
PFoutput
if(interactive()){ data(GH, package='RSEIS') buts = "SavePF" RSEIS::swig(GH, PADDLAB=buts) }if(interactive()){ data(GH, package='RSEIS') buts = "SavePF" RSEIS::swig(GH, PADDLAB=buts) }
Pick stations and components interactively. This is a routine used in swig.
SELstaDB(IDB, sel=1, newdev=TRUE, STAY=FALSE)SELstaDB(IDB, sel=1, newdev=TRUE, STAY=FALSE)
IDB |
list with component vectors, usta and ucomp |
sel |
vector of index to selected traces |
newdev |
logical, whether to create a new device. |
STAY |
logical, whether to keep device active. |
vector of index to list of stations and components
Jonathan M. Lees<[email protected]>
infoDB, makeDB
if(interactive()){ ### make a database from the files on disk ### DBnov = makeDB(fpath, fpat, kind=2, Iendian=1, BIGLONG=FALSE) ### IDB = infoDB(DBnov) ### or, as an example: data(GH, package='RSEIS') DBnov = list(usta = unique(GH$STNS), unique(GH$COMPS)) k = SELstaDB(IDB) }if(interactive()){ ### make a database from the files on disk ### DBnov = makeDB(fpath, fpat, kind=2, Iendian=1, BIGLONG=FALSE) ### IDB = infoDB(DBnov) ### or, as an example: data(GH, package='RSEIS') DBnov = list(usta = unique(GH$STNS), unique(GH$COMPS)) k = SELstaDB(IDB) }
Update an Earthquake location following a relocation.
UPdateEQLOC(PF, sol, vel, stas = NULL)UPdateEQLOC(PF, sol, vel, stas = NULL)
PF |
Pickfile List |
sol |
solution vector (lat, lon, z, t0) |
vel |
1D velocity model |
stas |
station list (name, lat, lon, z) |
After re-picking or changing the model or the station corrections, update the event location in the pickfile.
Pickfile List
Jonathan M. Lees<[email protected]>
EQXYresid, NLSlocate,PFoutput
data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE) sol = c(AQ$EQ$lat, AQ$EQ$lon, AQ$EQ$z, AQ$EQ$t) upf = UPdateEQLOC(g1, sol , vel, stas=g1$STAS)data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE) sol = c(AQ$EQ$lat, AQ$EQ$lon, AQ$EQ$z, AQ$EQ$t) upf = UPdateEQLOC(g1, sol , vel, stas=g1$STAS)
Hypocenter Determination with error checking and adjustments.
Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100, REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.1, toly = 0.1, tolz = 0.5, RESMAX = c(.4,.5), maxITER = c(7, 5, 7, 4), PLOT=FALSE)Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100, REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.1, toly = 0.1, tolz = 0.5, RESMAX = c(.4,.5), maxITER = c(7, 5, 7, 4), PLOT=FALSE)
Ldat |
list, must inlude: lat, lon ,err, sec, cor (see details) |
EQ |
list, must inlude: lat,lon,z, t |
vel |
list, 1D velocity structure |
distwt |
distance weighting factor |
lambdareg |
regularization parameter for damping |
REG |
logical, TRUE=use regularization |
WTS |
logical, TRUE==use weighting |
STOPPING |
logical, TRUE=use stopping criteria |
tolx |
numeric, tolerance in km in x direction |
toly |
numeric, tolerance in km in y direction |
tolz |
numeric, tolerance in km in z direction |
RESMAX |
vector, residual max for P and S, default=c(4,5) |
maxITER |
vector, Maximum number of iterations for each section of the location routine, default=c(7,5,7,4) |
PLOT |
logical, plot results during iterations |
This is a wrapper for XYlocate, only here the lat-lon of the stations is passed and the code does the projection internally.
There are 3 main loops, each controled by differing input params: first event is located only in XY keeping the depth fixed (7 iterations). Then an initial free solution is estimated using robust elimination of residual based on RESMAX (5 iterations). Finally a set of 7 iterations is applied providing the final estimate, along with error bars, elliposids, etc.
In the event no good solution is derived, the regularization parameter is doubled and a loop with 4 iterations is applied, and the result returned.
list:
EQ |
Hypocenter lcoation |
ERR |
Error Analysis |
its |
number of iteration |
Ksolutions |
list of matrices, each with intermediate x,y,z,t locations |
The schedule may be adjusted by duplicating this function and changing the maxit parameters.
Jonathan M. Lees<[email protected]>
Lee and Stewart
XYlocate, Klocate, DoRLocate
library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE)library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) wstart = which.min(Ldat$sec) EQ = list(lat=Ldat$lat[wstart], lon=Ldat$lon[wstart], z=6, t=Ldat$sec[wstart] ) AQ = Vlocate(Ldat,EQ,vel, distwt = 10, lambdareg =100 , REG = TRUE, WTS = TRUE, STOPPING = TRUE, tolx = 0.01, toly = 0.01 , tolz = 0.05, maxITER = c(7,5,7,4) , RESMAX = c(0.1, 0.1), PLOT=FALSE)
Error Bars in X and Y
XYerror.bars(x, y, xlo = 0, xhi = 0, ylo = 0, yhi = 0, pch = 1, col = 1, barw = 0.1, add = FALSE, ...)XYerror.bars(x, y, xlo = 0, xhi = 0, ylo = 0, yhi = 0, pch = 1, col = 1, barw = 0.1, add = FALSE, ...)
x |
X-values |
y |
Y-values |
xlo |
X Lower limit of error bars |
xhi |
X Upper limit of error bars |
ylo |
Y Lower limit of error bars |
yhi |
Y Upper limit of error bars |
pch |
plotting character |
col |
color |
barw |
width of the bar (inches) |
add |
logical, add=FALSE starts a new plot |
... |
other plotting parameters |
graphical side effects
Jonathan M. Lees<[email protected]>
set.seed(0) zup = rnorm(10) x = 1:10 y = 2*x+5+zup ydown = rnorm(10) ydown = ydown-min(ydown)+.2 yup = rnorm(10) yup = yup-min(yup)+.2 zup = rnorm(10) xup = zup-min(zup)+.5 xdown = rnorm(10) xdown = xdown-min(xdown)+.2 #### example with different error on either side: XYerror.bars(x, y, y-ydown, y+yup, x-xdown, x+xup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE)set.seed(0) zup = rnorm(10) x = 1:10 y = 2*x+5+zup ydown = rnorm(10) ydown = ydown-min(ydown)+.2 yup = rnorm(10) yup = yup-min(yup)+.2 zup = rnorm(10) xup = zup-min(zup)+.5 xdown = rnorm(10) xdown = xdown-min(xdown)+.2 #### example with different error on either side: XYerror.bars(x, y, y-ydown, y+yup, x-xdown, x+xup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE)
Non-linear hypocenter location with UTM geographical projection. Used for locating earthquakes in local or regional settings.
XYlocate(Ldat, EQ, vel, maxITER = 10, distwt = 10, lambdareg = 100, FIXZ = FALSE, REG = TRUE, WTS = TRUE, STOPPING = TRUE, RESMAX = c(.4,.5), tolx = 0.005, toly = 0.005, tolz = 0.01, PLOT = FALSE)XYlocate(Ldat, EQ, vel, maxITER = 10, distwt = 10, lambdareg = 100, FIXZ = FALSE, REG = TRUE, WTS = TRUE, STOPPING = TRUE, RESMAX = c(.4,.5), tolx = 0.005, toly = 0.005, tolz = 0.01, PLOT = FALSE)
Ldat |
list, must inlude: x,y,err, sec, cor (see details) |
EQ |
list, must inlude: x,y,z, t |
vel |
list, 1D velocity structure |
maxITER |
Maximum number of iterations |
distwt |
distance weighting factor |
lambdareg |
regularization parameter for damping |
FIXZ |
logical, TRUE = fix depth, i.e. only calculate x,y,t |
REG |
logical, TRUE=use regularization |
WTS |
logical, TRUE==use weighting |
STOPPING |
logical, TRUE=use stopping criteria |
RESMAX |
vector, residual max for P and S, default=c(4,5) |
tolx |
numeric, tolerance in km in x direction |
toly |
numeric, tolerance in km in y direction |
tolz |
numeric, tolerance in km in z direction |
PLOT |
logical, plot results during iterations |
Input pick list must have at x,y,z, sec, cor, err elements for each station. If no station correction is available it is set to zero. If no uncertainty (err) is available, it is set to 0.05 sec. Each station must have a finite x-y coordinate and arrival time in seconds. Events are located relative to the minute.
Routine uses the svd in a sequence of linear inversions to estimate the nonlinear location.
List:
EQ |
list, Earthquake hypocenter and time |
its |
number of iterations |
rms |
rms residual |
wrms |
wheighted rms residual |
used |
vector, index of used equations |
guesses |
list of x,y,z,t intermediate locations when converging |
This routine should be called by a wrapper (Vlocate) that applies the algorithm several times and changes parameters based on the quality.
If RESMAX is used and the robust approach yields fewer than 4 equations, the best (smallest) four residuals will be used to determiine the event location.
Jonathan M. Lees<[email protected]>
Vlocate
library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) #### get station X-Y values in km XY = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon, proj) ### add to Ldat list Ldat$x = XY$x Ldat$y = XY$y wstart = which.min(Ldat$sec) EQ = list(x=XY$x[wstart], y=XY$y[wstart], z=6, t=Ldat$sec[wstart] ) maxITER = 7 ###print(EQ) AQ = XYlocate(Ldat,EQ,vel, maxITER = maxITER, distwt = 1, lambdareg =10 , FIXZ = FALSE, REG = TRUE, WTS = TRUE, STOPPING = TRUE, RESMAX = c(0.1,0.1), tolx = 0.001, toly = 0.001 , tolz = 0.5, PLOT=FALSE) ######## update the new location AXY = GEOmap::XY.GLOB(AQ$EQ$x, AQ$EQ$y, proj) AQ$EQ$lat = AXY$lat AQ$EQ$lon = AXY$lon if(AQ$EQ$lon>180) { AQ$EQ$lon = AQ$EQ$lon-360 } plot(c(Ldat$x, AQ$EQ$x) , c(Ldat$y,AQ$EQ$y), type='n' , xlab="km", ylab="km" ) points(Ldat$x, Ldat$y, pch=6) points(AQ$EQ$x, AQ$EQ$y, pch=8, col='red') points(EQ$x, EQ$y, pch=4, col='blue') legend("topright", pch=c(8,4, 6), col=c("red", "blue", "black"), legend=c("Final location", "Initial guess", "Station")) print(AQ) ##### try a different case with an extremely wrong start EQ$x = 10 EQ$y = 2library(RSEIS) data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D w1 = which(!is.na(g1$STAS$lat)) sec = g1$STAS$sec[w1] N = length(sec) Ldat = list( name = g1$STAS$name[w1], sec = g1$STAS$sec[w1], phase = g1$STAS$phase[w1], lat=g1$STAS$lat[w1], lon = g1$STAS$lon[w1], z = g1$STAS$z[w1], err= g1$STAS$err[w1], yr = rep(g1$LOC$yr , times=N), jd = rep(g1$LOC$jd, times=N), mo = rep(g1$LOC$mo, times=N), dom = rep(g1$LOC$dom, times=N), hr =rep( g1$LOC$hr, times=N), mi = rep(g1$LOC$mi, times=N) ) MLAT = median(Ldat$lat) MLON = median(Ldat$lon) proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON) #### get station X-Y values in km XY = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon, proj) ### add to Ldat list Ldat$x = XY$x Ldat$y = XY$y wstart = which.min(Ldat$sec) EQ = list(x=XY$x[wstart], y=XY$y[wstart], z=6, t=Ldat$sec[wstart] ) maxITER = 7 ###print(EQ) AQ = XYlocate(Ldat,EQ,vel, maxITER = maxITER, distwt = 1, lambdareg =10 , FIXZ = FALSE, REG = TRUE, WTS = TRUE, STOPPING = TRUE, RESMAX = c(0.1,0.1), tolx = 0.001, toly = 0.001 , tolz = 0.5, PLOT=FALSE) ######## update the new location AXY = GEOmap::XY.GLOB(AQ$EQ$x, AQ$EQ$y, proj) AQ$EQ$lat = AXY$lat AQ$EQ$lon = AXY$lon if(AQ$EQ$lon>180) { AQ$EQ$lon = AQ$EQ$lon-360 } plot(c(Ldat$x, AQ$EQ$x) , c(Ldat$y,AQ$EQ$y), type='n' , xlab="km", ylab="km" ) points(Ldat$x, Ldat$y, pch=6) points(AQ$EQ$x, AQ$EQ$y, pch=8, col='red') points(EQ$x, EQ$y, pch=4, col='blue') legend("topright", pch=c(8,4, 6), col=c("red", "blue", "black"), legend=c("Final location", "Initial guess", "Station")) print(AQ) ##### try a different case with an extremely wrong start EQ$x = 10 EQ$y = 2
Set up matrix for hypocenter inversion
XYSETUP(STAS, init, vel)XYSETUP(STAS, init, vel)
STAS |
station information from pickfile |
init |
initial event location |
vel |
list, velocity |
This sets up the matrix used for nonlinear inversion. The code does not include information on the weighting. Station corrections are included.
The STAS are an internal component of the pickfile.
matrix
Need scheme for weighting according to errors in picks and distance weighting.
Jonathan M. Lees<[email protected]>
setPROJ, GLOB.XY,NLSlocate
## start with the location of the closest station data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D STAS = GH$pickfile$STAS w1 = STAS$phase == 'P' initz = 6 t0a = GH$pickfile$LOC$sec XY = XYSETUP(STAS, c(STAS$lat[w1],STAS$lon[w1], initz, STAS$sec[w1]-t0a ) , vel )## start with the location of the closest station data(GH, package='RSEIS') g1 = GH$pickfile data(VELMOD1D, package='RSEIS') vel= VELMOD1D STAS = GH$pickfile$STAS w1 = STAS$phase == 'P' initz = 6 t0a = GH$pickfile$LOC$sec XY = XYSETUP(STAS, c(STAS$lat[w1],STAS$lon[w1], initz, STAS$sec[w1]-t0a ) , vel )
Removes extraneous other-phase from a pick file. If Ypix were made initially as a rough pick, this removes them.
Y2Pphase(twpx, phase)Y2Pphase(twpx, phase)
twpx |
WPX list |
phase |
character, phase to exchange to P |
Initially many events may be picked using GPIX button. These should be removed after the P-phases have been determined with PickWin.
WPX returned without other-phases
Jonathan M. Lees<[email protected]>
PPIX, GPIX, YPIX, PickWin
data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) twpx$phase[twpx$phase=='P'] = 'Y' #### now twpx is like a Ypix from swig ### switch to P newwpx = Y2Pphase(twpx, "Y" )data(GH, package='RSEIS') WW = RSEIS::uwpfile2ypx(GH$pickfile) twpx = latlonz2wpx(WW, GH$pickfile$STAS ) twpx$phase[twpx$phase=='P'] = 'Y' #### now twpx is like a Ypix from swig ### switch to P newwpx = Y2Pphase(twpx, "Y" )