Seismic Analysis and Display

Description

Multiple interactive codes to view and analyze seismic data, via spectrum analysis, wavelet transforms, particle motion, hodograms. Includes general time-series tools, plotting, filtering, interactive display.

Note

Seismic Sections

JGET.seis view.seis swig Mine.seis VELOCITY.SEISN DISPLACE.SEISN ZOOM.SEISN wlet.drive SENSORsensitivity PLOT.MATN PLOT.SEISN PLOT.TTCURVE PLOT.ALLPX plotevol MTMdisp MTMplot NEW.getUWSTAS NEWPLOT.WPX INSTFREQS INSTresponse GLUE.GET.seis GLUEseisMAT FILT.SEISN FILT.spread CHOP.SEISN get.corner grotseis

Velocity-Travel Time:

Put1Dvel Ray.time1D setLQUAKE selAPX Get1Dvel Comp1Dvel Comp1Dvels travel.time1D

Particle Motion:

hodogram PMOT.drive complex.hodo addpoints.hodo idpoints.hodo DO.PMOT.ARR partmotnet prep1wig prepSEIS EmptySEIS GAZI

Time series:

xcor2 wlet.drive wlet.do wiggle.env plotwlet STLTcurve SPECT.drive rsspec.taper evolfft GETARAIC PSTLTcurve getphaselag2 envelope hilbert LocalUnwrap lagplot applytaper autoreg butfilt choosfilt MTM.drive

Date-Time Functions:

yeardate YRsecdif Zdate recdatel recdate tojul getjul getmoday secdifL secdif secdifv JtimL Jtim fromjul

Graphics:

plocator ilocator meshgrid ymarginfo zlocator winmark vline screens RESCALE pwlet2freqs addtix circle circ letter.it jpostscript JBLACK JGRAY HOZscale gaddtix Gcols jlegend tomo.colors

Misc:

BKpfile2ypx brune.doom brune.func brune.search comp.env contwlet deconinst detail.pick rdistaz rDUMPLOC EmptyPickfile ETECTG finteg fixcompname fixcomps fixUWstasLL fmod FRWDft getb1b2 getNcard getpfile getseisinfo getvertsorder gpoly GreatDist gwpix2ypx hilow hypot integ1 INVRft itoxyz jadjust.length jpolyval jstats local.file logspace makefreq mirror.matrix Mmorlet mtapspec peaks PICK.DOC pickit plt.MTM0 PLTpicks PPIX PreSet.Instr ReadSet.Instr readUW.OSTAS scal2freqs SEARCHPIX setstas setwelch shade.col SNET.drive T12.pix Thresh.J TOCART trapz tung.pulse unpackAcard uwpfile2ypx

Author(s)

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

See Also

RPGM, RFOC

Examples

data("GH")
swig(GH)

ASCII RSEIS data dump

Description

Write RSEIS list to a file in ASCII format.

Usage

ASCII.SEISN(GH, sel = 1, HEAD = TRUE, destdir='.' )

Arguments

Details

Used for data exchange for users who do not want to use RSEIS. The header consists of one line start date (yr, jd, hr, min, sec) and sample rate (dt).

Value

Side effects - files are created.

Author(s)

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

Examples

## Not run: 
##### this example creates an ascii version of the
#####   seismic data for exchange purposes
data("GH")
tempd = tempdir()
sel <-  which(GH$COMPS == "V" & GH$STNS=="CE1"   )
ASCII.SEISN(GH, sel = 1, HEAD = TRUE, destdir=tempd)

## End(Not run)

Butterworth filter with Augmentation

Description

Design and apply butterworth low/high/band pass filters with augmentation of the signal on either end to suppress edge effects.

Usage

AUGMENTbutfilt(a, fl = 0, fh = 0.5, deltat = 1, type = "BP",
proto = "BU", npoles = 5, chebstop = 30, trbndw = 0.3,
RM = FALSE, zp = TRUE, pct = 0.1)

Arguments

Details

Creation of butfilt is a described by the following arguments:

LP

low pass

HP

high pass

BP

band pass

BR

band reject

BU

Butterworth

BE

Bessel

C1

Chebyshev type 1

C2

Chebyshev type 2

Arguments chebstop , trbndw are ignored for non-chebyshev filters. LP and HP filters are seet by specifying fl for HP filters and fh for LP filters, the other argumentin each case is ignored.

Mean values should be removed prior to calling this function, and then set RM=FALSE. This is true especially if tapering is applied prior to filtering.

Zero phase filter is achived by running filter back and forth. Otherwise a single pass is returned. This should be equivalent to package signal filtfilt (from MATLAB).

Augmentation involves copying the first and last percent of the signal, reversiing the time and adding to the signal on each end. This is then filtered, and removed after filter is complete. It is assumed that the important part of the signal is in the center of the time series and the edges are less critical. Then the augmented part has the same statistical content as the edges of the signal (presumably noise) and will not affect the filtered signal considerably. This is then thrown away prior to return.

Value

Filtered time series with the augmentation removed after filter.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

butfilt

Examples

data(CE1)

ts1  <-  CE1$y
zz  <-  AUGMENTbutfilt(ts1, fl=1, fh=15,  deltat=CE1$dt, type="LP" ,  proto="BU",
npoles=5 )


##############    second example with plotting


data(KH, package ='RSEIS' )
w = KH$JSTR[[1]]
dt = KH$dt[1]

x = seq(from=0, by=dt, length=length(w));
plot(x,w, type='l')

par(mfrow=c(2,1) )
        
i=1
       fl = 1/50
fh= 1/2
       ftype = 'BP'
       ##########  normal band pass filter
        
 zz = butfilt(w, fl, fh,  dt, ftype ,  "BU")
     f.stamp =   filterstamp(fl=fl, fh=fh, type=ftype)

plot(x, zz, type='l', xlab='s', ylab='amp', main= f.stamp)
title(sub='butfilt')
       
       ####  
   zz1 = AUGMENTbutfilt(w, fl, fh,  dt,    type=ftype ,  proto="BU", zp=TRUE, pct=0.2 )
     f.stamp =   filterstamp(fl=fl, fh=fh, type=ftype)
plot(x, zz1, type='l', xlab='s', ylab='amp', main= f.stamp)
title(sub='AUGMENTbutfilt')

Replot Function for SELBUT

Description

Replot Function for SELBUT

Usage

BUTREPLOT(opts, ncol = 5, HOZ = TRUE, TOP = TRUE,
 cols = "white", main = "", xlim = c(0, 1),
        ylim = c(0, 1), newplot = TRUE)

Arguments

Details

Used internally in SELBUT as a replotting function

Value

list

Author(s)

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

See Also

SELBUT, swig

Examples

STDLAB <- c("DONE", "QUIT", "zoom.out", "zoom.in",
  "SELBUT", "FILT","UNFILT", "PSEL", "SGRAM",
       "WLET", "SPEC", "XTR" )
BUTREPLOT(STDLAB)

Single Seismogram

Description

Single Seismogram from Coso California

Usage

data(CE1)

Format

list(x=0, y=0, dt=0, name ="", Tpick=0, mark ="", deltat=0)

References

Lees, J.M., 2004. Scattering from a fault interface in the Coso geothermal field. Journal of Volcanology and Geothermal Research, 130(1-2): 61-75.

Examples

data(CE1)
plot(CE1$x, CE1$y, type='l')

CHOP SEISmic structure

Description

Take a seismic structure and return a time limited version

Usage

CHOP.SEISN(GH, sel = 1:4, WIN = NULL)

Arguments

Value

Seismic trace structure

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

data("GH")
sel <-  which(GH$COMPS=="V")

KF  <-  CHOP.SEISN(GH, sel=sel, WIN = c(0 , 5)  )

swig(KF, SHOWONLY=0)

Seismic Component Order

Description

Set seismic component order

Usage

COMPorder(STNS, COMPS)

Arguments

Details

Sets up components so they are ordered according to V, N, E. used internally in swig.

Value

order vector

Author(s)

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

Compare a pair of 1D models

Description

plot a pair of 1D velocity Models for comparison

Usage

Comp1Dvel(v, v2, col=c('blue', 'brown'), ...)

Arguments

Value

Graphical Side Effects

Author(s)

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

See Also

Plot1Dvel, Get1Dvel, travel.time1D

Examples

VEL <- list()
VEL$'zp' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vp' <- c(1.1,2.15,3.2,4.25,5.3,6.25,6.7,6.9,7,7.2)
VEL$'ep' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'zs' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vs' <- c(0.62,1.21,1.8,2.39,2.98,3.51,3.76,3.88,3.93,4.04)
VEL$'es' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'name'  <-  '/data/wadati/lees/Site/Hengil/krafla.vel'

VELNish <-  list()
VELNish$'zp' <- c(0,0.1,0.6,1.1,21.1)
VELNish$'vp' <- c(2.8,3.4,4.1,4.7,4.7)
VELNish$'ep' <- c(0,0,0,0,0)
VELNish$'zs' <- c(0,0.1,0.6,1.1,21.1)
VELNish$'vs' <- c(1.6,2,2.4,2.7,2.7)
VELNish$'es' <- c(0,0,0,0,0)
VELNish$'name' <- 'Nish'

Comp1Dvel(VEL, VELNish)

Compare 1D models

Description

Plot 1D velocity Models for comparison.

Usage

Comp1Dvels(INV, depth = 1:50)

Arguments

Details

takes several velocity models, finds the range of all, makes a plot so that all models fit on figure.

Value

Graphical Side Effects

Author(s)

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

See Also

Plot1Dvel, Comp1Dvel, Get1Dvel

Examples

VEL <- list()
VEL$'zp' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vp' <- c(1.1,2.15,3.2,4.25,5.3,6.25,6.7,6.9,7,7.2)
VEL$'ep' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'zs' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vs' <- c(0.62,1.21,1.8,2.39,2.98,3.51,3.76,3.88,3.93,4.04)
VEL$'es' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'name' <- '/data/wadati/lees/Site/Hengil/krafla.vel'

VELNish <- list()
VELNish$'zp' <- c(0,0.1,0.6,1.1,21.1)
VELNish$'vp' <- c(2.8,3.4,4.1,4.7,4.7)
VELNish$'ep' <- c(0,0,0,0,0)
VELNish$'zs' <- c(0,0.1,0.6,1.1,21.1)
VELNish$'vs' <- c(1.6,2,2.4,2.7,2.7)
VELNish$'es' <- c(0,0,0,0,0)
VELNish$'name' <- 'Nish'

Comp1Dvels(c("VEL", "VELNish"))

Days per Year

Description

Calculate the number of days per calendar year

Usage

DAYSperYEAR(yr)

Arguments

Value

Author(s)

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

Examples

years <- seq(from=1850, to=2010, by=1)

DAYSperYEAR(years)

Decimate a set of traces

Description

Decimate, or reduce the sample rate of a set of traces stored in event RSEIS format

Usage

DECIMATE.SEISN(TH, sel=1:length(TH$JSTR), dec=5 ,
 type="LP", proto="BU" , fl=2,  fh=10, RM=FALSE, zp=TRUE )

Arguments

Details

Reduces the number of samples by skipping every "dec" sample.

To achieve smoothing prior to sampling, low pass filter may be applied to avoid spikes or other sampling issues.

If type is FALSE, no filter is applied and samples are taken from the input.

Value

an RSEIS list.

Note

The dt, n and t2 are modified in info.

Author(s)

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

See Also

butfilt, downsample

Examples

data(GH)

dec = 250/50

#####  resample all traces by reducing from 250 to 50 samples/s
DH = DECIMATE.SEISN(GH, sel=1:length(GH$JSTR), dec=dec ,
 type="LP", proto="BU" , fl=2,  fh=50, RM=FALSE, zp=TRUE )

#####  compare

#####   times in
### starting second should be the same
GH$info$sec[1:5]  - DH$info$sec[1:5] 
####  number of samples should be reduced
cbind(GH$info$n[1:5] , DH$info$n[1:5]  )
### ending seconds should be close but not identical
cbind(GH$info$t2[1:5] , DH$info$t2[1:5]  )

cbind(GH$info$dt[1:5] , DH$info$dt[1:5]  )

cbind( sapply(GH$JSTR, 'length'), sapply(DH$JSTR, 'length') )



####  for visual comparison:
### par(mfrow=c(2,1) )
## g =  swig(GH, sel=which(GH$COMPS=="V" ), SHOWONLY=0 )
## d =  swig(DH, sel=which(DH$COMPS=="V" ), SHOWONLY=0 )

  

Displacement seismogram

Description

Removes seismic instrument response and integrates to displacement.

Usage

DISPLACE.SEISN(TH, sel = 1:length(TH$JSTR), inst = 1,
Kal = Kal,waterlevel = 1e-08,  FILT = list(ON = FALSE,
 fl = 1/30, fh = 7, type = "HP", proto = "BU",RM=FALSE, zp=TRUE))

Arguments

Details

Instrument responses are lists of poles and zeros for each instrument defined.

Value

Same as input list with new traces representing displacement versus velocity

Author(s)

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

See Also

VELOCITY.SEISN, deconinst, butfilt

Examples

data(KH)

Kal <- PreSet.Instr()

DH  <-   DISPLACE.SEISN(KH, sel = 1 , inst = 1,
Kal = Kal, FILT = list(ON = FALSE, fl = 1/200, fh = 7,
type = "BP", proto = "BU"))

if(interactive()){
  SOUT  <-   swig(DH, PADDLAB=c("CENTER", "fspread", "HALF", "PREV")  )
}

Distance Cross section

Description

Plot time series vertically at specified distances. Produces a seismic cross section with correct spacing between traces.

Usage

DISTxsec(GH, dist, TIM.WIN = c(0, 3600), sel, trace.width = 10,
col = "black", text.col = "blue", text.font = 2, text.size = 0.8,
 add = FALSE, plot = TRUE)

Arguments

Details

Distances should be a vector for each trace in the RSEIS list.

Value

vector of x-y coordinates of the plot.

Author(s)

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

See Also

swig, prepSEIS

Examples

#### example using data in the RSEIS package
data(GH)
####  get the source location
lat.org  = GH$pickfile$LOC$lat
lon.org  = GH$pickfile$LOC$lon
####  get the station locations
g1 =GH$stafile
#### find the distance to each station
gd = rdistaz(lat.org, lon.org, g1$lat, g1$lon )

##### optional, filter the data
sel= which( GH$COMPS == 'V')
### filter traces
Fdef  <-  list(ON=TRUE, fl=1, fh=1, type="HP", proto="BU", RM=TRUE, zp=TRUE )
KF  <-  FILT.SEISN(GH, FILT=Fdef)

### match the stations in GH to the station distances
m1 = match(GH$STNS  , g1$name) 
dist.GH = gd$dist[m1]
TIM.WIN = range(GH$ex) 


#######  prepare plot, but do not add traces
A = DISTxsec(KF, dist.GH, TIM.WIN, sel, trace.width = 0.5 , add=FALSE,
plot=FALSE )
#####  add traces
B = DISTxsec(KF, dist.GH, TIM.WIN, sel, trace.width = 0.5 , add=TRUE,
plot=TRUE, col='black'  , text.col='red', text.size=1  )

Particle Motion Analysis with arrows

Description

Plot particle motion arrows

Usage

DO.PMOT.ARR(E, N)

Arguments

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

PMOT.drive

Examples

data(GH)
XLIM = c(1226, 1322 )

e = GH$JSTR[[1]][XLIM[1]:XLIM[2]]
n = GH$JSTR[[2]][XLIM[1]:XLIM[2]]


   xx = range(e, na.rm =TRUE)
    yy = range(n, na.rm =TRUE)
 sx = range(c(xx, yy))

   x  = RPMG::RESCALE(e, 0, 1, sx[1], sx[2])
    y  = RPMG::RESCALE(n, 0, 1, sx[1], sx[2])



plot(range(x), range(y) , type='n')
lines(x, y, col=grey(0.8) )
DO.PMOT.ARR(x, y)

Gabor transform

Description

Gabor transform with simple spectrum

Usage

DOsgram(Xamp, DT = 0.008, multi = 1, scale.def = 0,
 TWIN = 2, TSKIP = 0.2, PCTTAP = 0.05, PLOT=TRUE)

Arguments

Details

This is a non-interactive version of SPECT.drive.

Value

list

Author(s)

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

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

See Also

evolMTM, evolfft, evolAR, plotevol

Examples

data(KH)
###   swig(KH)

Xamp  <-  KH$JSTR[[1]]
Xamp  <-  Xamp[57914:72989]

Nfft <- 1024   ###  fft length
Ns <- 512      ###  number of samples in a window
Nov <- 480    ###  number of samples of overlap per window
fl <- 0        ###  low frequency to return
fh <- 12     ###  high frequency to return


EV  <-  DOsgram(Xamp, DT = 0.008, multi = 1, scale.def = 0,
TWIN = 2, TSKIP = 0.2, PCTTAP = 0.05)

Epoch Day

Description

Number of days since Origin Year

Usage

EPOCHday(yr, jd = 1, origyr = 1972)

Arguments

Details

Either jd or mo, dom can be provided

Value

List:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

EPOCHyear, recdate

Examples

tyears   <-  1973:2009

E1  <-  EPOCHday(tyears, jd=1,   origyr=1972       )

EPOCHyear(E1$jday,   origyr=1972  )

Epoch Year

Description

Get year and julian day given number of days since origin

Usage

EPOCHyear(iday, origyr = 1972)

Arguments

Value

List:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

EPOCHday, recdate

Examples

tyears   <-  1973:2009

E1  <-  EPOCHday(tyears, jd=1,   origyr=1972       )

EPOCHyear(E1$jday,   origyr=1972  )

#######  here is an example using year Month and day of month
###  use March 19 for each year:
ii  <-  tojul(tyears, 3, 19)-tojul(tyears, 1, 1)

E1  <-  EPOCHday(tyears, jd=ii,  origyr=1972       )

EPOCHyear(E1$jday,   origyr=1972  )

Event Detection

Description

Event Detection for a seismic section

Usage

ETECTG(GH, sel = sel, FRWD = 8, BKWD = 8, sbef = 1,
saft = 6, DFRWD = 0.5, DBKWD = 0.5, thresh = 2,
Tthresh2 = 7, stretch = 1000, flo = 0.1, fhi = 5,
 PLOT = FALSE, Kmin = 7, perc = 0.05, kind = 1, DOARAIC = FALSE)

Arguments

Details

Very complicated picking routine - designed for volcanic regions with emergent arrivals. Works with lots of tuning.

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Create an empty RSEIS pickfile structure

Description

Creates a structure list with no data

Usage

EmptyPickfile(GH)

Arguments

Value

RSEIS pickfile list

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

EmptySEIS

Examples

data(GH)
EmptyPickfile(GH)

Create an empty RSEIS structure

Description

Creates a structure list with no data

Usage

EmptySEIS()

Value

RSEIS list

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

EmptyPickfile

Examples

EmptySEIS()

Fake Data for Examples.

Description

Create a list of artifical seismic traces to illustrate examples that require a database or long sequences.

Usage

FAKEDATA(amp, OLDdt = 0.01, newdt = 0.1, yr = 2000,
JD = 5, mi = 0, sec = 0, Ntraces = 48, seed = 200,
noise.est = c(1, 100), verbose = FALSE)

Arguments

Details

The input signal can be any time series, or even a made up signal. This is just to give the look of the result something like real data. The noise level is extracted from the man and std of the real data at the samples indicated by noise.est.

The sampling rate (dt, sec/sample ) is increased mainly for speed and plotting. This may be skipped for certain functions involving spectrum analysis.

The signal is distributed randomly in each hour along the total span of the requested period, i.e. each hour has one instance of the signal.

The date is arbitrary, of course.

Value

List of data in a format similar to the output of GET.seis.

Author(s)

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

See Also

GET.seis

Examples

#####  get a time series
data(KH)

amp = KH$JSTR[[1]]
OLDdt = KH$dt[1]
####  downsample to:
newdt = 0.1

JK = FAKEDATA(amp, OLDdt=OLDdt, newdt = 0.1, yr = 2000,
         JD = 4, mi = 12, sec = 0,  Ntraces = 3,
seed=200, noise.est=c(1, 100) , verbose=TRUE  )

 op <- par(no.readonly = TRUE)
par(mfrow=c(length(JK), 1) )
for(i in 1:length(JK) )
{
    DATTIM = paste(c(unlist(JK[[i]]$DATTIM), JK[[i]]$N), collapse=' ')
    
    plotGH( JK[[i]] )
    mtext(DATTIM, side=3, at=JK[[i]]$DATTIM$t2/2)
    
}
 par(op)

Filter Traces

Description

Filter Traces in a seismic structure

Usage

FILT.SEISN(TH, sel = 1:length(TH$JSTR),
FILT = list(ON = TRUE, fl = 0.5, fh = 7, type = "HP",
 proto = "BU", RM=FALSE, zp=TRUE ), TAPER = 0.1, POSTTAPER = 0.1, AUGMENT=FALSE)

Arguments

Details

RSEIS Seismic structure is filtered, trace by trace. If AUGMENT is TRUE, traces are augmented at beginning and end, filtered and then truncated to suppress edge effects. In that case no tapering is applied post fitler.

Value

RSEIS Seismic structure, traces are filtered and a proc is added to the trace history.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

butfilt

Examples

##   Fdef = choosfilt()
Fdef  <-  list(ON=FALSE, fl=0.5, fh=7.0, type="BP", proto="BU", RM=FALSE, zp=TRUE )
data("GH")
sel <-  which(GH$COMPS=="V")

sel  <-  1:3
 KF  <-  FILT.SEISN(GH, sel = sel, FILT=Fdef)
swig(KF, sel=sel, SHOWONLY=0)

Filter trace with a spread of filters

Description

Show a time series and a spread of user defined filters to show signal at a variety of bandwidths.

Usage

FILT.spread(x, y, dt, fl = fl, fh = fh, sfact = 1,
 WIN = NULL, PLOT = TRUE, TIT = NULL, TAPER = 0.05,
 POSTTAPER=0.05, RM=FALSE, zp=TRUE )

Arguments

Details

Use the TAPER and POSTTAPER to reduce the edge effects prior to and after filtering.

Value

list:

Author(s)

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

See Also

butfilt, PLOT.MATN

Examples

data(KH)
dt  <-  KH$dt[1]

y  <-   KH$JSTR[[1]]

x  <-   seq(from=0, by=dt, length=length(y))



fl <- rep(1/100, 5)
fh <- 1/c(1,2,5,10,20)

FILT.spread(x, y, dt, fl = fl, fh = fh, sfact = 1,
   WIN = NULL, PLOT = TRUE, TIT = NULL, TAPER = 0.05)

Forward fourier Transform

Description

Forward fourier Transform

Usage

FRWDft(g, n, tstart, dt)

Arguments

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

INVRft

Examples

zil  <-  rnorm(300)
fss  <-  FRWDft( zil, length(zil), 0, 0.004)   

Get azimuthal particle motion

Description

Do particle motion analysis

Usage

GAZI(ADAT, dt = 1, ex = seq(0, 100),
comp = c(4, 5, 6), sta = "ZZZ", az = 0,
len = 50, shift = 10, prev = 1, fileid = "", picks = NA, labs = NA)

Arguments

Value

list(aex=aex[1:jall], rateig=rateig[1:jall], aaz=aaz[1:jall], ai=ai[1:jall], figaz=figaz, azpar=azpar, incpar=incpar )

Examples

data("GH")


temp  <-  cbind(GH$JSTR[[4]], GH$JSTR[[5]], GH$JSTR[[6]])

pmolabs <- c("Vertical", "North", "East")


 G   <-  GAZI(temp, dt =GH$dt[4]  , comp = pmolabs, sta = GH$STNS[4] ,
az = 0, len =75, shift = 10, prev = 1)

Reads various seismic file formats

Description

This fuction calls binary routines to read in ‘segy’, ‘sac’.

Usage

GET.seis(fnames, kind = 1, Iendian=1, BIGLONG=FALSE ,
HEADONLY=FALSE,  PLOT = -1, RAW=FALSE)

JGET.seis(fnames, kind = 1, Iendian=1, BIGLONG=FALSE ,
HEADONLY=FALSE,  PLOT = -1, RAW=FALSE)

Arguments

Details

"kind" can be numeric or character: options are 'RDS', 'RDATA', 'SEGY', 'SAC', corresponding to (-1, 0, 1, 2).

Uses readBin to extract data in SAC/SEGY format. User must know what kind of machine the data was created on for I/O purposes.

If data was created on a little endian machine but is being read on big endian machine, need to call the endian "swap" for swapping.

Iendian can be a vector if input files have different endian-ness.

If data was created on a machine with LONG=4 bytes, be sure to call the program with BIGLONG=FALSE.

The data returned is a list of lists, each element is one trace not necessarily related to the other traces in the list.

Once the data is read in, use prepSEIS to reformat the data into a list more amenable to further analysis in RSEIS.

See examples below for different cases.

Value

List containing the seismic data and header information. Each trace consists of a list with:

Note

The easiext way to process datais to convert the data to an R-format type, using either save (kind=0) or saveRDS (kind=-1). If these are used then I/O is simple.

OLDER:

Information in the file names is ignored, so be sure to modify headers prior to using this method of extracting meta-data. (Or modify the meta data from the file names after reading in the data.)

For SEGY files, in LINUX-UNIX, use: rename, segymod (PASSCAL) to modify the headers

JGET.seis extracts digital seismic data from binary files stored in the file system. The program uses readBin for I/O and passes data back to R. Currently SAC, SEGY formats are installed but it is easy to extend. AH format is available for LINUX systems, but there were problems compiling in WINDOWS and MACOS so this feature was removed.

A filter for mseed format is currently being developed. Could use package 'IRISSeismic'

Author(s)

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

See Also

plotJGET, JSAC.seis , prepSEIS, Mine.seis

Examples

data(GH)

DD = data.frame(GH$info)

#### get only vertical traces
WV = which( GH$COMPS=='V' )
L1 = length(WV)

GIVE = vector(mode='list')

for(j in 1:L1 )
{
    i = WV[j]
    AA = DD[i,]
    GIVE[[j]] = list(fn = AA$fn, sta =GH$STNS[i] , comp = GH$COMP[i], 
                     dt = AA$dt, DATTIM = AA, N = AA$n1, units = NA, 
                     coords = NA, amp = GH$JSTR[[i]]  )
}
#### par(mfrow=c(length(GIVE) , 1) )
#   for(i in 1:length(GIVE) ) { plotGH(GIVE[[i]])   }
tdir = tempdir()
for(i in 1:length(GIVE) )
{
    sig = GIVE[[i]]
  d1 =  dateStamp(sig$DATTIM, sep='_')
   nam1 =  paste(d1,sig$sta, sig$comp, sep='_')
    nam2 = paste0(nam1, '.RDS')
    nam3 = paste(tdir, nam2, sep='/')
    saveRDS(file=nam3, sig) 
    }
########################  Now  read files and make the DataBase:
LF  = list.files(path=tdir,pattern='.RDS',  full.names=TRUE)


Gseis = GET.seis(LF, kind = -1, Iendian=1, BIGLONG=FALSE ,
HEADONLY=FALSE,  PLOT = -1, RAW=FALSE)

zed <- prepSEIS(Gseis)

####   plot the data, and interact with the data
swig(zed, sel=which(zed$COMPS=='V'),  SHOWONLY=0)

if(interactive()){ plotJGET(Gseis) }

###   for data created on UNIX (SUN) but read on linux:
###   S1   <-  GET.seis(Lname, kind = 1, Iendian="swap", BIGLONG=FALSE, PLOT = -1)

###   for data created on linux (32 bit)  but read on linux 64 bit:
###  S1   <-  GET.seis(Lname, kind = 1, Iendian="little", BIGLONG=FALSE, PLOT = -1)

###   for SEGY data created on linux (64 bit)  but read on linux 32 bit:
###  S1   <-  GET.seis(Lname, kind = 1, Iendian="little", BIGLONG=TRUE, PLOT = -1)


###   for SAC data created on MAC-OS (64 bit)  but read on linux 32 bit:
###  S1   <-  GET.seis(Lname, kind = 2, Iendian="swap", BIGLONG=TRUE, PLOT = -1)

Auto-Regressive AIC estimate of arrival time

Description

Auto-Regressive AIC for arrival estimate, signal detection

Usage

GETARAIC(z4, DT = 0.008, Mar = 8, O1 = 2, O2 = 0.2, WW = 2, T1 = 1, PLOT = FALSE)

Arguments

Details

Method of Sleeman for automatic phase determination.

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Sleeman

See Also

PSTLTcurve

Examples

data(CE1)
plot(CE1$x, CE1$y, type='l')

Xamp = CE1$y[CE1$x>4.443754 & CE1$x<6.615951]
Mar=8
z4 = Xamp
DT = CE1$dt
 T1 = 50

O1 = 10*DT
O2 =  10*DT
WW = 10*DT
Nz4 = length(z4)


 araict = GETARAIC(Xamp, DT=CE1$dt, Mar=8, T1=T1, O1=O1, O2=O2, WW=WW,   PLOT=TRUE)

Earthquake Seismic Data

Description

Example of seismic data structure. Geothermal Earthquake.

Usage

 data(GH) 

Format

List, consisting of:

JSTR

list of digital seismic data traces

STNS

vector of stations

dir

directory

ifile

original file names

COMPS

Component names, V N E, e.g.

OCOMPS

Old Component names

dt

vector of delta-t, sampling time intervals

KNOTES

Notes for plotting on panels

info

List, detailed information about traces, including

dat

not used

nn

Number of traces

ex

time axis for plotting

pcol

colors for plotting

ok

which traces are okay

wintim

window span time, seconds

ftime

alphanumeric time stamp

pickfile

pickfile, see below

velfile

velocity model list

stafile

station information list including lat, lon, z

aname

source name for loading

UWFILEID

event ID number

The info list consists of:

fn

file name

name

identification name

yr

start year

jd

start julianday

mo

month

dom

day of month

hr

hour

mi

minute

sec

second

msec

millisecond

dt

delta-t

t1

time 1

t2

time 2

off

offset

n1

number of samples

n2

not used

n3

not used

n

number of samples

The pickfile consists of:

LOC

list(yr, jd, mo, dom, hr, mi, sec, lat, lon, z, mag, gap, delta , rms, hozerr)

MC

list(az1, dip1, az2, dip2, dir, rake1, dipaz1, rake2, dipaz2, F=list(az, dip), G=list(az, dip), U=list(az, dip), V=list(az, dip), P=list(az, dip), T=list(az,dip),sense,M=list( az1, d1, az2, d2, uaz, ud, vaz, vd, paz, pd , taz, td), UP=TRUE, icol=1, ileg, fcol='red', CNVRG, LIM =c(0,0,0,0))

STAS

list(tag, name, comp, c3, phase, sec, err, pol, flg , res)

LIP

vector, length=6

H

list(yr,mo,dom,hr,mi,sec,lat,lon,z,mag)

N

name card

E

list(rms,meanres,sdres,sdmean, sswres,ndf,fixflgs,sterrx,sterry,sterrz,sterrt,mag,sterrmag)

filename

file name

PICKER

Name of Picker

UWFILEID

numeric ID

winID1

win format ID

comments

Vector of comments

OSTAS

Old station names

References

Lees, J.M., 2004. Scattering from a fault interface in the Coso geothermal field. Journal of Volcanology and Geothermal Research, 130(1-2): 61-75.

Examples

data(GH)

GLUE.GET.seis

Description

Once a database has been mined this program re-arranges the seismograms and creates a structure used in other programs.

Usage

GLUE.GET.seis(GG)

Arguments

Value

structure of seismograms glued together

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

Mine.seis

GLUEseisMAT

Description

Find duplicated stations in a matrix and fill in the traces that are continuations, return the new matrix and the vector duplicates

Usage

GLUEseisMAT(GFIL)

Arguments

Value

New List of data and headers with same sensors/components glued together

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

Mine.seis

Read 1D velocity model

Description

Read in a velocity model

Usage

Get1Dvel(infile, PLOT = TRUE)

Arguments

Details

Reads Velocity model from a text file

Value

LIST:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

Plot1Dvel, Comp1Dvel, Comp1Dvels, travel.time1D

Examples

data(VELMOD1D)

Get1Dvel(VELMOD1D, PLOT=TRUE)

 

Distance Along Great Circle Arc

Description

Distance Along Great Circle Arc in degrees, kilometers

Usage

GreatDist(LON1, LAT1, LON2, LAT2, EARTHRAD= 6371)

Arguments

Value

LIST:

Author(s)

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

Examples

###  get distance between London, England and Santiago, Chile
london  <-  c(51.53333, -0.08333333)
santiago  <-  c(-33.46667, -70.75)

GreatDist(london[2], london[1], santiago[2], santiago[1])

Instrument Frequencies

Description

Vector of frequencies

Usage

INSTFREQS(b, a, w)

Arguments

Details

h = jpolyval(b,s) / jpolyval(a,s)

Value

h = jpolyval(b,s) / jpolyval(a,s)

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

K  <-  PreSet.Instr()
b  <-  K[[1]]$zeros
a  <-  K[[1]]$poles
INSTFREQS(b, a, 2*pi*12)

Instrument Respnse Function

Description

Extract Instrument Response from Poles and Zeros

Usage

INSTresponse(Kal, key, ff, tt = tt, plotkey = NULL)

Arguments

Details

response is fourier transform of delta function run through the filter

Value

List:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Adapted from Ken Creager's Matseis

See Also

deconinst

Examples

########  set list of possible instruments:
Kal  <-  PreSet.Instr()
###  get instrument reponse for first in list:
resp1   <-  INSTresponse(Kal, 1, c(0,100)  , tt=c(1,0.008), plotkey=TRUE)
###  plots amplitude and phase

Inverse Fourier Transform

Description

Inverse Fourier Transform

Usage

INVRft(G, n, tstart, dt)

Arguments

Details

G is a vector spectrum evaluated at positive and negative frequencies as defined by makefreq. tstart, dt and n define the output time vector as described above.

g is the Inverse Fourier Transform of G scaled by dt. time shift theorem has been used to account for time not starting at t=0.

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

makefreq, FRWDft, INSTresponse

Examples

zil  <-  rnorm(300)
fss  <-  FRWDft( zil, length(zil), 0, 0.004)   
INVRft(fss$G, length(zil), 0, 0.004)

Gray scale Color Palette

Description

generate a gray scale color palette

Usage

JBLACK(n, acol=rgb(0,0,0))

Arguments

Details

Creates a black color palette suitable for replacing rainbow for B/W color plots. This is inserted in case user needs to completely elliminate color from a plot that uses color palettes for fixing colors.

Value

n characters used for color palette

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

shade.col, rainbow, colors

Examples

pal <- JBLACK(100)

Gray scale Color Palette

Description

generate a gray scale color palette

Usage

JGRAY(n)

Arguments

Details

Creates a grey scale color palette suitable for replacing rainbow for grey shade plots.

Value

n characters used for color palette

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

shade.col, rainbow, colors

Examples

pal <- JGRAY(100)
data(volcano)
image(volcano, col=pal)

JSAC.seis

Description

Read SEGY/SAC format binary data

Usage

JSAC.seis(fnames, Iendian = 1 , HEADONLY=FALSE,
BIGLONG=FALSE, PLOT = -1, RAW=FALSE)
JSEGY.seis(fnames, Iendian = 1 , HEADONLY=FALSE,
BIGLONG=FALSE, PLOT = -1, RAW=FALSE)

Arguments

Details

Uses readBin to extract data in SAC format. user must know what kind of machine the data was created on for I/O purposes.

For SEGY data the program is the same, although SEGY data does not have the problem of the BIGLONG so that is ignored.

For either code, a full header is returned, although the header for each format may be different.

Value

List containing the seismic data and header information. Each trace consists of a list with:

Note

SAC created on PC (windows) or LINUX machines typically will be in little endian format. SAC created on a SUN will be in big endian format. If you want to swap endian-ness , choose swap.

MAC uses different convention.

Iendian can be a vector if input files have different endian-ness.

SAC inserts -12345 for no data.

There are other issues regarding the size of long.

The units are often questionable and depend on the processing. The user should be careful and check to see that the proper conversions and multipliers have been applied.

Author(s)

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

See Also

Mine.seis, rseis2sac

Examples

#####  make some SAC files, then read them in
data(GH)
apath = tempdir()
##  setwd(apath)
## apath = 'TEMP'
J =  rseis2sac(GH, sel =1:5,  path = apath, BIGLONG =FALSE )

####### next read them in
Lname  <-  list.files(path=J , pattern='SAC', full.names=TRUE)

S1   <-  JSAC.seis(Lname, Iendian = .Platform$endian, BIGLONG =FALSE ,  PLOT = -1)

####  check just the first one
i = 1
  plotGH(S1[[i]])

Decimal Julian Day

Description

convert JD, HR, MIN SEC to Decimal Julian Day

Usage

Jtim(jj, hr = hr, mi = mi, sec = sec, yr=NULL, origyr=NULL)
JtimL(j)

Arguments

or

Details

Using a NULL value for yr gives the fractional julian day in a year. If yr is a legitimate year, and the origyr is provided, then the EPOCH number of days from origyr are added onto the fractional julian day. The default for origyr is 1972 for most of seismology.

If the dates span a new year, sometimes it is useful to use the earliest year as the origyr.

Value

Julian day

Author(s)

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

See Also

secdif

Examples

Jtim( 9 , hr= 14 , mi= 53 ,sec= 16.7807606880087 )

Jtim( 9 , hr= 14 , mi= 53 ,sec= 16.7807606880087, yr=2019, origyr=1972 )


########   or,
j = list(jd=9 , hr= 14 , mi= 53 ,sec= 16.7807606880087)

JtimL(j)

Volcano Seismic Data

Description

Seismic data from erupting Reventador Volcano. Vertical component only.

Usage

 data(KH) 

Format

KH = list( LOC=list(yr=0, jd=0, mo=0, dom=0, hr=0, mi=0, sec=0, lat=0, lon=0, z=0, mag=0, gap=0, delta=0 , rms=0, hozerr=0), MC=list(az1=0, dip1=0, az2=0, dip2=0, dir=0, rake1=0, dipaz1=0, rake2=0, dipaz2=0, F=list(az=0, dip=0), G=list(az=0, dip=0), U=list(az=0, dip=0), V=list(az=0, dip=0), P=list(az=0, dip=0), T=list(az=0, dip=0),sense=0,M=list( az1=0, d1=0, az2=0, d2=0, uaz=0, ud=0, vaz=0, vd=0, paz=0, pd =0, taz=0, td=0), UP=TRUE, icol=1, ileg="", fcol='red', CNVRG="", LIM =c(0,0,0,0) ),

STAS=list(tag="", name="", comp="", c3="", phase="", sec=0, err=0, pol="", flg=0 , res=0),

LIP=vector(length=6),

H=list(yr=0,mo=0,dom=0,hr=0,mi=0,sec=0,lat=0,lon=0,z=0,mag=0),

N=list(name=""),

E=list(rms=0,meanres=0,sdres=0,sdmean=0,sswres=0,ndf=0,fixflgs=0, sterrx=0,sterry=0,sterrz=0,sterrt=0,mag=0,sterrmag=0),

filename="",

PICKER="", UWFILEID="",winID1="",comments="", OSTAS="")

References

Lees, J. M., J. B. Johnson, M. Ruiz, L. Troncoso, M. Welsh, Reventador Volcano 2005: Eruptive Activity Inferred from Seismo-Acoustic Observation Journal of Volcanology and Geothermal Research in Press, 2007.

Examples

data(KH)
######  set SHOWONLY=FALSE for interactive
swig(KH, SHOWONLY=0)

Unwrap spectrum phase

Description

unwrap the phase spectrum so it does not wrap around

Usage

LocalUnwrap(p, cutoff = cutoff)

Arguments

Value

Unwrapped spectrum

Note

Algorithm minimizes the incremental phase variation by constraining it to the range [-pi,pi]

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

x  <-  1:512
amp <- sin(1*2*pi*x/16) + sin(2*2*pi*x/16)  + sin(3*2*pi*x/16)


spc <- fft(amp)

plot(Mod(spc), type='l')

angle  <-  Arg(spc)

plot(angle, type='l')

unang  <-  LocalUnwrap(angle, cutoff =pi )
plot(unang, type='l')

Interactive MTM driver

Description

MTM analysis of signals

Usage

MTM.drive(a, f1 = f1, f2 = f2, len2 = 1024, COL = 2, PLOT = FALSE,
PADDLAB = NULL, GUI = TRUE)

Arguments

Value

Graphical Side effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

See Also

plt.MTM0

Examples

data("GH")
sel <-  which(GH$COMPS=="V")

amp  <-  list()
dees  <-  list()
stamps  <-  list()

for( i in 1:3)
{
amp[[i]]  <-  GH$JSTR[[sel[i]]]
dees[i]  <-  GH$dt[sel[i]]
stamps[i]  <-  paste(GH$STNS[sel[i]], GH$COMPS[sel[i]])
}

 a  <-  list(y=amp, dt=dees, stamps=stamps)


 f1  <-  0.1

f2  <-  floor(0.33*(1/a$dt[[1]]))

speccol   <-  c('red', 'blue', 'purple')


MTM.drive(a, f1, f2, COL=speccol, PLOT=TRUE)

MTMdisp

Description

Display MTM displacement spectrum.

Usage

MTMdisp(a, f1 = f1, f2 = f2, len2 = 1024, PLOT = FALSE)

Arguments

Details

Uses Multi-taper estimate of spectrum and divides the spectrum by 1/(2*pi*f) to get integration of velocity seismogram.

Value

Returns displacement spectrum. Graphical Side effect.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

See Also

mtapspec

Examples

data(CE1)
xvel <- list(y=CE1$y[CE1$x>5.443754 & CE1$x<5.615951], dt=CE1$dt)

len2 <- next2(length(xvel$y))
Spec <- MTMdisp(xvel, f1=.01, f2=25, len2=len2, PLOT=FALSE )

Evolutive MTM Spectrum

Description

Time varying Auto-Regressive Spectrum (Gabor Transform) using MTM

Usage

MTMgabor(a, dt = 0, ppoint=95 , numf = 1024, Ns = 0, Nov = 0, fl = 0, fh = 10)

Arguments

Details

This is a spectrogram function similar to the Gabor Transform but uses the MTM (multi-taper method) for spectrum estimation. This is a non-interactive version of MTM.drive.

Value

List

Author(s)

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

References

Percival and Walden;

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

Percival, Donald B.,Walden, Andrew T. (1993):Spectral Analysis for Physical Applications,Cambridge University Press, Cambridge, 583p.

See Also

evolfft, evolMTM, MTM.drive, GETARAIC, doGABOR.AR, DOsgram, doGABOR.MTM

Examples

data(KH)
###   swig(KH)

Xamp <- KH$JSTR[[1]]
Nfft <- 1024   ###  fft length
Ns <- 512      ###  number of samples in a window
Nov <- 480    ###  number of samples of overlap per window
fl <- 0        ###  low frequency to return
fh <- 12     ###  high frequency to return
 dt  <-  KH$dt[1]

####  shorten the signal here, just for speed on the example:
sig = Xamp[37501:75001]

EV <-   MTMgabor(sig, dt = dt, numf =Nfft , Ns = Ns, Nov = Nov, fl = fl, fh= fh)

PE   <- plotevol(EV, log=1, fl=0.01, fh=fh, col=rainbow(100),
              ygrid=FALSE, STAMP="", STYLE="ar")

Plot Multi-taper Spectrum

Description

Plots output of MTM specturm

Usage

MTMplot(a, f1 = f1, f2 = f2, len2 = 1024, PLOT = FALSE, PADDLAB = NULL, GUI = TRUE)

Arguments

Details

Uses Lees' MTM code.

Value

list(len2=len2, f=f, f1=f1, f2=f2, displ=displ, ampsp=amp, flag=flag )

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

MTM.drive, MTMdisp, plt.MTM0

Mine a seismic data base to extract secions of time limited data

Description

Mine a seismic data base to extract secions of time limited data

Usage

Mine.seis(at1, at2, DB, grepsta, grepcomp, kind = 1, Iendian=1,
BIGLONG=FALSE, CHOP=TRUE, verbose=FALSE, chtoken=NULL, statoken=NULL, RAW=FALSE)

Arguments

Details

The data base is a list or dataframe containing the files names, the beginning time (t1) and ending time (t2) for each file in the data base. Mine.seis uses grep on the file names to extract specific files from the DB list.

Mine.seis needs to know what format the data was created in: little/big endian and the size of the LONG.

If data was created on a little endian machine but is being read on big endian machine, need to call the endian "swap" for swapping.

If data was created on a machine with LONG=4 bytes, be sure to call the program with BIGLONG=FALSE.

Use sysinfo to findout the system parameters for the local system. You need to know, however, what machine the binary files were created on.

In some situation the chanel name and the station name are not embedded in the file headers - in that case use the token from the file name.

Value

List of seismograms cut from the database

Note

The headers in the digital (segy or SAC) data files may not necessarily match the file names. Note that program JGET.seis extracts the station name and component name from the digital header and does not use the file name. It may be prudent to force the file names and header files to match prior to using Mine.seis. For SEGY files, in LINUX-UNIX, use: rename, segymod (PASSCAL) to modify the headers.

For SAC files, use sac software.

For R-based codes save the files in a format that has the relevant information (DAT format).

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

makeDB, GLUEseisMAT, JGET.seis, JSAC.seis, JSEGY.seis, sysinfo

Examples

data(GH)

DD = data.frame(GH$info)

#### get only vertical traces
WV = which( GH$COMPS=='V' )
L1 = length(WV)

GIVE = vector(mode='list')

for(j in 1:L1 )
{
    i = WV[j]
    AA = DD[i,]
    GIVE[[j]] = list(fn = AA$fn, sta =GH$STNS[i] , comp = GH$COMP[i], 
                     dt = AA$dt, DATTIM = AA, N = AA$n1, units = NA, 
                     coords = NA, amp = GH$JSTR[[i]]  )
}
#### par(mfrow=c(length(GIVE) , 1) )
#   for(i in 1:length(GIVE) ) { plotGH(GIVE[[i]])   }
tdir = tempdir()
for(i in 1:length(GIVE) )
{
    sig = GIVE[[i]]
  d1 =  dateStamp(sig$DATTIM, sep='_')
   nam1 =  paste(d1,sig$sta, sig$comp, sep='_')
    nam2 = paste0(nam1, '.RDS')
    nam3 = paste(tdir, nam2, sep='/')
    saveRDS(file=nam3, sig) 
    }
########################  Now  read files and make the DataBase:
LF  = list.files(path=tdir,pattern='.RDS',  full.names=TRUE)
DB = FmakeDB(LF, kind=-1)
IDB = infoDB(DB)


SAMPseis <- Mine.seis(IDB$at1, IDB$at2, DB, IDB$usta[1:3], IDB$ucomp[1], kind = -1   )

w <- swig(SAMPseis, SHOWONLY=0)

Morlet Wavelet

Description

Make Morlet Wavelet

Usage

Mmorlet(UB = -4, LB = 4, N = 256, plot = FALSE)

Arguments

Details

create a morlet function based on the matlab style routines

Value

time series list:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

scal2freqs, Rwave

Examples

mm  <- Mmorlet(-8, 8, 256)

get UW station file

Description

Match Picks with stations and return station structure

Usage

NEW.getUWSTAS(PICS)

Arguments

Details

matches Picks with stations

Value

STAS structure

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Plot Window Picks (WPX)

Description

adds picks to existing seismic section

Usage

NEWPLOT.WPX(t0, STNS, COMPS, YPX, FILL = FALSE, FORCE = TRUE, cex = cex, srt = srt)

Arguments

Details

Used in conjunction with swig program

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

########  no example available now

Delta-O18 isotpe record

Description

Data from Delta-O18 Isotope record of climate change. Periodicities of this data show the Milancovic cycles.

Usage

 data(OH) 

Format

OH = list( LOC=list(yr=0, jd=0, mo=0, dom=0, hr=0, mi=0, sec=0, lat=0, lon=0, z=0, mag=0, gap=0, delta=0 , rms=0, hozerr=0), MC=list(az1=0, dip1=0, az2=0, dip2=0, dir=0, rake1=0, dipaz1=0, rake2=0, dipaz2=0, F=list(az=0, dip=0), G=list(az=0, dip=0), U=list(az=0, dip=0), V=list(az=0, dip=0), P=list(az=0, dip=0), T=list(az=0, dip=0),sense=0,M=list( az1=0, d1=0, az2=0, d2=0, uaz=0, ud=0, vaz=0, vd=0, paz=0, pd =0, taz=0, td=0), UP=TRUE, icol=1, ileg="", fcol='red', CNVRG="", LIM =c(0,0,0,0) ),

STAS=list(tag="", name="", comp="", c3="", phase="", sec=0, err=0, pol="", flg=0 , res=0),

LIP=vector(length=6),

H=list(yr=0,mo=0,dom=0,hr=0,mi=0,sec=0,lat=0,lon=0,z=0,mag=0),

N=list(name=""),

E=list(rms=0,meanres=0,sdres=0,sdmean=0,sswres=0,ndf=0, fixflgs=0,sterrx=0,sterry=0,sterrz=0,sterrt=0,mag=0,sterrmag=0),

filename="",

PICKER="", UWFILEID="",winID1="",comments="", OSTAS="")

Note

The sample unit here is set to 0.3 which is 10000 times the correct sample rat.

References

Lees, J. M. and J. Park (1995): Multiple-taper spectral analysis: A stand-alone C-subroutine: Computers & Geology: 21, 199-236.

Examples

data(OH)
xx <-  swig( OH, sel=which(OH$COMPS == "V"), SHOWONLY=0)

XTR button to RSEIS

Description

Convert output of XTR button to RSEIS list.

Usage

P2GH(P1)

Arguments

Details

Running swig out after a selection of a window and the XTR button, one can create an RSEIS structure for further use in swig.

Value

RSEIS list

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig, prepSEIS

Examples

if(interactive()){
data(GH)

#### click twice and  select the XTR button
P1 <- swig(GH)

LH <- P2GH(P1)
L1 <- swig(LH)

}

Convert PDEs to List

Description

Converta list of individual PDE events to a list of lat, lon, z...etc

Usage

PDE2list(PDF)

Arguments

Details

uses getmem

Value

list

Author(s)

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

See Also

getmem, getPDEcsv, parse.pde,getPDEscreen

Documentation for swig

Description

Prints brief documentation for buttons in swig

Usage

PICK.DOC(w)

Arguments

Details

Buttons are defined in advance

Value

printed side effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

if(interactive() ) PICK.DOC(6:23)

plot all phase arrival picks

Description

plot all phase arrival picks

Usage

PLOT.ALLPX(t0, STNS, COMPS, YPX, PHASE = NULL, POLS = TRUE,
 FILL = FALSE, FORCE = TRUE, cex = cex, srt = srt)

Arguments

Details

for use in conjunction with PLOT.SEISN program

Value

Graphical Side Effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

PLOT.SEISN, swig

Examples

data(GH)

WPX = data.frame(GH$pickfile$STAS)
T0 = data.frame(GH$info)[1,]

sel = which(GH$COMPS=='V')
 PLOT.SEISN(GH, sel=sel)

PLOT.ALLPX(T0, GH$STNS, GH$COMPS, WPX, PHASE='P',FORCE=TRUE)

plot a matrix of several seismograms

Description

Matrix of several seismograms

Usage

PLOT.MATN(ascd, tim=1, dt=1, T1=0,  WIN=c(0,1), labs="",
notes=notes, sfact=1,ampboost=0,  shift=NULL, LOG="",
COL='red', add=1, AXES=1, units=NULL, VS=FALSE)

Arguments

Details

Plots a matrix of seismograms that each have the same starting time. For the AXES argument, 0 = no axes, AXES=1 plot scale for largest amplitude band and a multiplier for all others, AXES=2 left side, AXES=3 right side, AXES=4 alternate sides

Value

Graphical side effects and,

Author(s)

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

See Also

swig, matsquiggle, dowiggles, varsquiggle

Examples

dt  <-  0.001

 t  <-  seq(0, 6, by=0.001)

thefreqs  <-  seq(from=10, to=100, by=10)
theamps  <-   runif(length(thefreqs))

#   sample rate is 1000 Hz, 0.001 seconds 601 samples
 x  <-  NULL

for(i in 1:length(thefreqs))
{
x  <-  cbind(x, theamps[i]*sin(2*pi*thefreqs[i]*t))


}

PLOT.MATN(x,  dt = dt)

Plot Seismic Section

Description

Seismic traces are plotted on a panel horizontally.

Usage

PLOT.SEISN(GH, tim = 1, dt = 1, sel =c(1:4) , WIN =c(1,0) ,
labs=c("CE1") ,
notes = "CE1.V", subnotes=NA, tags ="CE1.V" ,
sfact = 1, LOG = "", COL = 'red', add = 1, pts = FALSE,
YAX = 1, TIT = NULL, SHIFT = NULL,COLLAPSE=FALSE,   rm.mean = TRUE, UNITS = "volts",
MARK = TRUE, xtickfactor = 1, vertline=NA )

Arguments

,

Details

panel of N traces are plotted. For YAX, default is YAX=1, plot an axis with no units label and scale all the traces to

Value

Graphical Side effect. list(n=nn, dy=dy, minS=minS, maxS=maxS, meanS=meanS, DX=range(tim[tflag]) )

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

data("GH")
m  <-  match( GH$STNS,    GH$stafile$name)
LATS  <-  GH$stafile$lat[m]
LONS  <-  GH$stafile$lon[m]
dees  <-  rdistaz( GH$pickfile$LOC$lat, GH$pickfile$LOC$lon, LATS, LONS)

sel  <-  which(GH$COMPS=="V")
sel  <-  sel[order(dees$dist[sel])]


###   set up good colors
pcols  <-  seiscols(GH)

### select only vertical components

PLOT.SEISN(GH, sel=sel)

GH$units  <-  rep("m/s", times=length(GH$KNOTES))
GH$pcols  <-  pcols


######  simple plot of GH structure
YN  <-  PLOT.SEISN(GH, WIN=c(5,12))

######  a color must be provided for all traces.
######  simple plot of GH structure, with selection and colors

YN  <-  PLOT.SEISN(GH, WIN=c(5,12), sel=sel, COL=rainbow(length(sel)) )


####   alternating Y axes
 YN  <-  PLOT.SEISN(GH, WIN=c(5,12) ,  dt=GH$dt[sel], sel=sel, sfact=1 ,
notes=GH$KNOTES[sel], YAX =3, UNITS = TRUE ,labs = GH$units[sel],
COL=pcols    , TIT="test")


#### Y  axes on same side
 YN  <-  PLOT.SEISN(GH, WIN=c(5,12) ,  dt=GH$dt[sel], sel=sel, sfact=1 ,
notes=GH$KNOTES[sel], YAX =2, UNITS = TRUE ,labs = GH$units[sel],
COL=pcols    , TIT="test")

Plot Seismic Section, travel time curve

Description

Seismic traces are plotted on a panel horizontally, with spacing according to distance from source.

Usage

PLOT.TTCURVE(GH, STAXY = NULL, DIST = c(0, 10), DY = 0.1,
tim = 1, dt = 1, sel = c(1:4), WIN = c(1, 0), labs = c("CE1"),
 notes = "CE1.V", tags = "CE1.V", sfact = 1, COL = "red",
add = 1, pts = FALSE, YAX = FALSE, TIT = NULL, SHIFT = NULL,
 rm.mean = TRUE, UNITS = "volts", MARK = TRUE)

Arguments

Value

Graphical Side effect.

list(n=nn, dy=dy, minS=minS, maxS=maxS, meanS=meanS, DX=range(tim[tflag]), DY=DY, DIST=DIST )

Note

This program is similar to PLOT.SEISN but traces are plotting with increasing distance from a set point. The distances are calculated prior to execution and passed as a vector or structure.

Author(s)

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

See Also

PLOT.SEISN

Plot picks on seismic record

Description

Add lines at phase arrival times

Usage

PLTpicks(picks, labs = NA, cols = NA)

Arguments

Details

picks = vector of times relative to the start of the plot (seismogram)

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

ex  <-  seq(from=0, to=4*pi, length = 200)

y  <-  sin(ex)
picks = c(0.5*pi, 2.3*pi)
plot(ex, y, type='l')


PLTpicks(picks, labs =c("P","P")  , cols =c('red','green')  )

PLTpicks(picks+2, labs =c("S","PKIKP") , cols ='blue'  )

Interactive Particle Motion Plot

Description

Plot Hodogram and show seismic particle motion

Usage

PMOT.drive(temp, dt, pmolabs = c("Vertical", "North", "East"), STAMP = "", baz = 0)

Arguments

Details

Input matrix should V, N, E. Baz is not implemented yet.

Value

Graphical Side Effect.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

data("GH")
sel <-  which(GH$STNS == "CE1")


YMAT  <-  cbind(GH$JSTR[[sel[1]]][1168:1500],
GH$JSTR[[sel[2]]][1168:1500],
GH$JSTR[[sel[3]]][1168:1500])

dt  <-   GH$dt[ sel[1] ]
ftime  <-  Zdate(GH$info, sel[1], 1)

if(interactive()){
PMOT.drive(YMAT, dt, pmolabs = c("Vertical", "North", "East"),
STAMP =ftime )

}

P-picking

Description

Add Pick Marks and Labels

Usage

PPIX(zloc, YN = NULL, col = 1, lab = "")

Arguments

Value

Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Short Term/Long Term curve

Description

ST/LT ratio curve for sutomated picking routines

Usage

PSTLTcurve(y, dt = 0.008, fwlen = 125, bwlen = 125,
 perc = 0.05, stretch = 1000, MED = 255, PLOT = FALSE)

Arguments

Value

list(flag=1, ind=ix, eye=eye, mix=mix, SNR=SNR, s2=s2, rat=therat)

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

require(stats)

data(CE1)
plot(CE1$x, CE1$y, type='l')

z <- CE1$y[ CE1$x>5.352622 & CE1$x<5.589836]
x <- CE1$x[ CE1$x>5.352622 & CE1$x<5.589836]

G <- PSTLTcurve(z, dt = CE1$dt, fwlen = 10,
bwlen = 10, perc = 0.05,
stretch = 10, MED = 11, PLOT = FALSE)

### get time from beginning of trace
tpick <- x[G$ind]
abline(v=x[G$ind], col='red', lty=2)

Plot 1D Velocity Model

Description

Plot 1D velocity model showing P-wave and S-wave layered models.

Usage

Plot1Dvel(v, tit = NULL, col=c('blue', 'brown'), ...)

Arguments

Details

Velocity model consists of a list of P and S depths and layer velocity values. See example below.

Value

Graphical Side effect

Note

Errors are not required, although future versions may include the plotting of error bars.

Author(s)

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

See Also

Get1Dvel, Comp1Dvel, Comp1Dvels, travel.time1D

Examples

VEL <- list()
VEL$'zp' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vp' <- c(1.1,2.15,3.2,4.25,5.3,6.25,6.7,6.9,7,7.2)
VEL$'ep' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'zs' <- c(0,0.25,0.5,0.75,1,2,4,5,10,12)
VEL$'vs' <- c(0.62,1.21,1.8,2.39,2.98,3.51,3.76,3.88,3.93,4.04)
VEL$'es' <- c(0,0,0,0,0,0,0,0,0,0)
VEL$'name' <- '/data/wadati/lees/Site/Hengil/krafla.vel'

Plot1Dvel(VEL, tit = 'This is an Example' )

Set up Standard Instrument Responses

Description

A set of standard known instrument responses.

Usage

PreSet.Instr()

Value

List of instrument responses. Each is a list:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

ReadSet.Instr

Examples

MYset <- PreSet.Instr()
MYset[[1]]

Dump a velocity model to an ascii file

Description

Dump a velocity model to an ascii file

Usage

Put1Dvel(vel, outfile)

Arguments

Value

Side Effects

Author(s)

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

See Also

Get1Dvel, travel.time1D

Seismic 1D Travel Time and raypath

Description

Travel time and raypath from source to reciever in 1D local model.

Usage

Ray.time1D(indelta, inhpz, instaz, inlay, ztop, vel)

Arguments

Details

Uses local 1D velocity model, not appropriate for spherical earth.

Value

list:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

travel.time1D, Get1Dvel

Examples

data(VELMOD1D)

v <-  VELMOD1D
indelta=23.;
 inhpz=7.;
 instaz=0.;
nz = length(v$zp)

tees  <-  travel.time1D(indelta, inhpz, instaz, nz , v$zp , v$vp)
rays  <-  Ray.time1D(indelta, inhpz, instaz, nz , v$zp , v$vp)



plot(rays$rnod[1:rays$nnod] , -rays$znod[1:rays$nnod],type="n",
 xlab="distance, km" , ylab="Depth, km")

abline(h=-v$zp, lty=2, col=grey(0.80) )
points(rays$rnod[1:rays$nnod] , -rays$znod[1:rays$nnod], pch=8, col='green')
lines(rays$rnod[1:rays$nnod] , -rays$znod[1:rays$nnod])
points(rays$rnod[rays$nnod] , -rays$znod[rays$nnod], pch=6, col='red', cex=2)
#####  to coordinate this in space, need to rotate about
#####    the line between source and receiver locations

Read Instrument Response in IRIS SEED format

Description

Read Instrument Response, poles and zeros, in IRIS SEED format.

Usage

ReadInstr(fn)

Arguments

Details

RSEIS currently has a function (ReadSet.Instr) to read pole/zero files, but it seems to expect a format different from what one gets from IRIS. This one is compatible with pole/zero files produced by rdseed when converting seed files from the DMC to SAC files.

Value

List of poles and zeros compatible for swig decon

Author(s)

Jake Anderson<ajakef@gmail.com>

See Also

ReadSet.Instr

Examples

######  create a SAC format response file:
temp.file= tempfile("PZ")
cat(file=temp.file, c(
"ZEROS 4",
"-999.0260  0.0000",
"POLES 6",
"-0.1480  0.1480",
"-0.1480  -0.1480",
"-314.1600  0.0000",
"-9904.8000  3786.0000",
"-9904.8000  -3786.0000",
"-12507.0000  0.0000",
"CONSTANT 4.540182e+20"),  sep='\n')


RESP <- ReadInstr(temp.file)

Read Instrument Response file

Description

Read in an instrument response file, or

Usage

ReadSet.Instr(file)

Arguments

Details

If file is a path to a file it is read in and processed. If file is a vector of character strings from a file that has already been read in, the file is processed directly. The tag names (ZEROS, POLES, SENSE, CONSTANT) can be upper,lower or mixed case. Alternative to SENSE = sensitivity, and CONSTANT=norm or knorm.

Value

list:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

###  in this case a file has already been read in:
CMG <- c(
      "ZEROS 2",
      "0.0000E+00 0.0000E+00",
      "0.0000E+00 0.0000E+00",
      "POLES 3",
      "-0.1480E+00 0.1480E+00",
      "-0.1480E+00 -0.1480E+00",
      "-50.0 0.0",
      "CONSTANT 1.0",
      "SENSE 800")

ReadSet.Instr(CMG)

Search Pix

Description

Search through pick strucutre to select phase arrivals

Usage

SEARCHPIX(KPX, IPX, tol = 0.5)

Arguments

Details

returns index vector of picks that satisfy: wn = which( abs(t2-t1) < tol)

Value

index vector

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

data(GH, package='RSEIS')
IPX = data.frame( uwpfile2ypx(GH$pickfile ) )
 
#######  take for example on pick
KPX = IPX[6, ]

SEARCHPIX(KPX, IPX, tol = 0.5)

Convert a SEIS list to a list of seismograms

Description

Convert a SEIS list to a list of seismograms each independent.

Usage

SEIS2list(GH)

Arguments

Details

The list returned is useful for editing or modifying the seismic data prior to swig.

Value

List of seismograms.

Author(s)

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

See Also

plotGH, swig

Examples

data(GH)
gg = SEIS2list(GH)
##   for(i in 1:length(gg) )
i = 1

{
plotGH(gg[[i]]); Sys.sleep(0.2) 
}

Minimum time in an RSEIS list

Description

Return date/time of trace with earliest date/time.

Usage

SEISNtime(GH)

Arguments

Value

Author(s)

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

Examples

data(GH)
SEISNtime(GH)

Select Buttons

Description

Select buttons interactively.

Usage

SELBUT(OPTS, onoff = 1, ocols = "white", default = "opt")

Arguments

Details

Used in swig. OPtions can be added, subtracted, deleted, or completely filled out based on interactive choice.

Value

character list of chosen options.

Author(s)

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

See Also

swig

Examples

if(interactive()){ 
STDLAB <- c("DONE", "QUIT", "zoom.out", "zoom.in", "SELBUT",
"FILT","UNFILT", "PSEL", "SGRAM", "WLET", "SPEC", "XTR" )
onoff = rep(0, length(STDLAB))
onoff[1:5] <- 1
SELBUT(STDLAB, onoff=onoff)

}

Pick stations and components interactively

Description

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

Usage

SELSTA(GH, 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

swig

Examples

data(GH)

SELSTA(GH, sel=1:7 , newdev = TRUE, STAY = FALSE)

Sensor Sensitivity from a known set of seismo/acoustic sensor

Description

From published sensitivities of seismic and acoustic sensors.

Usage

SENSORsensitivity(K = 1)

Arguments

Value

Sensitivity

Note

Current choices are: c("40T", "3T", "L28", "LD", "EL", "MC", "EL(SANGAY)")

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Johnson, J.B., R.C. Aster, M.C. Ruiz, S.D. Malone, P.J. McChesney,J.M. Lees, and P.R. Kyle, Interpretation and utility of infrasonic records from erupting volcanoes, J. Volc. Geoth. Res., 121 (1-2), 15-63, 2003.

Examples

SENSORsensitivity(3)
SENSORsensitivity(5)

stereonet representation of particle motion

Description

stereonet representation of particle motion

Usage

SNET.drive(intempmat, pmolabs = c("Vertical", "North", "East"), STAMP = "")

Arguments

Details

Interactive driver for partmotnet.

Value

Graphical Side effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

partmotnet

Examples

data("GH")

temp  <-  cbind(GH$JSTR[[1]], GH$JSTR[[2]], GH$JSTR[[3]])

atemp  <- temp[1168:1500, ]
SNET.drive(atemp, pmolabs = c("Vertical", "North", "East"), STAMP = "")

Interactive Spectrogram Driver

Description

Interactive Spectrogram Driver

Usage

SPECT.drive(Xamp, DT = 0.008, NEW = TRUE, STAMP = NULL ,
 freqlim=c(0, 20, 0, 20), winparams=c(4096,256, 204 ))

Arguments

Details

Interactive buttons are set internally. The parameters freqlim and winparams can be changed - these are simply the starting parameters for the initial display.

For winparams, the parameters are set to be appropriate for sample rates of typical seismic data, 100-125 samples per second. The number of points in the FFT are initially set to 4096 and the time window is set to 256. The overlap is calculated by subtracting 20 percent of the time window, so the overlap is 80 percent. Of course, the number of samples in a window must be less than the length of input time series.

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

plotevol, RPMG

Examples

data(CE1)

#########  Xamp = CE1$y[CE1$x>5.443754 & CE1$x<5.615951]

Xamp = CE1$y
plot(Xamp, type='l')


DT = CE1$dt

if(interactive() ) {
SPECT.drive(Xamp, DT = DT, NEW = TRUE, STAMP = NULL) }

Short term, long term average ratio

Description

Calculate the short term, long term average ratios of the squared amplitude in a time series.

Usage

STALTA(y, fwlen = 125, bwlen = 125)

Arguments

Details

Calculates the ratio of the forward/backard mean square sum.

Value

vector of ratios

Note

All filtering or pre and post analysis should be done outside of ratio curve estimate.

Author(s)

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

See Also

STLTcurve, PSTLTcurve

Examples

### easy example find P and S-wave arrivals, low noise
data(GH)
i = 6
z = GH$JSTR[[i]]

z.curve = STALTA(z, fwlen = 10, bwlen = 325)

ex = seq(from=0, length=length(z), by=GH$dt[i])
par(mfrow=c(2, 1) )
plot(ex, z, type='l')
plot(ex, z.curve, type = 'l' )

aa = peaks(z.curve, span = 11, do.pad = TRUE)
wa = which( aa & z.curve>50 )

abline(v=wa*GH$dt[i]  , col='red')
par(mfg=c(1,1) )
abline(v=wa*GH$dt[i]  , col='red')

Short-term/Long-term Average curve

Description

Get short-term average long-term verage ratio curve for picking

Usage

STLTcurve(y, dt = 0.008, fwlen = 125, bwlen = 125,
stretch = 1000, MED = 255, PLOT = FALSE)

Arguments

Details

Uses C-code and fast tanking algorithm written at UW

Value

sample to significant change in ratio curve

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

PSTLTcurve

Examples

data(CE1)

y  = CE1$y

DT = CE1$dt



sy = STLTcurve(y, dt=DT, fwlen =  25,  bwlen  = 25,
stretch=1000, MED=255, PLOT=FALSE)

par(mfrow=c(2,1))

plot(CE1$x, CE1$y, type='l')
plot(CE1$x,sy$rat, type='l')

Calculate Different Spectrum Types in Physical Units

Description

Spectrum is a wrapper function for stats::fft and RSEIS::mtapspec. For a given method (multi-taper spectrum or fft spectrum) and spectrum type (power, energy, amplitude, or phase), it returns the spectrum in physical units (obeying Parseval's theorem) and the corresponding frequency axis.

Usage

Spectrum(x, dt, one_sided = TRUE, type = 1, method = 1)

Arguments

Details

Phase spectrum is currently enabled only for method = 1 (fft). All possible energy and power spectra obey Parseval's relation (sum(s)*df ~= mean(x^2) for power; sum(s)*df ~= sum(x^2)*dt for energy). Parseval's relation may not be exact due to approximations used in making the spectrum one-sided or in the multi-taper method.

Input units are assumed to be volts and seconds; if other input units are used, adjust output units accordingly.

Value

List with following elements.

Author(s)

Jake Anderson

See Also

RSEIS::mtapspec stats::fft

Examples

## example time series
x = rnorm(1000)
dt = 0.01

## power spectrum, multi-taper method, one-sided
S = Spectrum(x, dt, type = 1, method = 2, one_sided = TRUE)
sum(S$spectrum) * S$df ## frequency-domain power
mean(x^2) ## time-domain power

## energy spectrum, fft method, two-sided
S = Spectrum(x, dt, type = 2, method = 1, one_sided = FALSE)
sum(S$spectrum) * S$df ## frequency-domain energy
sum(x^2) * dt ## time-domain energy

Get T1, T2

Description

Modify opick data frame and add T2=T1+dur

Usage

T12.pix(A)

Arguments

Details

Given t1 and duration, returns to structure, t2=t1+dur.

Value

pick data.frame with t2 as a member.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Taper Traces

Description

Taper traces in a seismic structure using a cosine function on the ends.

Usage

TAPER.SEISN(TH, sel = 1:length(TH$JSTR),  TAPER = 0.1 )

Arguments

Details

Seismic structure

Value

Seismic structure

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

butfilt

Examples

data("GH")
sel <-  which(GH$COMPS=="V")

sel  <-  1:3
 KF  <-  TAPER.SEISN(GH, sel = sel, TAPER=0.1)
swig(KF, sel=sel, SHOWONLY=0)

Convert to Cartesian coordinates

Description

Convert to cartesian coordinates

Usage

TOCART(az, nadir)

Arguments

Value

LIST

Author(s)

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

See Also

tocartL

Examples

TOCART(132, 69)

Threshold Adjuster

Description

determine cut off for ratio curve

Usage

Thresh.J(y, thresh)

Arguments

Details

Attempts to automatically optimize the threshold for automated picking. Used deep in picking algorithm.

Value

list(J=J, L=L)

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Sample Velocity Model

Description

Seismic Velocity Model for Coso California

Usage

data(VELMOD1D)

Format

LIST:

zp

vector of Tops of Layers, P-wave, (km)

vp

vector of velocities of Layers, P-wave,(km/s)

ep

errors for velocities, P-wave,(km/s)

zs

vector of Tops of Layers, S-wave, (km)

vs

vector of velocities of Layers, S-wave,(km/s)

es

errors for velocities, S-wave,(km/s)

name

character, name of model

descriptor

character vector description of model

Details

Velocity model from a text file

References

Wu, H., and J. M. Lees (1999), Three-dimensional P- and S-wave velocity structures of the Coso Geothermal Area, California, from microseismic traveltime data, J. Geophys. Res. 104, 13,217-13,233.

Examples

data(VELMOD1D)
Get1Dvel(VELMOD1D, PLOT=TRUE)

Velocity Seismogram

Description

Removes seismic instrument response and corrects for sensitivity of seismoc instrument, returning units of m/s rather than volts.

Usage

VELOCITY.SEISN(TH, sel = 1:length(TH$JSTR), inst = 1,
Kal = Kal,waterlevel = 1e-08,  FILT = list(ON = FALSE,
fl = 1/30, fh = 7, type = "HP", proto = "BU"))

Arguments

Details

Instrument responses are lists of poles and zeros for each instrument defined.

Value

Same as input list with new traces representing velocity versus volts

Author(s)

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

See Also

DISPLACE.SEISN, deconinst

Examples

Kal  <-  PreSet.Instr()
data(KH)

inst  <-  rep(0, length(KH$STNS))

VH  <-     VELOCITY.SEISN(KH, sel = 1, inst = 1,
Kal = Kal, FILT = list(ON = FALSE, fl = 1/30, fh = 7,
type = "HP", proto = "BU"))

Window a GH structure and extract data

Description

Window a time slice of seismic data and extract from a GH structure.

Usage

WINGH(GH, sel = 1, WIN = c(0,1) )

Arguments

Details

Preserves the data structure of the GH list. The purpose of this function is to extract a small subset of data from a larger data set (or longer time series) for subsequent processing.

Value

New GH structure.

Author(s)

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

See Also

swig

Examples

if(interactive()){
data(GH)

swig(GH, sel=which(GH$COMPS=="V" ))


jh = WINGH(GH, sel = which(GH$COMPS=="V" ), WIN = c(3.821281, 12.861820) )

swig(jh)
## compare with:
swig(GH, sel=which(GH$COMPS=="V" ), WIN = c(3.821281, 12.861820))



}

Extract data to RSEIS file

Description

swig Button Extract seismic data in RSEIS and save in GH format for exchange.

Usage

X2RSEIS(nh, g)

Arguments

Details

This function is used internally in RSEIS as a button in swig. The program should be run in a directory that has write permission.

The data is saved as a GH list.

Value

No value, writes to disk

Author(s)

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

See Also

XTR, X2SAC, swig

Examples

if(interactive()){
###  get data:
GH <-  Mine.seis(at1, at2, DB, NULL , NULL,
        kind = 1, Iendian=1)
w <- swig(GH, PADDLAB=c("X2SAC","X2RSEIS","YPIX" ) )

}

Extract Data to SAC format

Description

swig Button Extract seismic data in RSEIS and save in SAC format for exchange.

Usage

X2SAC(nh, g)

Arguments

Details

This function is used internally in RSEIS as a button in swig. The program should be run in a directory that has write permission.

Value

No value, writes to disk

Author(s)

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

See Also

XTR, X2RSEIS, swig

Examples

if(interactive()){
###  get data:
GH <-  Mine.seis(at1, at2, DB, NULL , NULL,
        kind = 1, Iendian=1)
w <- swig(GH, PADDLAB=c("X2SAC","X2RSEIS","YPIX" ) )

}

Buttons for swig

Description

defining functions for swig

Usage

XTR(nh, g)
NEXT(nh, g)
PREV(nh, g)
HALF(nh, g)
MARK(nh, g)
DOC(nh, g)
REFRESH(nh, g)
RESTORE(nh, g)
ZOOM.out(nh, g)
ZOOM.in(nh, g)
RIGHT(nh, g)
LEFT(nh, g)
SCALE(nh, g)
PSEL(nh, g)
FLIP(nh, g)
PTS(nh, g)
FILT(nh, g)
UNFILT(nh, g)
SPEC(nh, g)
WWIN(nh, g)
SGRAM(nh, g)
WLET(nh, g)
XTR(nh, g)
Pinfo(nh, g)
TSHIFT(nh, g)
RMS(nh, g)
LocStyle(nh, g)
CENTER(nh, g)
fspread(nh, g)
Xwin(nh, g)

Arguments

Details

Buttons can be defined on the fly.

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

Examples

if(interactive()){

MYFUNC<-function(nh, g)
  {
    print("pressed MYFUNC")
    g$sel
    d  <-  data.frame(list(stations=nh$STNS[g$sel],
                   components=nh$COMPS[g$sel]))
     print(d)        
    g$action  <-  "replot"
    invisible(list(global.vars=g))	
  }

STDLAB <- c("DONE", "QUIT", "SELBUT" , "PSEL", "MYFUNC" )
data(GH)
JJ  <-  swig(GH, sel=1:10, STDLAB=STDLAB)

}

PICK Buttons for swig

Description

defining functions for swig

Usage

YPIX(nh, g)
WPIX(nh, g)
NOPIX(nh, g)
REPIX(nh, g)
DELpix(nh, g)
PickWin(nh, g)
pADDPIX(nh, g, phase)
Ppic(nh, g)
Spic(nh, g)
Apic(nh, g)
POLSWITCH(nh, g, dir)
Pup(nh, g)
Pnil(nh, g)
Pdown(nh, g)
FILLPIX(nh, g)
RIDPIX(nh, g)
SEEPIX(nh, g)
ROT.RT(nh, g)
JustV(nh, g)
JustE(nh, g)
JustN(nh, g)
JustF(nh, g)
SHOW3(nh, g)

Arguments

Details

Buttons can be defined on the fly.

YPIX

Multiple picks on a panel

WPIX

window picks (start and end)

NOPIX

remove the picks

REPIX

un-remove the picks

DELpix

Delete pix near clicks

PickWin

Pick window for 3 component picking

pADDPIX

add picks

Ppic

P-wave arrival (only one per station)

Spic

S-wave arrival (only one per station)

Apic

acoustic-wave arrival (only one per station)

POLSWITCH

flip polarity

Pup

Polarity Up

Pnil

Polarity nil

Pdown

Polarity down

FILLPIX

Fill the pick from bottom to top of panel

RIDPIX

remove pick

SEEPIX

print current picks to screen

ROT.RT

Rotate to radial and transverse (need event and station locations

JustV

Display only vertical components

JustE

Display only east components

JustN

Display only north components

JustF

Display only infrasound (F) components

SHOW3

Display All 3 components

iNEXT

Used internally in PickWin to move to next station

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

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)
JJ  <-  swig(GH, sel=1:10, STDLAB=STDLAB)

}

Return difference in seconds

Description

Difference between two Date/Times (Julian Day)

Usage

YRsecdif(jd1, hr1, mi1, sec1, jd2, hr2, mi2, sec2, yr1 = 0, yr2 = 0)

YRsecdifL(T1, T2)

Arguments

Details

Returns T2-T1, year is used.

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

secdifL, secdif

Examples

T1  <-  list(jd=12, hr=13, mi=23, sec=21, yr=1964 )
T2  <-  list(jd=14, hr=23, mi=23, sec=2, yr=1976)

YRsecdif(T1$jd, T1$hr, T1$mi, T1$sec,  T2$jd, T2$hr, T2$mi, T2$sec,
1964, 1976)

####  or

YRsecdifL(T1, T2)

ZOOM SEISMIC Panel

Description

Zoom interactively on Seismic panel data.

Usage

ZOOM.SEISN(GH, sel = 1:length(GH$dt), WIN = NULL)

Arguments

Value

Seismic trace structure

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

data("GH")
sel <-  which(GH$COMPS=="V")

KF  <-  ZOOM.SEISN(GH, sel=sel, WIN = c(0 , 5)  )

if(interactive()){  swig(KF)
}

Date functions

Description

Make character vector from dates

Usage

Zdate(info, sel=1, t1=0, sep='_')
dateList(datevec)
dateStamp(datelist, sep='_')

Arguments

Details

Format date stamp for plotting and identification. Used for STAMP.

Value

character strings

Note

If using Zdate to create a file name, becareful about the separator. A colon in the file name on PC and MAC systems can be confusing for the OS.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig, dateStamp, ghstamp, filedatetime

Examples

data("GH")

sel <-  which(GH$COMPS == "V")

ftime  <-  Zdate(GH$info, sel[1:5], 1)

dvec  <-  c(2009, 134, 5, 14, 10, 32, 24.5, 0)
A  <-  dateList(dvec)
dateStamp(A, sep=".")

dateStamp(A, sep="_")

Add one pick to WPX file

Description

Add one pick to WPX file

Usage

addWPX(WPX, ppx)

Arguments

Details

Adds one pick to end of list.

Value

WPX list

Note

Uses, the last pick as a reference.

Author(s)

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

See Also

catWPX

Examples

s1 <- setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) 
s2 <- setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(1)) 

s3 <- addWPX(s1, s2)

Add points to a hodogram plot

Description

Add points to a hodogram plot

Usage

addpoints.hodo(nbaz, dt, sx, flag = 1:10, pch = 3, col = 1)

Arguments

Value

Graphical Side Effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

PMOT.drive, idpoints.hodo

add tix to plot

Description

Add tick marks to edge of plot

Usage

addtix(side = 3, pos = 0, tck = 0.005, at =c(0, 1), labels = FALSE, col = 2, ...)

Arguments

Value

Graphical Side Effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

par

Apply taper to seismic trace

Description

Apply taper to ends of a time series for spectrum analysis.

Usage

applytaper(f, p = 0.05)

Arguments

Details

10 percent taper is 5 percent on each end.

Value

Tapered time series.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

data(CE1)
Xamp  <-  CE1$y[CE1$x>5.443754 & CE1$x<5.615951]
Tamp  <-  applytaper(Xamp, p = 0.05)

Epoch Time Window

Description

Set a time window in Epoch days for extraction from a DB file

Usage

attime12(t1, t2 = t1, origyr = 1972, pre = 0, post = 0)

Arguments

Details

If t2 is missing, t1=t2.

Value

Author(s)

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

Examples

j1  <-  list(yr = 2005, jd= 214 , hr= 7 , mi= 1 ,sec= 0.5235)
j2  <-  list(yr=2005, jd= 214 , hr= 7 , mi= 1 ,sec= 0.5235+6)

at <-  attime12(j1, t2=j1, origyr=2005, pre=100, post=100)

### given an RSEIS format list
data(GH)
AT =  SEISNtime(GH)
ats = attime12(AT[[1]], t2 = AT[[2]],
   origyr =AT[[1]]$yr , pre = 0, post= 0)

Auto-Regressive Spectrum Estimate

Description

Auto-Regressive Spectrum Estimate

Usage

autoreg(a, numf = 1024, pord = 100, PLOT = FALSE, f1 = 0.01, f2 = 50)

Arguments

Value

LIST:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

fft, mtapspec, plt.MTM0

Examples

data(CE1)
Xamp <- CE1$y[CE1$x>5.443754 & CE1$x<5.615951]

ZIM <- autoreg(Xamp , numf=length(Xamp) , pord = 100, PLOT=FALSE,  f1=.01, f2=50)

Brune Modeling

Description

Modeling the Brune spectrum with Graphical Diagnostics

Usage

brune.doom(amp, dt = 1, f1 = 0.01, f2 = 15, PLOTB = FALSE, tit = "")

Arguments

Value

List:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees, J. M. and G. T. Lindley (1994): Three-dimensional Attenuation Tomography at Loma Prieta:Inverting t* for Q, J. Geophys. Res., 99(B4), 6843-6863.

Examples

data(CE1)
plot(CE1$x, CE1$y, type='l')
Xamp = CE1$y[CE1$x>5.443754 & CE1$x<5.615951]

BF = brune.doom( Xamp, CE1$dt ,f1=.5, f2=12 ,  PLOTB = TRUE)

Brune Earquake Model

Description

Calculate Forward Brune model

Usage

brune.func(freq, omega0, tstar0, fc, alpha, gamma)

Arguments

Details

Brune model.

Value

returns displacement spectrum from given parameters

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees, J. M. and G. T. Lindley (1994): Three-dimensional Attenuation Tomography at Loma Prieta:Inverting t* for Q, J. Geophys. Res., 99(B4), 6843-6863.

See Also

brune.doom

Search for Brune fit to displacement spectrum

Description

Model of the spectrum of a seismic arrival. Uses Brune's Model.

Usage

brune.search(infreq, inspec, f1, f2, omega0, fcorn, tstar0, gamma)

Arguments

Details

see paper by Lees and Lindley

Value

list(omega0=omega0,tstar0=tstar3[3] , fc=fcorn, alpha=0, gamma=gam3[3])

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees and Lindley

See Also

MTM

Examples

data(CE1)

####   set frequency range for modeling
f1  <-  0.01
f2  <-  14.0

##  set up data and parameters
amp  <-   CE1$y
len2  <-  2*next2(length(amp))
a  <-  list(y=amp, dt=CE1$dt)

Spec  <- MTMdisp(a, f1=f1, f2=f2, len2=len2, PLOT=FALSE )
 
lspec  <-  Spec$displ

###  get initial estimate of parameters
xc  <-  get.corner(  Spec$f , lspec, dt, f1, f2, PLOT=FALSE)

jmod  <-  brune.search(Spec$f, lspec, f1, f2,
                 xc$omega0, xc$corn, xc$tstar0, 2.0)

Butterworth filter

Description

Design and apply butterworth low/high/band pass filters.

Usage

butfilt(a, fl=0, fh=0.5, deltat=1, type="BP", proto="BU",
 npoles=5, chebstop=30.0, trbndw=0.3, RM=FALSE, zp=TRUE )

Arguments

Details

Creation of butfilt is a described by the following arguments:

LP

low pass

HP

high pass

BP

band pass

BR

band reject

BU

Butterworth

BE

Bessel

C1

Chebyshev type 1

C2

Chebyshev type 2

Arguments chebstop , trbndw are ignored for non-chebyshev filters. LP and HP filters are seet by specifying fl for HP filters and fh for LP filters, the other argumentin each case is ignored.

Mean values should be removed prior to calling this function, and then set RM=FALSE. This is true especially if tapering is applied prior to filtering.

Zero phase filter is achived by running filter back and forth. Otherwise a single pass is returned. This should be equivalent to package signal filtfilt (from MATLAB).

Value

Filtered time series.

Author(s)

originally written in FORTRAN by David Harris, converted to C and modified by Jonathan M. Lees<jonathan.lees@unc.edu>

References

Harris, D., 1990: XAPiir: A recursive digital filtering package. United States: N. p., Web. doi:10.2172/6416972.

See Also

AUGMENTbutfilt

Examples

data(CE1)

ts1  <-  CE1$y
zz  <-  butfilt(ts1, fl=1, fh=15,  deltat=CE1$dt, type="LP" ,  proto="BU",
npoles=5 )



###  try plotting:

### the above, by default, is zero phase.
#####  next filter with non-zero-phase
z2  <-  butfilt(ts1, fl=1, fh=15,  deltat=CE1$dt, type="LP" ,  proto="BU",
npoles=5, zp=FALSE )
ex = seq(from=0, by=CE1$dt, length=length(ts1))

plot(ex,  ts1, type='l')
lines(ex, zz, col='red')
lines(ex, z2, col='blue')


plot(ex[ex<0.5],  ts1[ex<0.5], type='l')
lines(ex[ex<0.5], zz[ex<0.5], col='red')
lines(ex[ex<0.5], z2[ex<0.5], col='blue')

Concatenate two WPX lists

Description

Concatenate (combine) two WPX lists.

Usage

catWPX(WPX, ppx)

Arguments

Details

Adds second list to the end of the first list.

Value

WPX list

Author(s)

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

See Also

addWPX, setWPX, checkWPX, cleanWPX, clusterWPX, repairWPX, saveWPX

Examples

s1  <-  setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) 

s2  <-  setWPX(name="BYE", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) 

s3  <-  catWPX(s1, s2)

Check WPX

Description

Check and verify WPX list for compliance.

Usage

checkWPX(wpx)

Arguments

Details

Perform several checks on completeness, length of components, station names, component names and date-times of the WPX lists.

Value

• 0 no problems

• 1 list incomplete

• 2 names incomplete

• 3 components incomplete

• 4 dates incomplete

Note

No action taken in the event an error occurs - see repairWPX to fix problems.

Author(s)

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

See Also

addWPX,catWPX, saveWPX,cleanWPX,clusterWPX,repairWPX,setWPX

Examples

s1 <- setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) 

s1$col=NULL

INteractive CHoice of Filter

Description

Choose Butterworth filter from a selection

Usage

choosfilt(thefilts = thefilts, ncol = 5)

Arguments

Details

Used for interactive choices in swig. See example below.

Value

filter parameter list:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

butfilt, RPMG

Examples

 thefilts <-
          list(flo=
               c(0.02, 0.02, 0.02, 0.02, 0.02,   0.02,
                 0.02, 0.02, 0.02,  0.02, 0.02,  0.02,
                 0.02,
                 1/2, 1/50,1/100, 1/100,1,1,
                 0.2, 15, 5, 2,1,
                 100),
               fhi=
               c(1/10, 1/6, 1/5, 1/4, 1/3, 1/2,
                 0.2,  0.5, 1.0,  2.0, 3.0,  4.0,
                 7.0,
                 8, 1/2.0,1/5.0,1/10.0,10,5,
                 7.0, 100, 100, 100,10,
                 100),
             type =
               c("LP","LP", "LP", "LP", "LP", "LP",
                 "LP","LP", "LP", "LP", "LP", "LP",
                 "LP",
                 "BP", "BP","BP","BP","BP","BP",
                 "HP", "HP","HP", "HP","HP",
                 "None"))

if(interactive() ) choosfilt(thefilts = thefilts, ncol = 5)

Draw a circle

Description

Draw a circle

Usage

circ()

Details

Draw a circle on new plot.

Value

Graphical Side Effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

net

Examples

circ()

Clean WPX

Description

Return an empty (clean) WPX.

Usage

cleanWPX()

Details

Returns an empty list with NA's and 0's

Value

WPX list

Note

Used internally.

Author(s)

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

See Also

addWPX, catWPX, checkWPX,repairWPX,clusterWPX,saveWPX, setWPX

Examples

s0 <- cleanWPX()

Clean up Pickfile structure

Description

Given a pickfile, clean out stations that do not ocnform

Usage

cleanpickfile(P)

Arguments

Details

stations with name="" are removed

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

EmptyPickfile

Examples

P <- EmptyPickfile()
cleanpickfile(P)

Plot a seimic trace colored in time

Description

Plot a seimic trace colored in time. useful for coordinating other plots to specific times along a seismic trace.

Usage

colorwig(x1, y1, COL = rainbow(100))

Arguments

Value

Graphical Side Effects

Author(s)

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

Examples

data(KH)

x <- KH$ex[KH$ex>95& KH$ex<125]
y <- KH$JSTR[[1]][KH$ex>95& KH$ex<125]


colorwig(x , y , rainbow(100))

Combine SEIS lists

Description

Combine 2 SEIS format lists into one list suitable for swig.

Usage

combineSEIS(IH, IV)

Arguments

Details

This will take two SEIS lists and merge them into one.

Value

SEIS list suitable for swig.

Author(s)

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

See Also

swig, Mine.seis, prepSEIS

Examples

#######  say you have 2 databases - extract from each:
####GH = Mine.seis(at1, at2, DB1, grepsta, grepcomp, kind = -1)
####JH = Mine.seis(at1, at2, DB2, grepsta, grepcomp, kind = -1)
#### merge the 2 structures

data(KH)

MH = KH

BH = combineSEIS(KH, MH)
###### plot and interact
swig(BH, SHOWONLY=TRUE )

Compare Envelopes

Description

calculate and plot signal envelopes.

Usage

comp.env(ex, Y, PLOT = TRUE, stamps = stamps)

Arguments

Details

Takes in an common x predictor and compares the envelopes of each column in the Y matrix. All the Y's must have the same length as ex.

Value

Graphical Side effects. returns envelope series.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

data("GH")

temp <- cbind(GH$JSTR[[1]], GH$JSTR[[2]], GH$JSTR[[3]])

atemp <- temp[1168:1500, ]
ex  <- seq(from=0,length=length(temp[1168:1500, 1]), by=GH$dt[1])

comp.env(ex, atemp, PLOT = TRUE, stamps = c("1","2", "3") )

HodoGram Plot

Description

HodoGram Plot

Usage

complex.hodo(nbaz, dt = dt, labs = c("Vertical", "North", "East"),
 COL = rainbow(100), STAMP = "")

Arguments

Value

sx = list graphical side effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

data("GH")

temp <- cbind(GH$JSTR[[1]][1168:1500], GH$JSTR[[2]][1168:1500],
GH$JSTR[[3]][1168:1500])

pmolabs <- c("Vertical", "North", "East")


sx <- complex.hodo(temp, dt=GH$dt[1]  ,labs=pmolabs,
STAMP="Example",  COL=rainbow(100) )

Contour Wavelet Transform

Description

Contour Wavelet Transform

Usage

contwlet(baha, Ysig, dt, clev = 0.75, NLEV = 12,
 zscale = 1, zbound = NULL, col = col, ygrid = FALSE,
 WUNITS = "Volts", PEAX = NULL)

Arguments

Value

Graphical Side Effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

plotwlet, wlet.do, wlet.drive

Convert Seismic data to RSEIS

Description

Convert Seismic in SAC or SEGY format to RSEIS native format.

Usage

convert2Rseis(FLS, NEWDIR = ".", kind = 1, Iendian = "little", BIGLONG =
FALSE, NEWsta = "", NEWcomp = "")

Arguments

Details

Converts the data to R format so it can be loaded with the load command. After this conversion, files should be loaded in subsequent calls by using kind=0.

Value

Side effects - creates new files on local system

Note

JGET.seis extracts digital seismic data from binary files stored in the file system. The program uses readBin for I/O and passes data back to R. Currently SAC, SEGY formats are installed but it is easy to extend. AH format is available for LINUX systems, but there were problems compiling in WINDOWS and MACOS so this feature was removed. A filter for mseed format is currently being developed.

Author(s)

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

See Also

JGET.seis, JSAC.seis , Mine.seis

Examples

Iendian = .Platform$endian
data(GH)

##########  create some SAC files:
apath = tempdir()
J =  rseis2sac(GH, sel = 1:5, path = apath, BIGLONG =FALSE )
####  get SAC file file names:
Lname  <-  list.files(path=J , pattern='SAC', full.names=TRUE)

#####   convert each file to a saved RSEIS file, saved in apath
####  reading in SAC files, kind=2
convert2Rseis(Lname, NEWDIR = apath, kind = 2, Iendian = Iendian, BIGLONG =
FALSE )
####  check if files are there
list.files(path=apath)

DateHour to List

Description

Convert a julian day+time to an RSEIS date list.

Usage

convertATT(at1, yr)

Arguments

Details

Calculates the data-list that RSEIS uses in calculations. The Month and Day-of-month are also returned.

Value

List with date and time

Author(s)

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

See Also

recdate, recdatel, dateList, dateStamp, filedatetime, rangedatetime, yeardate, Zdate, as.POSIXct

Examples

yr = 2014
j = 233.1234
convertATT(j, yr)

Moveout Correction

Description

Shift traces accoring to given moveout times

Usage

correct.moveout(GH, sel = 1, tims = 0)

Arguments

Details

Each trace listed in sel gets shifted forward or backward according to time in tims. This is useful for shifting traces according to a given moveout curve.

Value

RSEIS list structure returned with adjusted traces

Author(s)

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

Deconvolve instrument response from seismic data

Description

Deconvolve instrument response from seismic data

Usage

deconinst(data, sintr, KAL, key, Calibnew, waterlevel = 1e-08)

Arguments

Details

To avoid problems with dividing by very small numbers, water level is set =1.e-8

Value

deconvolved signal

Note

Calibnew(1)==3 then use a cos (hanning) taper

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

PreSet.Instr, ReadSet.Instr, INSTresponse

Examples

Kal <- PreSet.Instr()
amp <- rnorm(1024)
Calibnew <- c(1,1.0, 0.0 )

dy  <- deconinst(amp, 0.008, Kal,1, Calibnew, waterlevel=1.e-8)

Delete picks to WPX file

Description

Delete pick to WPX file

Usage

deleteWPX(WPX, ind=1)

Arguments

Details

Deletes one pick to end of list.

Value

WPX list

Note

Uses, the last pick as a reference.

Author(s)

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

See Also

addWPX, catWPX

Examples

s1 <- setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(5)) 
s2 <- setWPX(name="HI", yr=2011, jd=231, hr=4, mi=3, sec = runif(1)) 

s3 <- addWPX(s1, s2)

s4 <- deleteWPX(s3, ind=2:3)

Detail Pick on 3-component seismogram

Description

Pops up three components and prepares menu items for picking

Usage

detail.pick(y, ex, dt, TIT = "")

Arguments

Details

Creates interactive session for picking seismograms. Is called from swig.

Value

KSAVE = list(x=xsave, y=ysave)

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

swig

Examples

data(CE1)

  detail.pick(CE1$y, CE1$x, CE1$dt, TIT = "")
   

Remove trend from time series signal

Description

Remove trend from time series signal

Usage

detrend(x)

Arguments

Details

Removes the trend from a signal.

Value

vector with linear trend removed.

Author(s)

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

See Also

mean

Examples

dt <- 0.001

 t <- seq(0, 6, by=0.001)


y <- 5*sin(2*pi*10*t)

plot(t,y, type='l')

y <- y  + 3 * t
plot(t,y, type='l')


dy <- detrend(y)

plot(t,dy, type='l')

Distances from an RSEIS list

Description

Calculate euclidian distances from an RSEIS seismic data list, stations and event location.

Usage

distseisnXY(GH, sta=list(nam="", x=0 , y=0 , z=0) , LOC=list(x=0, y=0 , z=0))

Arguments

Value

Note

Locations of stations and source should be projected.

Author(s)

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

Examples

data(GH)
###   assume the lat lon in GH are x, y (projected)

staxy <-  list(nam=GH$stafile$name, x=GH$stafile$lon,
               y=GH$stafile$lat, z=GH$stafile$z)
 LOC <- list(x=GH$pickfile$LOC$lon, y=GH$pickfile$LOC$lat,
             z=GH$pickfile$LOC$z)
distseisnXY(GH, sta =staxy, LOC = LOC)

Gabor Transform with AR spectrum method

Description

Gabor Transform with AR spectrum method

Usage

doGABOR.AR(Xamp, DT = 0.008, multi = 1, scale.def = 0, TWIN = 2, TSKIP =
0.2, PCTTAP = 0.05, pord=100, PLOT=TRUE)

Arguments

Details

This is a spectrogram function similar to the Gabor Transform but uses the AR method for spectrum estimation.

Value

list

Note

The main difference between this and other similar calls is the way the windows are determined.

Author(s)

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

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

Percival, Donald B.,Walden, Andrew T. (1993):Spectral Analysis for Physical Applications,Cambridge University Press, Cambridge, 583p.

See Also

evolfft, evolMTM, MTM.drive, GETARAIC, doGABOR.AR, DOsgram, doGABOR.MTM

Examples

data(KH)
###   swig(KH)

Xamp <- KH$JSTR[[1]]
Xamp <- Xamp[57914:72989]


EV <- doGABOR.AR(Xamp, DT = KH$dt[1] , multi = 1, scale.def = 0,
TWIN = 2, TSKIP = 0.2, PCTTAP = 0.05)

Evolutive MTM Spectrum

Description

Time varying Auto-Regressive Spectrum (Gabor Transform) using MTM. This is a driver for MTMgabor.

Usage

doGABOR.MTM(Xamp, DT = 0.008, ppoint=95 , multi = 1,
 scale.def = 0, TWIN = 2, TSKIP = 0.2, PCTTAP = 0.05,  PLOT=TRUE)

Arguments

Details

This is a spectrogram function similar to the Gabor Transform but uses the MTM (multi-taper method) for spectrum estimation. This is a non-interactive version of MTM.drive.

Value

list output of MTMgabor:

Note

The main difference between this and other similar calls is the way the windows are determined.

Author(s)

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

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

Percival, Donald B.,Walden, Andrew T. (1993):Spectral Analysis for Physical Applications,Cambridge University Press, Cambridge, 583p.

See Also

MTMgabor, evolfft, evolMTM, MTM.drive, GETARAIC, doGABOR.AR, DOsgram

Examples

data(KH)
###   swig(KH)

Xamp = KH$JSTR[[1]]
Xamp = Xamp[57914:64914]

EV = doGABOR.MTM(Xamp, DT = KH$dt[1], multi = 1, scale.def = 0,
TWIN = 1, TSKIP = .1, PCTTAP = 0.05)

Dummy Button Function

Description

This is a dummy button function showing how buttons can be created on the fly

Usage

doMYBUTTS(butt = "", clicks = NULL, x = NULL)

Arguments

Plot wiggles

Description

Plot wiggles

Usage

dowiggles(AMAT, dt, dx)

Arguments

Value

graphical side effects

Author(s)

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

See Also

wiggleimage, matsquiggle

Examples

S1 = symshot1()

dowiggles(S1$smograms,S1$dt, S1$x)

Interpolate time series at higher sample rate.

Description

Interpolate a times series with a higher/lower sample rate for processes that are sensitive to low samples.

Usage

downsample(sig, dt=0.001, newdt=0.01, PLOT=FALSE )

Arguments

Details

Linear interpolation is performed between samples. If the newdt is an integer multiple of the old dt, The samples will not be modified.

Value

time series vector with new sample rate.

Author(s)

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

Examples

data(KH)
sig = KH$JSTR[[1]]
####  reduce samples from 125 (0.008)  to 25Hz (0.04) 
newdt = KH$dt[1]*5
sig2 = downsample(sig, dt = KH$dt[1], newdt = newdt )

L0 = length(sig)
L1 = length(sig2)

op <- par(no.readonly = TRUE)
par(mfrow=c(2,1) )
  plot.ts(ts(sig, deltat=KH$dt[1] ), xlab='s',
ylab='Amplitude', main=paste('Orignal', L0) )
grid()
     plot.ts(ts(sig2, deltat=newdt ), xlab='s',
ylab='Amplitude', main=paste('Downsample', L1) )
grid()
par(op)

Edit Data Base

Description

Edit, or remove items from an RSEIS data base after it has been read in.

Usage

editDB(DB, w)
pathDB(DB,  path1="", path2="")

Arguments

Details

The DB is a list. The program cycles through the elements of the list and removes all lnes that correspond to the idecies given in w.

Value

Returns a DB list

Note

A problem arises if the makeDB program reads in, or tries to read in files that have not data base header information. This program can eliminate these from the data base.

Author(s)

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

See Also

makeDB, infoDB

Examples

##########   create a data set and a DB
tdir = tempdir()
data(GH)

DD = data.frame(GH$info)
WV = which(GH$COMPS=='V')

L1 = length(WV)

######
GIVE = vector(mode='list')

for(j in 1:L1)
{
 i = WV[j]   
AA = DD[i,]

GIVE[[j]] = list(fn = AA$fn, sta =GH$STNS[i] , comp = GH$COMP[i], 
            dt = AA$dt, DATTIM = AA, N = AA$n1, units = NA, 
            coords = NA, amp = GH$JSTR[[i]]  )

}

#######  save files in the tempdir
for(i in 1:length(GIVE) )
{

    sig = GIVE[[i]]
  d1 =  dateStamp(sig$DATTIM)
   nam1 =  paste(d1,sig$sta, sig$comp, sep='_')
    nam2 = paste0(nam1, '.RDS')
    nam3 = paste(tdir, nam2, sep='/')
    saveRDS(file=nam3, sig) 
    }
LF  = list.files(path=tdir,pattern='.RDS', full.names=TRUE)

#######   make the database
cosoDB = FmakeDB(LF, kind=-1)

######  change the DB path:
path1<-tdir
path2<-"."

#######  change the path name of the trace files
newDB <- pathDB(cosoDB, path1, path2 )

Envelope Function with Hilbert Transform

Description

Envelope Function with Hilbert Transform

Usage

envelope(x)

Arguments

Details

Uses the hilbert transform to get the envelope function.

Value

vector of the absolute of the hilbert transform

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

data(CE1)
ev <- envelope(CE1$y)
plot(CE1$x, CE1$y, type='l')
lines(CE1$x,ev, col='red')

Evolutive Auto-Regressive Spectrum

Description

Time varying Auto-Regressive Spectrum (Gabor Transform)

Usage

evolAR(a, dt = 0, numf = 1024, pord = 100, Ns = 0, Nov = 0, fl = 0, fh = 10)

Arguments

Details

This is a spectrogram function similar to the Gabor Transform but uses the Auto-Regressive method for spectrum estimation.

Value

List

Author(s)

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

See Also

evolfft, evolMTM, MTM.drive, GETARAIC

Examples

data(KH)
###   swig(KH)

Xamp <- KH$JSTR[[1]]

dt <- KH$dt[1]
plot(seq(from=0, length=length(Xamp), by=dt), Xamp, type='l')
##  limit the trace, somewhat
Xamp <- Xamp[12670:22669]
plot(seq(from=0, length=length(Xamp), by=dt), Xamp, type='l')


Nfft<-1024   ###  fft length
Ns<-512      ###  number of samples in a window
Nov<-480    ###  number of samples of overlap per window
fl<-0        ###  low frequency to return
fh<-12     ###  high frequency to return

EV <- evolAR(Xamp, dt = dt, numf =Nfft , pord = 100, Ns = Ns,
       Nov = Nov, fl = fl, fh = fh)


PE <- plotevol(EV, log=1, fl=0.01, fh=fh,
      col=rainbow(100), ygrid=FALSE,
      STAMP="", STYLE="ar")

Evolutive Multi-taper Spectrum

Description

Time varying Multi-taper Spectrum (Gabor Transform)

Usage

evolMTM(a, dt = 0, numf = 1024, Ns = 0, Nov = 0, fl = 0, fh = 10)

Arguments

Details

This is a spectrogram function similar to the Gabor Transform but uses the MTM method for spectrum estimation.

Value

List

Author(s)

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

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

See Also

evolfft, MTM.drive

Examples

data(KH)
###   swig(KH)

Xamp <- KH$JSTR[[1]]

dt <- KH$dt[1]
plot(seq(from=0, length=length(Xamp), by=dt), Xamp, type='l')
##  limit the trace, somewhat
Xamp <- Xamp[12670:22669]
plot(seq(from=0, length=length(Xamp), by=dt), Xamp, type='l')


Nfft<-4096   ###  fft length
Ns<-512      ###  number of samples in a window
Nov<-480    ###  number of samples of overlap per window
fl<-0        ###  low frequency to return
fh<-12     ###  high frequency to return

EV <- evolMTM(Xamp, dt = dt, numf = Nfft, Ns = Ns, Nov = Nov, fl = fl, fh
= fh)


PE <- plotevol(EV, log=1, fl=0.01, fh=fh, col=rainbow(100), ygrid=FALSE,
STAMP="", STYLE="ar")


##   compare with:
## EVf <- evolfft(Xamp, dt = dt, Nfft =Nfft , Ns =Ns , Nov =Nov , fl =fl, fh = fh)

##  PE <- plotevol(EVf, log=1, fl=fl, fh=fh, col=rainbow(100), ygrid=FALSE,STAMP="", STYLE="fft")

Spectrogram fft

Description

Spectrogram using simple fft (Gabor Transform)

Usage

evolfft(a, dt = 0, Nfft = 0, Ns = 0, Nov = 0, fl = 0, fh = 10, pcttap =
0.05, adjust=TRUE )

Arguments

Details

This is a duplication of the spectrogram function in matlab which applies Welsh's Method. Each mini-window is tapered with a cosine window.

Value

List

Note

Parameter adjust is by default TRUE so that the choice of Ns, Nov, and kcol will be optimized, more or less. Set this logical to FALSE to force the function to use user input parameters.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

evolMTM, evolAR, MTM.drive

Examples

data(CE1)

####  plot signals
plot(CE1$x, CE1$y, type='l')

###  set parameters
Nfft<-1024   ###  fft length
Ns<-250      ###  number of samples in a window
Nov<-240     ###  number of samples of overlap per window
fl<-0        ###  low frequency to return
fh<-1/(2*CE1$dt)     ###  high frequency to return

########  calculate the evolutive fft (Gabor Transform)
EV <- evolfft(CE1$y, dt =CE1$dt , Nfft = Nfft, Ns =Ns , Nov =Nov , fl =fl
, fh = 25)

### plot image, but it does not look too interesting
image(EV$DSPEC)


### plot Gabor transform with special function
PE <- plotevol(EV, log=0, fl=0.01, fh=100, col=rainbow(100), ygrid=FALSE,
STAMP="", STYLE="fft")

Create a character string from a date

Description

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

Usage

filedatetime(orgtim, tims=0, datesep="-", timesep="_", secsep="_")

Arguments

Value

Author(s)

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

Examples

data(GH)

g1 <- getGHtime(GH)
g2 <- unlist(g1)

filedatetime(g2, 1)

Make Filter Stamp

Description

Create an text stamp describing a filter

Usage

filterstamp(fl=1/2, fh=10, type="BP")

Arguments

Details

If the frequency is less than 1, the period is displayed. For now only 3 digits are displayed. If the first argument, fl, is a list the parameters are extracted from the list and the other arguments are ignored.

Value

Author(s)

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

See Also

FILT.spread

Examples

fl <- c(0.01, 2)
fh <- c(10, 20)
type <- "BP"
filterstamp(fl, fh, type)

FILT<-list(ON=TRUE, fl=1/2, fh=12, type="HP", proto="BU")
filterstamp(FILT)

FILT<-list(ON=TRUE, fl=1/2, fh=12, type="BP", proto="BU")
filterstamp(FILT)

FILT<-list(ON=TRUE, fl=1/2, fh=12, type="LP", proto="BU")
filterstamp(FILT)

Integration in Frequency Domain

Description

Integration of seismic signal in Frequency Domain. Used for converting velocity seismogram to displacement.

Usage

finteg(data, dt)

Arguments

Value

Integrated time series signal

Note

To avoid problems with dividing by very small numbers, water level is set =1.e-8

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

   ##  waterlevel=1.e-8


    dfor5 <- rnorm(1000)

  idfor5 <- finteg(dfor5, 0.008)
 

Fix NA values.

Description

Replace NA values in a time series with mean values between end points of missing segments, or first and last real values in case the NA's are at the beginning or ends of traces.

Usage

fixNA(y)

Arguments

Details

fixNA searches for stretches of NA 's in a time series and replaces the NA values with numeric values based ont he two end points of each section.

Value

numeric vector with no NA values.

Note

function is used primarily in filter applications.

Author(s)

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

See Also

butfilt

Examples

## source("~/Site/TA_DATA/CODE/fixNA.R")

### last samples are NA
zig = rnorm(25)
zig[10:15] = NA

noNA = fixNA(zig)

### first samples are NA
zig = rnorm(25)
zig[1:5] = NA
noNA = fixNA(zig)


zig = rnorm(25)
zig[1:5] = NA
zig[21:25] = NA

noNA = fixNA(zig)



zig = rnorm(25)
zig[1] = NA
zig[21:25] = NA
zig[10:12] = NA

noNA = fixNA(zig)
cbind(zig, noNA)

fixUWstasLL

Description

Matches station locations to pickfile stations

Usage

fixUWstasLL(STAS, stafile)

Arguments

Details

Matches station locations to pickfile stations

Value

structure of station lat, lon, z

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Fix component names for uniformity

Description

Fix component names for uniformity

Usage

fixcompname(comp)

Arguments

Details

Translate the component names to something uniform that can be used for sorting and other functions.

Value

one of "V", "N", "E"

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

fixcompname("SHV")
fixcompname("SHE")

Fix Station Component Names

Description

Convert components to common names: V N E

Usage

fixcomps(oldcomps, SEGY = FALSE)

Arguments

Details

Attemps to convert irregular component names to common format for later processing.

Value

character vector

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

fixcompname

Examples

comp <- c("v", "e")
fixcomps(comp)

given julian day and year get month/day

Description

given julian day and year get month/day

Usage

fromjul(jul, yy)

Arguments

Value

list(mo=mm, dom=dd)

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

tojul

Examples

iyear <- 2001
 jul <- 233
inine <-   tojul(iyear,1,1);
ijul <-    inine + jul - 1;
fromjul( ijul, iyear);

add tic marks

Description

Add tic marks to plot

Usage

gaddtix(side = 3, pos = 0, tck = 0.005, at = c(0, 1),
labels = NULL, col = 2, addline = FALSE, ...)

Arguments

Value

Graphical side effect

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

par

Examples

plot(c(0,1), c(0,1), type='n', ann=FALSE, axes=FALSE)

gaddtix(side=1, pos=0,   tck=-0.01, at=seq(from=0, to=.5, by=.2) ,
labels=seq(from=0, to=.5, by=.2), col=1)

Ricker Wavelet

Description

Generate a ricker wavelet of a specfied frequency and length

Usage

genrick(freq, dt, nw)

Arguments

Value

ricker wavelet as a vector.

Note

Original code by Leonard Lisapaly (leonardl@fisika.ui.ac.id), converted to R by J.M. Lees.

Author(s)

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

Examples

dt <- 0.01
freq <- 16
nlen <- 35

G <- genrick(freq, dt, nlen)

tee <- seq(from=0, by=dt, length=length(G))

plot(tee, G, type='l')

Get Corner Frequency: Linear Model

Description

Search for low frequency asymptote, corner frequency, and fall off slope of seismic spectrum.

Usage

get.corner(INfreq, INspec, dt, f1, f2, PLOT = FALSE, VERBOSE = FALSE)

Arguments

Details

This routine does not assume any particular mathematical model. It searches for a three parameters that describe two lines that mimic the displacement spectrum. The search is done via least squares.

Value

Model of 3 parameters, best fit.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

brune.doom

Examples

data(CE1)

##   set frequency range for modeling for this high frequency data
##   we use f2 = 50, but for volcano data should be f2<15

f1 <- 0.01
f2 <- 50.0  

##  set up data and parameters
amp <-  CE1$y
len2 <- 2*next2(length(amp))
a <- list(y=amp, dt=CE1$dt)

Spec <- MTMdisp(a, f1=f1, f2=f2, len2=len2, PLOT=FALSE )

lspec <- Spec$displ

###  get initial estimate of parameters
xc <- get.corner(  Spec$f , lspec, CE1$dt, f1, f2, PLOT=FALSE)

Get Slepian Tapers

Description

Return a matrix of Slepian tapers

Usage

get.slepians(npoints = 900, nwin = 5, npi = 3)

Arguments

Details

This function onlyu returns the tapers for inspection. To apply the tapers use the function mtapspec.

Value

Matrix: nwin vectors of npoints Slepian tapers

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Lees, J. M. and Park, J., 1995: Multiple-taper spectral analysis: A stand-alone C-subroutine, Computers and Geology, 21(2), 199-236.

See Also

mtapspec

Examples

nwin <- 5
npi <- 3
npoints <- 900
sleps <- get.slepians(npoints, nwin, npi)

matplot(sleps, type='l', xlab="Index", ylab="Taper Amplitude")
legend('topleft', legend=1:nwin, lty=1:nwin, col=1:nwin)

Error Card

Description

Location Error Card

Usage

getEcard(ECARD)

Arguments

Value