Rquake (version 2.5-1)

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

Value

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

Arguments

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

Author

Jonathan M. Lees<[e-mail]>

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.

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