Seismic Analysis of Earthquake Hypocenter determination

Description

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.

Details

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.

Note

Functions

CONTPF EQXYresid INITpickfile NLSlocate PFoutput RQ SavePF UPdateEQLOC XYSETUP Y2Pphase chak contPFarrivals doAmap gMAP getregionals prepPDE viewCHAC

Author(s)

Jonathan M. Lees<jonathan.lees.edu> Maintainer:Jonathan M. Lees<jonathan.lees.edu>

References

Lee, W.H.K., and S.W. Stewart, Principles and Applications of Microearthquake Networks, Academic Press, New York, 1981.

See Also

RSEIS

Examples

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

Description

1D Velocity Ecuador

Usage

data(ASW.vel)

Format

a list of velocities for hypocenter relocation

Source

Mario Ruiz

Examples

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"  ))

Jackknife earthquake location

Description

Perform jackknife on earthquake location by eliminating stations.

Usage

BLACKJACK(Ldat, vel)

Arguments

Details

Stations are eliminated, not rows.

Value

event list with pseudo values

Note

events are located with P and S-wave arrivals, but code here should eliminate just stations.

Author(s)

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

References

Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.

See Also

Vlocate, plotJACKLLZ

Examples

###### 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)

Button to Contour Pickfile Arrivals

Description

Button to Contour Pickfile Arrivals, used internally in swig.

Usage

CONTPF(nh, g, idev = 3)

Arguments

Details

Driver for contPFarrivals

Value

Side effects

Author(s)

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

See Also

contPFarrivals

Examples

if(interactive()){
######  interactive: addition of button in swig
data(GH, package='RSEIS')

buts = "CONTPF"
RSEIS::swig(GH, PADDLAB=buts, SHOWONLY=FALSE )
}

Distance wheighting

Description

Distance weighting for non-linear earthquake location.

Usage

DistWeight(dist, err, distwt)
DistWeightLL(lat, lon, elat, elon, err, distwt)
DistWeightXY(x, y, ex, ey, err, distwt)

Arguments

Details

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.

Value

vector of weights

Author(s)

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

Examples

DistWeight(1:10, .4, 20)

Calculate Residuals

Description

Given an earthquake hypocenter and a list of station information, retrieve the station residuals.

Usage

EQXYresid(XY, vel = list(), h1 = c(0, 0, 0, 0), PLOT = FALSE)

Arguments

Details

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.

Value

vector, right hand side of the least squares problem.

Author(s)

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

See Also

travel.time1D,UPdateEQLOC

Examples

####  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 Pand S travel times and derivatives

Description

Get Pand S travel times and derivatives

Usage

GETpsTT(phase, eqz = 6, staz = 0, delx = 1, dely = 1, deltadis = 6, vel)

Arguments

Details

Creates a vector of travel times, and a matrix and derivatives used for inversion.

Value

list:

Author(s)

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

See Also

many.time1D

Examples

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)

PICK Buttons for swig

Description

defining functions for swig

Usage

GPIX(nh, g)

Arguments

Details

Buttons can be defined on the fly.

GPIX

Multiple picks on a panel

Value

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.

Author(s)

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

See Also

swig, XTR

Examples

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)

}

First guess from a pick file

Description

Extract the lat lon from the pick file.

Usage

Gfirstguess(Ldat, type = "first")

Arguments

Details

Either the earliest arrival or the average station is returned. Used internally in the earthquake location program to provide a first guess.

Value

vector, lat, lon, z and tee

Author(s)

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

See Also

Klocate

Examples

data(GH, package='RSEIS')
WW = RSEIS::uwpfile2ypx(GH$pickfile)

twpx  = latlonz2wpx(WW, GH$pickfile$STAS )

g1 = Gfirstguess(twpx, type = "first")

Jackknife a list of events

Description

Jackknife a list of events

Usage

HiJACK(lps, sta, vel)

Arguments

Details

Driver for BLACKJACK

Value

jackknife pseudovalues for each event

Author(s)

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

References

Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.

See Also

BLACKJACK

Examples

#####  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 )

Initialize a pickfile

Description

Initialize a pickfile

Usage

INITpickfile(stas = NULL, src = NULL, WPX = NULL)

Arguments

Details

Initialize a pickfile with a set of picks extracted from swig.

Value

list, pickfile

Author(s)

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

See Also

EmptyPickfile

Examples

data(GH, package='RSEIS')
WW = RSEIS::uwpfile2ypx(GH$pickfile)

PF =  INITpickfile(stas=GH$stafile, src=NULL, WPX=WW )

Earthquake Hypocenter Location

Description

Earthquake Hypocenter Location

Usage

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))

Arguments

Details

Inversion is done with SVD.

Value

Event location in Lat-Lon-Z-T.

Note

Damped least squares.

Author(s)

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

See Also

swig, defaultVEL

Examples

######  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  )

  

List location data

Description

List location data

Usage

LDATlist(g1, w1)

Arguments

Value

side effects

Author(s)

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

Adjust times relative to least minute.

Description

Adjust times relative to least minute.

Usage

LeftjustTime(g1)

Arguments

Details

Reutrns the list with the times adjusted to the least minimum (left adjusted)

Value

list is returned.

Author(s)

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

See Also

recdate

Examples

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)

Mean Station Distance

Description

calculate the mean km distance of a set of Lat-lon pairs

Usage

MeanStaDist(Ldat)

Arguments

Details

Given a list with elements named lat and lon, find the mean station distance.

Value

scalar

Author(s)

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

See Also

setPROJ, GLOB.XY, dist

Examples

data(GH, package='RSEIS')
MeanStaDist(GH$pickfile$STAS)

Nonlinear Least Squares Location

Description

Nonlinear Least Squares Location using Gieger's method

Usage

NLSlocate(GH, vel = list(), init = c(0, 0, 0, 0), PLOT = FALSE)

Arguments

Details

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.

Value

vector, new location

Note

At this stage there are no weighting mechanisms or code to eliminate data that has residuals that are too large.

Author(s)

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

References

Lee, W.H.K., and S.W. Stewart, Principles and Applications of Microearthquake Networks, Academic Press, New York, 1981.

See Also

swig

Examples

data(GH, package='RSEIS')
###  location is:
eqsol = NLSlocate(GH, vel=GH$velfile,  PLOT=TRUE )

One Phase Pick Per Station

Description

Require only one pick per station of a specified phase.

Usage

OnePerSta(twpx, phase = "Y")

Arguments

Details

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.

Value

WPX list

Note

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.

Author(s)

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

See Also

cleanWPX, repairWPX

Examples

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

Description

Create a character string from a date for naming unique output files.

Usage

PCfiledatetime(orgtim, tims)

Arguments

Value

Author(s)

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

Examples

library(RSEIS)
data(GH, package='RSEIS')

g1 = getGHtime(GH)
g2 = unlist(g1)

PCfiledatetime(g2, 1)

Save WPX list

Description

Save a WPX list to a file on the local file system.

Usage

PCsaveWPX(twpx, destdir = NULL)

Arguments

Details

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.

Value

Side effects on file system. The name of the output file is returned.

Note

User must have write access to the destination directory.

Author(s)

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

See Also

RSEIS::addWPX, RSEIS::catWPX, RSEIS::checkWPX, RSEIS::cleanWPX, RSEIS::clusterWPX, RSEIS::repairWPX, RSEIS::setWPX

Examples

#####  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

Description

Write a pickfile to disk, after updating the earthquake location, in a variety of formats.

Usage

PFoutput(PF, stas = NULL, sol = NULL, format = 0, destdir=NULL)

Arguments

Details

Writes files to disk in local directory.

Value

Side effects: writes files to user's disk

Note

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.

Author(s)

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

See Also

SavePF, RSEIS

Examples

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() )

PICK Buttons for swig

Description

Picking functions for swig

Usage

Pick3(nh, g)

Arguments

Details

Buttons can be defined on the fly.

Pick3

Multiple picks on a panel

Value

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.

Author(s)

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

See Also

swig, PickWin

Examples

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)

}

Post Processing on EQrquake

Description

Post Processing on EQrquake

Usage

PostREQquake(XQ, proj)

Arguments

Details

Following event locations, plot.

Value

graphical side effects

Author(s)

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

Plotting error ellipsoids of many events

Description

Plotting error ellipsoids of many events

Usage

PostVquake(MANYeq, GX, GY, XY, proj, add=FALSE, ...)

Arguments

Details

Plots the event and the error ellipsoids

Value

Graphical side effects

Note

This is used to plot many event locations and their error ellipsoids

Author(s)

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

See Also

eqlipse

Range of Date Time

Description

Return the range of dates and times for any list with a date/time list

Usage

Qrangedatetime(D)

Arguments

Value

Author(s)

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

Examples

library(RSEIS)
data(GH, package='RSEIS')

v = Qrangedatetime(GH$info)

Rquake Button

Description

Driver for NLSlocate

Usage

RQ(nh, g, idev = 3)

Arguments

Details

Button to be called from within swig after picking.

Value

new hypocenter

Author(s)

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

See Also

NLSlocate, EQXYresid, XYSETUP, swig,chak

Examples

if(interactive()){
##### interactive  
data(GH, package='RSEIS')

 buts = c("GPIX","PPIX", "PickWin",
         "fspread", "gMAP", "RQ" , "CONTPF")

RSEIS::swig(GH, PADDLAB=buts)
}

Button to reset the choices of station and component

Description

Button to reset the choices of station and component in swig and Mine.seis

Usage

ReSet(nh, g)

Arguments

Details

Driver for SELstaDB

Value

Side effects

Author(s)

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

See Also

SELstaDB, Mine.seis

Examples

if(interactive()){
data(GH, package='RSEIS')

buts = "ReSet"
RSEIS::swig(GH, PADDLAB=buts)

}

Rows to Keep for inversion

Description

Selects which rows in the hypocenter determination to keep during non-linear itaerations based on robust rsidual elimination.

Usage

Rowz2Keep(Ldat, EQ, G1, RESMAX)

Arguments

Details

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.

Value

Index of good rows

Author(s)

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

See Also

XYlocate

Pick stations and components interactively

Description

Pick stations and components interactively. This is a routine used in swig.

Usage

SELstaDB(IDB, sel=1, newdev=TRUE, STAY=FALSE)

Arguments

Value

vector of index to list of stations and components

Author(s)

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

See Also

infoDB, makeDB

Examples

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)

}

Save Pick File Button

Description

Save a pick file from within swig

Usage

SavePF(nh, g)

Arguments

Details

Uses PFoutput to save a pickfile to disk.

Value

Side Effects

Note

Pickfile is saved as a native R file with wpx extension

Author(s)

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

See Also

PFoutput

Examples

if(interactive()){
data(GH, package='RSEIS')
buts = "SavePF"
RSEIS::swig(GH, PADDLAB=buts)

}

Update an Earthquake location

Description

Update an Earthquake location following a relocation.

Usage

UPdateEQLOC(PF, sol, vel, stas = NULL)

Arguments

Details

After re-picking or changing the model or the station corrections, update the event location in the pickfile.

Value

Pickfile List

Author(s)

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

See Also

EQXYresid, NLSlocate,PFoutput

Examples

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

Description

Hypocenter Determination with error checking and adjustments.

Usage

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)

Arguments

Details

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.

Value

list:

Note

The schedule may be adjusted by duplicating this function and changing the maxit parameters.

Author(s)

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

References

Lee and Stewart

See Also

XYlocate, Klocate, DoRLocate

Examples

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)

Set up matrix for hypocenter inversion

Description

Set up matrix for hypocenter inversion

Usage

XYSETUP(STAS, init, vel)

Arguments

Details

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.

Value

matrix

Note

Need scheme for weighting according to errors in picks and distance weighting.

Author(s)

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

See Also

setPROJ, GLOB.XY,NLSlocate

Examples

##  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  )

Error Bars in X and Y

Description

Error Bars in X and Y

Usage

XYerror.bars(x, y, xlo = 0, xhi = 0, ylo = 0,
yhi = 0, pch = 1, col = 1, barw = 0.1, add = FALSE, ...)

Arguments

Value

graphical side effects

Author(s)

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

Examples

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)

Locate Earthquake with UTM projection

Description

Non-linear hypocenter location with UTM geographical projection. Used for locating earthquakes in local or regional settings.

Usage

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)

Arguments

Details

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.

Value

List:

Note

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.

Author(s)

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

See Also

Vlocate

Examples

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 = 2

Convert Y-phase to P-phase

Description

Removes extraneous other-phase from a pick file. If Ypix were made initially as a rough pick, this removes them.

Usage

Y2Pphase(twpx, phase)

Arguments

Details

Initially many events may be picked using GPIX button. These should be removed after the P-phases have been determined with PickWin.

Value

WPX returned without other-phases

Author(s)

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

See Also

PPIX, GPIX, YPIX, PickWin

Examples

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" )

Check Location data

Description

Check to see if location data has the minimally correct list components.

Usage

checkLOCATEinput(Ldat, EQ, vel = NULL)

Arguments

Details

Input pick list must have at x,y,z, sec, cor, err elements for each station.

Value

logical: FALSE mean problem with data

Author(s)

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

See Also

XYlocate

Examples

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)

Cluster Analysis of Picks

Description

Given a pick file in WPX format, break the picks apart clustered accoring to single link cluster analysis.

Usage

clusterWPX(twpx, tol = 200, PLOT = FALSE)

Arguments

Details

If there is not significant separation of picks, only one cluster is returned. To avoid spurious clusters, increase the tolerance.

Value

list of WPX lists

Note

Cluster depends on what one considers a cluster.

Author(s)

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

See Also

RSEIS::addWPX, RSEIS::catWPX, RSEIS::checkWPX,RSEIS::cleanWPX, PCsaveWPX, RSEIS::setWPX, RSEIS::repairWPX

Examples

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)

Contour Pickfile Arrivals

Description

Contour plot of arrival times recorded in a pickfile list.

Usage

contPFarrivals(PF, stas, proj=NULL, cont=TRUE, POINTS=TRUE, image=FALSE ,
             col=RSEIS::tomo.colors(50), gcol="black",   phase="P", add=TRUE)

Arguments

Details

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.

Value

Graphical Side Effects

Author(s)

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

See Also

doAmap

Examples

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 File

Description

Coso Station Location file, 1989-1999

Usage

data(coso_sta_LLZ)

Format

Name, Lat, Lon, Z

Source

Personal Files

References

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.

Selection of pickfiles from Coso Geothermal Field

Description

Set of selected seismic arrival files with hypocenter locations.

Usage

data("cosopix")

Format

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

Details

Each element of this list is an individual earthquake record.

Examples

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

Description

Default Velocity Function is returned in the event no velocity function is available.

Usage

defaultVEL(kind = 1)

Arguments

Details

A set of default velocity functions are available.

Value

velocity list, P and S waves

Author(s)

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

See Also

fuj1.vel

Examples

v = defaultVEL(1)

Plot a map of station locations

Description

Plot a map of station locations

Usage

doAmap(stas, doproj = TRUE)

Arguments

Details

The range of the plot is expanded by 10 percent prior to plotting.

Value

list, GEOmap projection

Author(s)

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

See Also

gMAP,expandbound,GLOB.XY

Examples

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

Description

Error Elipse for Hypocenter Location

Usage

eqlipse(x, y, cov, wcols = c(1, 2), dof = 2, pct=0.05, ...)

Arguments

Details

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.

Value

Side effects, graphical

Author(s)

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

See Also

eqwrapup

Examples

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"  )

Earthquake Wrap Uo

Description

Calculate error and summary information on earthquake location.

Usage

eqwrapup(Ldat, EQ, vel, distwt=20, lambdareg = 0.0, verbose=FALSE)

Arguments

Details

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.

Value

List

Note

The Damping parameter (lambda) is set to zero. In the UW lquake program, lambda is set to 0.02.

Author(s)

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

See Also

Klocate, Glocate, getGAP

Examples

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)

Euler Rotation Angles

Description

Given three angles return rotation matrix.

Usage

euler_passive(phi, theta, psi)

Arguments

Details

Code borrowed from cpp code in package cda. used in rgl.ellipsoid.

Value

3 by 3 rotation matrix.

Author(s)

Jonathan M. Lees<jonathan.lees@unc.edu>, Baptiste Auguie<baptiste.auguie@gmail.com>

See Also

rgl.ellipsoid

Examples

options(rgl.useNULL = TRUE)
phi=30*pi/180 ; theta= 20*pi/180; psi = 6*pi/180
rr = euler_passive(phi,theta,psi)

Generic Map Button

Description

Generic Map Button

Usage

gMAP(nh, g, idev = 3)

Arguments

Details

This is a button used internally in swig

Value

Graphical Side Effects

Author(s)

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

See Also

swig

Examples

if(interactive()){
####  this is interactive
###  adds button to swig menu
data(GH, package='RSEIS')

buts = "gMAP"
RSEIS::swig(GH, PADDLAB = buts )

}

Get Eulers Angles

Description

Given a covariance matrix calculated with Vlocate, extract euler's angles for plotting in rgl

Usage

getEulers(R)

Arguments

Details

Extract the euler angles for plotting an ellipsoid. psi about X-axis, theta about Y axis, phi about Z-axis.

Value

vector, phi theta psi

Note

Used in conjunction with ROTcovQUAKE

Author(s)

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

See Also

ROTcovQUAKE

Examples

options(rgl.useNULL = TRUE)
R = matrix( runif(9), ncol=3)

getEulers(R)

Get Seismic Gap

Description

Given an earthquake and a set of stations, return the maximum angle subtended between adjacent stations relative to the epicenter.

Usage

getGAP(EQ, Ldat, PLOT = FALSE)

Arguments

Details

Theangles are calculated in cartesian coordinates with the epicenter at the origin using a UTM projection.

Value

numeric, gap in degrees

Author(s)

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

See Also

eqwrapup

Examples

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)

}

Extract regional events

Description

Extract regional events from a hypocenter list (catalog)

Usage

getregionals(KAT, Mlat, Mlon, rad = 1000, t1 = 1, t2 = 2)

Arguments

Details

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.

Value

Catalog

Author(s)

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

See Also

RSEIS::Mine.seis, RSEIS::swig

Examples

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 )

Travel time residuals

Description

Given an earthquake location and a set of arrival times, return a vector of residuals.

Usage

getresidTT(Ldat, EQ, stas, vel)

Arguments

Details

1D travel time calculation.

Value

vector of residuals

Author(s)

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

See Also

travel.time1D

Examples

#########  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)

Image Influence of stations

Description

Plot contours/image of Influence scores.

Usage

imageINFLUENCE(B, sta, proj)

Arguments

Details

Following jackknife - plot results. this function is called by plotJACKLLZ.

Value

side effects

Author(s)

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

References

Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.

See Also

plotJACKLLZ

Last Pix

Description

'RSEIS' Button: Restore Last WPX file from memory. Function is used internally in swig.

Usage

lastPIX(nh, g)
editPIX(nh, g)

Arguments

Value

New WPX list attached to g

Author(s)

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

Add Lat-Lon-Z to WPX list

Description

Given an existing list of seismic picks, add Latitude, Longitude and Elevation associated with the indicated station.

Usage

latlonz2wpx(twpx, stas)

Arguments

Details

The names of the stations are matched to the station names int he station file.

Value

Pick file with LLZ added as list members.

Author(s)

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

See Also

Klocate

Examples

data(GH, package='RSEIS')
WW = RSEIS::uwpfile2ypx(GH$pickfile)

twpx  = latlonz2wpx(WW, GH$pickfile$STAS )

Legitimate Pix

Description

Check WPX list for legitimate picks

Usage

legitWPX(twpx, quiet=TRUE)

Arguments

Details

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.

Value

integer: 0=not legitimate, 1=legitimate

Note

Currently only the name is tested for all(NA), but this might be changed int he future for a more sophisticated test.

Author(s)

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

See Also

PCsaveWPX

Examples

###  test fails

library(RSEIS)
jk = RSEIS::cleanWPX()
legitWPX(jk)

####  test passes:
data(GH, package='RSEIS')
gwpx = RSEIS::uwpfile2ypx(GH$pickfile)

legitWPX(gwpx)

Plot Earthquake location

Description

Plot Earthquake location

Usage

plotEQ(Ldat, AQ, add = FALSE, prep = FALSE,
TIT = "UTM Projected Stations", proj = NULL,
 xlim = NULL, ylim = NULL)

Arguments

Details

used internally in RElocateEQ

Value

graphical side effects

Author(s)

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

See Also

RElocateEQ

BoxPlot Jackknife of station locations

Description

BoxPlot Jackknife of station locations

Usage

plotJACKLLZ(hjack, sta, proj = NULL,  PLOT=1)

Arguments

Details

takes the output of the HiJack function and extracts the pseudovalues and influence information for boxplots.

Value

Graphical side effects and

Author(s)

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

References

Iversen, E. S., and J. M. Lees (1996), A statistical technique for validating velocity models, Bull. Seismol. Soc. Am. 86(6), 1853-1862.

See Also

HiJACK, BLACKJACK,imageINFLUENCE, Vlocate

Examples

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 )

Rip off Event location information

Description

Extract Event location information following Vlocate

Usage

ripper(AQ)

Arguments

Details

Extract lat-lon from event locations to track intermediate solutions and convergence

Value

2 by N matrix, lat-lon

Author(s)

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

See Also

plotEQ

Examples

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)