Browse Source

array response plotting

pyrocko
braunfuss 2 years ago
parent
commit
07b0ee3283
11 changed files with 273 additions and 343 deletions
  1. +18
    -15
      src/cluster/ClusterProgs.py
  2. +1
    -2
      src/common/CommonProgs.py
  3. +25
    -90
      src/common/Program.py
  4. +1
    -1
      src/process/ProcessProgs.py
  5. +1
    -12
      src/process/array_crosscorrelation_v4.py
  6. +30
    -18
      src/process/main.py
  7. +20
    -18
      src/process/music.py
  8. +109
    -133
      src/process/sembCalc.py
  9. +16
    -52
      src/process/semp.py
  10. +2
    -2
      src/process/waveform.py
  11. +50
    -0
      src/tools/palantiri_init.py

+ 18
- 15
src/cluster/ClusterProgs.py View File

@ -1,7 +1,8 @@
import os
import sys
from palantiri.common import Basic
from palantiri.common import Globals
from palantiri.common import Globals
def Usage():
return 'arraytool.py cluster event_name [automatic array clustering and print arrayconfig]'
@ -14,21 +15,23 @@ class Intern(object):
def error(self, text):
print('\nError: ' + text + '\n')
print('Usage: ' + Usage())
sys.exit('\n*** Program aborted ***')
sys.exit('\n*** Palantiri clustering aborted ***')
def checkProgramParameter(self, nMinParams, nMaxParams):
eventName = sys.argv[2]
if len(sys.argv) < nMinParams: self.error('event name missing')
if len(sys.argv) > nMaxParams: self.error('Too many parameters')
if len(sys.argv) < nMinParams:
self.error('event name missing')
if len(sys.argv) > nMaxParams:
self.error('Too many parameters')
if not Globals.checkEventDirParameter(eventName): # Exists event directory ?
if not Globals.checkEventDirParameter(eventName):
s = 'Invalid parameter - <' + eventName +'>'
s += '\n '
s += 'Cannot find directory ' + Globals.EventDir()
self.error(s)
s = 'Invalid parameter - <' + eventName + '>'
s += '\n '
s += 'Cannot find directory ' + Globals.EventDir()
self.error(s)
def start(config):
@ -36,16 +39,16 @@ def start(config):
intern = Intern()
if sys.argv[1] == 'cluster':
intern.checkProgramParameter(3, 4)
intern.checkProgramParameter(3, 4)
workDir = [Globals.EventDir(), 'tmp2', 'cluster'] # ???
workDir = ['cluster']
cmd = "palantiri_cluster" + ' -f ' + Globals.EventDir()
workDir = [Globals.EventDir(), 'tmp2', 'cluster']
workDir = ['cluster']
cmd = "palantiri_cluster" + ' -f ' + Globals.EventDir()
else:
return False
return False
Basic.changeDirectory(workDir) # create working directory
Basic.changeDirectory(workDir)
os.system(cmd)
return True

+ 1
- 2
src/common/CommonProgs.py View File

@ -7,8 +7,7 @@ def start():
if sys.argv[1] != 'new_version':
return False
at = os.path.join(os.getcwd(),'Common', 'NewVersion.py')
at = os.path.join(os.getcwd(), 'Common', 'NewVersion.py')
os.system(sys.executable + ' ' + at)
return True

+ 25
- 90
src/common/Program.py View File

@ -1,102 +1,37 @@
import os, sys
import os
import sys
from palantiri.common import Basic, Globals, Logfile
class MainObj (object):
def __init__ (self, externClass, version, runTimeLog=None, errorLog=None) :
class MainObj(object):
self.version = version
self.runTimeLog = runTimeLog
self.errorLog = errorLog
self.extern = externClass
def __init__(self, externClass, version, runTimeLog=None, errorLog=None):
def run (self) :
self.version = version
self.runTimeLog = runTimeLog
self.errorLog = errorLog
self.extern = externClass
Logfile.init (self.runTimeLog, self.errorLog)
Logfile.setStartMsg (self.version)
def run(self):
if not self.extern.init() : Logfile.abort ()
Logfile.init(self.runTimeLog, self.errorLog)
Logfile.setStartMsg(self.version)
try :
ret = self.extern.process ()
if not self.extern.init():
Logfile.abort()
if ret : msg = 'Program finished'
else : msg = 'Program finished with error'
try:
ret = self.extern.process()
except KeyboardInterrupt :
msg = 'Program aborted by Control C'; ret = False
if ret:
msg = 'Palantiri finished'
else:
msg = 'Palantiri finished with error - maybe Sauron looked back?'
self.extern.finish ()
Logfile.showLabel (msg)
return ret
except KeyboardInterrupt:
msg = 'Gandalf made Pippin drop the Palantiri by Control C'
ret = False
class ExternMainObj (MainObj) :
def __init__ (self) :
MainObj.__init__ (self, self, version=xx, runTimeLog=None, errorLog=None)
# own init operations
def init (self) :
# own operations
return True
def process (self) :
#own operations
return True
def finish (self) :
#own operations
pass
def startTest (action, workingDir) :
events = 'NEAR-COAST-OF-NORTHERN-CHILE_2014-04-03T02-43-14'
dir = os.getcwd()
dir = Globals.setEventDir (os.path.join (dir,Globals.EVENTS, events))
args = [sys.argv[0], action, '-f', dir]
dir = Basic.changeDirectory (workingDir) # create and/or change working director
return args
from time import sleep
VERSION_STRING = '<Version String>'
class TestObj (MainObj) :
def __init__ (self, processFkt) :
MainObj.__init__ (self, self, VERSION_STRING, 'TestProgram_run.log', 'TestProgram.log')
self.processFkt = processFkt
def init (self) :
Logfile.add ('Init reached'); return True
def process (self) :
Logfile.add ('Process reached'); return self.processFkt()
def finish (self) :
Logfile.add ('Finish reached')
# -------------------------------------------------------------------------------------------------
def process1 () :
return True
def process2 () :
Logfile.add ('Press Ctrl C')
sleep (1000)
return True
def testAll () :
test1 = TestObj (process1)
test1.run ()
test2 = TestObj (process2)
test2.run ()
return True
self.extern.finish()
Logfile.showLabel(msg)
return ret

+ 1
- 1
src/process/ProcessProgs.py View File

@ -16,7 +16,7 @@ class Intern(object):
def error(self, text):
print('\nError: ' + text + '\n')
print('Usage: ' + Usage())
sys.exit('\n*** Program aborted ***')
sys.exit('\n*** Gandalf made Pippin drop the Palantiri by ***')
def checkProgramParameter(self, nMinParams, nMaxParams):


+ 1
- 12
src/process/array_crosscorrelation_v4.py View File

@ -203,7 +203,6 @@ class Xcorr(object):
t2 = UTCDateTime(self.Origin.time)
print('cant find')
found = False
for tr in traces:
@ -528,8 +527,6 @@ class Xcorr(object):
reftriggeronset = refp+onset-self.mintforerun
if int(self.Config['autoxcorrcorrectur']) == 1:
try:
refmarkername = os.path.join(self.EventPath,
('%s-marker') % (os.path.basename(
self.AF)))
@ -552,21 +549,13 @@ class Xcorr(object):
thrOff = float(self.Config['refsta'])
plotTrigger(trP, cft, thrOn, thrOff)
selection = float(raw_input('Enter self picked\
phase in seconds: '))
selection = float(input('Enter self picked phase in seconds: '))
tdiff = selection-self.mintforerun
refname = os.path.basename(self.AF)+'-shift.mseed'
trP.stats.starttime = trP.stats.starttime - selection
trP.write(os.path.join(self.EventPath, refname),
format='MSEED')
except Exception:
selection = 0.
refname = os.path.basename(self.AF)+'-shift.mseed'
trP.stats.starttime = trP.stats.starttime - selection - self.mintforerun
trP.write(os.path.join(self.EventPath, refname),
format='MSEED')
'''
tdiff = 0
trigger = trP.stats.starttime


+ 30
- 18
src/process/main.py View File

@ -16,7 +16,7 @@ from palantiri.tools.config import Event
from palantiri.process import ttt, sembCalc, waveform, times, deserializer
from palantiri.process.array_crosscorrelation_v4 import Xcorr, cmpFilterMetavsXCORR
from pyrocko import util, io
import subprocess
import numpy as num
if sys.version_info.major >= 3:
import _pickle as pickle
@ -35,7 +35,6 @@ ch.setFormatter(formatter)
logger.addHandler(ch)
evpath = None
def initModule():
global evpath
@ -595,9 +594,13 @@ def processLoop(traces=None, stations=None, cluster=None):
if switch == 0:
ff1 = filter.flo()
ff2 = filter.fhi()
if switch == 1:
ff1 = filter.flo2()
ff2 = filter.fhi2()
switchs = "l0"
else:
f1 = str('ff1 = filter.flo%s()'%filterindex+1)
ff1 = eval(f1)
f2 = str('ff2 = filter.fhi%s()'%filterindex+1)
ff2 = eval(f2)
switchs = "h1"
ps_wdf_emp = os.path.join(Folder['semb'],
"fobjpickle_process_emp\
_%s_%s%s"
@ -660,9 +663,11 @@ def processLoop(traces=None, stations=None, cluster=None):
ff1 = filter.flo()
ff2 = filter.fhi()
switchs = "l0"
if switch == 1:
ff1 = filter.flo2()
ff2 = filter.fhi2()
else:
f1 = str('filter.flo%s()'% str(filterindex+1))
ff1 = eval(f1)
f2 = str('filter.fhi%s()'% str(filterindex+1))
ff2 = eval(f2)
switchs = "h1"
ps_wdf = os.path.join(Folder['semb'], "fobjpickle_process_%s_%s%s" % (arrayname, ff1, ff2))
if cfg.Bool('load_wdf') is True:
@ -696,10 +701,13 @@ def processLoop(traces=None, stations=None, cluster=None):
workdepth), 'rb')
TTTGridMap, mint, maxt = pickle.load(f)
f.close()
if switch == 0:
step = cfg.step()
if switch == 1:
step = cfg.step_f2()
else:
s1 = str('cfg.step_f%s()')% str(filterindex+1)
step = eval(s1)
if cfg.UInt('forerun') > 0:
ntimes = int((cfg.UInt('forerun') +
cfg.UInt('duration')) / step)
@ -759,7 +767,7 @@ def processLoop(traces=None, stations=None, cluster=None):
nstats = stations_per_array
if switch == 0:
switchs = "l0"
if switch == 1:
else:
switchs = "h1"
arraySemb, weight, array_center = sembCalc.doCalc(
@ -792,7 +800,7 @@ def processLoop(traces=None, stations=None, cluster=None):
f.close()
if switch == 0:
switchs = "l0"
if switch == 1:
else:
switchs = "h1"
if cfg.pyrocko_download() is True:
@ -897,9 +905,11 @@ def processLoop(traces=None, stations=None, cluster=None):
ff1 = filter.flo()
ff2 = filter.fhi()
switchs = "l0"
if switch == 1:
ff1 = filter.flo2()
ff2 = filter.fhi2()
else:
f1 = str('filter.flo%s()'%filterindex+1)
ff1 = eval(f1)
f2 = str('filter.fhi%s()'%filterindex+1)
ff2 = eval(f2)
switchs = "h1"
ps_wdf = os.path.join(Folder['semb'],
"fobjpickle_process_%s_%s%s_\
@ -933,9 +943,11 @@ def processLoop(traces=None, stations=None, cluster=None):
ff1 = filter.flo()
ff2 = filter.fhi()
switchs = "l0"
if switch == 1:
ff1 = filter.flo2()
ff2 = filter.fhi2()
else:
f1 = str('filter.flo%s()'%filterindex+1)
ff1 = eval(f1)
f2 = str('ff2 = filter.fhi%s()'%filterindex+1)
ff2 = eval(f2)
switchs = "h1"
if cfg.Bool('bootstrap_array_weights') is False:
arraySemb, weight, array_center = sembCalc.doCalc(


+ 20
- 18
src/process/music.py View File

@ -10,6 +10,26 @@ from scipy import misc
from scipy import pi
def _spectrum(
metric,
antennas,
out,
thlo,thstep,thsz,
phlo,phstep,phsz
):
# Lower-level spectrum calculator with preprocessed arguments and
# pass-by-reference output array, for easier implementation with
# cython and being farmed out to multiple processes. (The problem is
# embarassingly parallel.
assert out.shape == (thsz,phsz)
for i in range(thsz):
th = thlo + i*thstep
for j in range(phsz):
ph = phlo + j*phstep
out[i,j] = _pmusic(metric,antennas,th,ph)
def sph2cart(sph):
"""
Convert one or more spherical coordinates to cartesian.
@ -275,21 +295,3 @@ def covar(samples):
def _pmusic(metric,antennas,theta,phi):
steer = sp.exp( 1j*antennas.dot(-aoa2prop_scalar(theta,phi)) )
return 1.0 / steer.conj().dot(metric).dot(steer).real
def _spectrum(
metric,
antennas,
out,
thlo,thstep,thsz,
phlo,phstep,phsz
):
# Lower-level spectrum calculator with preprocessed arguments and
# pass-by-reference output array, for easier implementation with
# cython and being farmed out to multiple processes. (The problem is
# embarassingly parallel.
assert out.shape == (thsz,phsz)
for i in range(thsz):
th = thlo + i*thstep
for j in range(phsz):
ph = phlo + j*phstep
out[i,j] = _pmusic(metric,antennas,th,ph)

+ 109
- 133
src/process/sembCalc.py View File

@ -92,7 +92,7 @@ def solve_timeshifts(maxp, nostat, nsamp, ntimes, nstep, dimX, dimY, Gmint,
low = -1.*cfg.Float('shift_max')
high = cfg.Float('shift_max')
if cfg.Bool('correct_shifts_empirical_manual_station_wise') is False:
if cfg.Bool('correct_shifts_empirical_station_wise') is False:
bounds = [(low, high)]
else:
for ref in refshifts:
@ -328,13 +328,15 @@ def writeSembMatricesSingleArray(SembList, Config, Origin, arrayfolder, ntimes,
dimX = cfg.dimX()
dimY = cfg.dimY()
if switch == 0:
winlen = cfg.winlen()
step = cfg.step()
if switch == 1:
winlen = cfg.winlen_f2()
step = cfg.step_f2()
else:
s1 = str('cfg.step_f%s()'% str(switch+1))
step = eval(s1)
w1 = str('cfg.winlen_f%s()'% str(switch+1))
winlen = eval(w1)
latv = []
lonv = []
@ -419,9 +421,11 @@ def collectSemb(SembList, Config, Origin, Folder, ntimes, arrays, switch,
if switch == 0:
winlen = cfg.winlen()
step = cfg.step()
if switch == 1:
winlen = cfg.winlen_f2()
step = cfg.step_f2()
else:
s1 = str('cfg.step_f%s()'% str(switch+1))
step = eval(s1)
w1 = str('cfg.winlen_f%s()'% str(switch+1))
winlen = eval(w1)
latv = []
lonv = []
@ -869,9 +873,11 @@ def collectSembweighted(SembList, Config, Origin, Folder, ntimes, arrays,
if switch == 0:
winlen = cfg.winlen()
step = cfg.step()
if switch == 1:
winlen = cfg.winlen_f2()
step = cfg.step_f2()
else:
s1 = str('cfg.step_f%s()'% str(switch+1))
step = eval(s1)
w1 = str('cfg.winlen_f%s()'% str(switch+1))
winlen = eval(w1)
latv = []
lonv = []
@ -1007,12 +1013,15 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
else:
dimX = cfg.dimX()
dimY = cfg.dimY()
if switch == 0:
winlen = cfg.winlen()
step = cfg.step()
if switch == 1:
winlen = cfg.winlen_f2()
step = cfg.step_f2()
else:
s1 = str('cfg.step_f%s()'% str(switch+1))
step = eval(s1)
w1 = str('cfg.winlen_f%s()'% str(switch+1))
winlen = eval(w1)
new_frequence = cfg.newFrequency()
forerun = cfg.Int('forerun')
@ -1280,9 +1289,15 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
for tracex, trl in zip(sorted(calcStreamMap.keys()), synthetic_traces):
if cfg.Bool('dynamic_filter') is False:
if switch == 0:
trl.bandpass(4,cfg_f.flo(), cfg_f.fhi())
elif switch == 1:
trl.bandpass(4,cfg_f.flo2(), cfg_f.fhi2())
ff1 = cfg_f.flo()
ff2 = cfg_f.fhi()
else:
f1 = str('cfg_f.flo%s()'% str(switch+1))
ff1 = eval(f1)
f2 = str('cfg_f.fhi%s()'% str(switch+1))
ff2 = eval(f2)
trl.bandpass(4,ff1, ff2)
calcStreamMap[tracex] = trl
@ -1441,6 +1456,7 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
################ Array Response calcualtion and visualization ########
if cfg.Bool('array_response') is True:
# Warning: for this functionality the obspy module is needed
from obspy.signal import array_analysis
from obspy.core import stream
from obspy.core.util import AttribDict
@ -1455,41 +1471,61 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
sl_s=0.03
# sliding window properties
# frequency properties
if switch == 0:
frqlow = cfg_f.flo()
frqhigh = cfg_f.fhi()
elif switch == 1:
frqlow = cfg_f.flo2()
frqhigh = cfg_f.fhi2()
prewhiten=0
# restrict output
semb_thres=-1e9
vel_thres=-1e9
forerun_r = util.time_to_str(util.str_to_time(Origin['time'])-forerun)
duration_r = util.time_to_str(util.str_to_time(Origin['time'])+duration)
stime = UTCDateTime(forerun_r)
etime = UTCDateTime(duration_r)
stream_arr = stream.Stream()
for trace in calcStreamMap.keys():
obspy_compat.plant()
tmins = []
coords = []
for trs in calcStreamMap.keys():
tr = calcStreamMap[trs]
tmins.append(tr.tmin)
if switch == 0:
frqlow = cfg_f.flo()
frqhigh = cfg_f.fhi()
else:
f1 = str('cfg_f.flo%s()'%switch+1)
frqlow = eval(f1)
f2 = str('cfg_f.fhi%s()'%switch+1)
frqhigh = eval(f2)
tr.bandpass(4, frqlow, frqhigh)
trace = obspy_compat.to_obspy_trace(tr)
center_lats = []
center_lons = []
for il in FilterMetaData:
if str(il) == str(trace):
calcStreamMap[trace].stats.coordinates = AttribDict({
'latitude': il.lat,
'elevation': il.ele,
'longitude': il.lon})
center_lats.append(float(il.lat))
center_lons.append(float(il.lon))
center_lat = num.mean(center_lats)
center_lon = num.mean(center_lons)
for il in FilterMetaData:
if str(il) == str(trs):
trace.stats.coordinates = AttribDict({
'latitude': float(il.lat),
'elevation': float(il.ele),
'longitude': float(il.lon)})
x, y = orthodrome.latlon_to_ne(center_lat, center_lon, float(il.lat), float(il.lon))
coords.append([float(il.lon), float(il.lat), float(il.ele)])
stream_arr.append(trace)
tmin = num.max(tmins)
forerun_r = util.time_to_str(tmin)
duration_r = util.time_to_str(tmin+150)
stime = UTCDateTime(forerun_r)
etime = UTCDateTime(duration_r)
stream_arr.append(calcStreamMap[trace])
results = array_analysis.array_processing(stream_arr, nsampr, nstepr,
sll_x, slm_x, sll_y, slm_y,
sl_s, semb_thres, vel_thres,
frqlow, frqhigh, stime,
etime, prewhiten)
timestemp = results[0]
relative_relpow = results[1]
absolute_relpow = results[2]
out = results
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colorbar import ColorbarBase
from matplotlib.colors import Normalize
@ -1508,15 +1544,15 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
# choose number of fractions in plot (desirably 360 degree/N is an integer!)
N = 36
N2 = 30
abins = np.arange(N + 1) * 360. / N
sbins = np.linspace(0, 3, N2 + 1)
abins = num.arange(N + 1) * 360. / N
sbins = num.linspace(0, 3, N2 + 1)
# sum rel power in bins given by abins and sbins
hist, baz_edges, sl_edges = \
np.histogram2d(baz, slow, bins=[abins, sbins], weights=rel_power)
num.histogram2d(baz, slow, bins=[abins, sbins], weights=rel_power)
# transform to radian
baz_edges = np.radians(baz_edges)
baz_edges = num.radians(baz_edges)
# add polar and colorbar axes
fig = plt.figure(figsize=(8, 8))
@ -1530,12 +1566,12 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
# circle through backazimuth
for i, row in enumerate(hist):
bars = ax.bar(left=(i * dw) * np.ones(N2),
height=dh * np.ones(N2),
width=dw, bottom=dh * np.arange(N2),
bars = ax.bar((i * dw) * num.ones(N2),
height=dh * num.ones(N2),
width=dw, bottom=dh * num.arange(N2),
color=cmap(row / hist.max()))
ax.set_xticks(np.linspace(0, 2 * np.pi, 4, endpoint=False))
ax.set_xticks(num.linspace(0, 2 * num.pi, 4, endpoint=False))
ax.set_xticklabels(['N', 'E', 'S', 'W'])
# set slowness limits
@ -1546,6 +1582,31 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
plt.show()
from obspy.imaging.cm import obspy_sequential
from obspy.signal.array_analysis import array_transff_wavenumber
# set limits for wavenumber differences to analyze
klim = 40.
kxmin = -klim
kxmax = klim
kymin = -klim
kymax = klim
kstep = klim / 100.
# compute transfer function as a function of wavenumber difference
transff = array_transff_wavenumber(stream_arr, klim, kstep, coordsys='lonlat')
# plot
plt.pcolor(num.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
num.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
transff.T, cmap=obspy_sequential)
plt.colorbar()
plt.clim(vmin=0., vmax=1.)
plt.xlim(kxmin, kxmax)
plt.ylim(kymin, kymax)
plt.show()
traces = calcStreamMap
if cfg.Int('dimz') != 0:
@ -1741,7 +1802,7 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
k = backSemb
TTTGrid = False
if cfg.Bool('correct_shifts_empirical_run') is True and cfg.Bool('correct_shifts_empirical_manual') is False and flag_rpe is True:
if cfg.Bool('correct_shifts_empirical_run') is True and flag_rpe is True:
trs_orgs = []
calcStreamMapshifted = calcStreamMap.copy()
@ -1792,91 +1853,6 @@ def doCalc(flag, Config, WaveformDict, FilterMetaData, Gmint, Gmaxt,
if TTTGrid:
if cfg.Bool('correct_shifts_empirical_run') is True and cfg.Bool('correct_shifts_empirical_manual') is True and flag_rpe is True:
start_time = time.time()
RefDict_stations = OrderedDict()
for s in range(0, nostat):
refshifts[s] = refshifts[s]*0.
RefDict_stations[s] = [stations[s].lat, stations[s].lon]
if sys.version_info.major >= 3:
fobjrefshift_stations = open(rp+'_stations', 'wb')
else:
fobjrefshift_stations = open(rp+'_stations', 'w')
pickle.dump(RefDict_stations, fobjrefshift_stations)
fobjrefshift_stations.close()
step_emp = cfg.Float('step_emp')
if cfg.UInt('forerun_emp') > 0:
ntimes = int((cfg.UInt('forerun_emp') + cfg.UInt('duration_emp'))/step_emp)
else:
ntimes = int((cfg.UInt('duration_emp')) / step_emp)
nsamp = int(winlen)
nstep = float(step)
Gmint = cfg.Int('forerun')
max_shifts = cfg.Float('shift_max')
semblance_max = 0. #TODO auf station basis
if cfg.Bool('correct_shifts_empirical_manual_station_wise') is True:
for niter in range(0,1000):
for s in range(0, nostat):
refshifts[s] = num.random.uniform(-max_shifts,max_shifts)
k = semblance(maxp, nostat, nsamp, ntimes, nstep, dimX, dimY, Gmint,
new_frequence, minSampleCount, latv, lonv, traveltimes,
traces, calcStreamMap, timeev, Config, Origin, refshifts,
nstats,
bs_weights=bs_weights, flag_rpe=True)
partSemb = k
partSemb = partSemb.reshape(ntimes, 1)
for a in range(0, ntimes):
semb_max = max(partSemb[a])
if semb_max > semblance_max:
semblance_max = semb_max
refshifts_best = refshifts
RefDict = OrderedDict()
for j in range(0, nostat):
RefDict[j] = refshifts_best[j]
if sys.version_info.major >= 3:
fobjrefshift = open(rp, 'wb')
else:
fobjrefshift = open(rp, 'w')
pickle.dump(RefDict, fobjrefshift)
fobjrefshift.close()
print("--- %s seconds ---" % (time.time() - start_time))
return partSemb, weight, array_center
else:
for shft in num.arange(-max_shifts, max_shifts, 0.01):
k = semblance(maxp, nostat, nsamp, ntimes, nstep, dimX, dimY, Gmint,
new_frequence, minSampleCount, latv, lonv, traveltimes,
traces, calcStreamMap, timeev+shft, Config, Origin, refshifts,
nstats,
bs_weights=bs_weights, flag_rpe=True)
partSemb = k
partSemb = partSemb.reshape(ntimes, 1)
for a in range(0, ntimes):
semb_max = max(partSemb[a])
if semb_max > semblance_max:
semblance_max = semb_max
array_shift = shft
RefDict = OrderedDict()
for j in range(0, nostat):
RefDict[j] = shft
if sys.version_info.major >= 3:
fobjrefshift = open(rp, 'wb')
else:
fobjrefshift = open(rp, 'w')
pickle.dump(RefDict, fobjrefshift)
fobjrefshift.close()
print("--- %s seconds ---" % (time.time() - start_time))
return partSemb, weight, array_center
else:
start_time = time.time()
if cfg.UInt('forerun') > 0:
ntimes = int((cfg.UInt('forerun') + cfg.UInt('duration'))/step)


+ 16
- 52
src/process/semp.py View File

@ -15,6 +15,7 @@ from pyrocko import obspy_compat
import math
from scipy.signal import coherence
from collections import OrderedDict
from palantiri.process import music
trace_txt = 'trace.txt'
travel_txt = 'travel.txt'
@ -133,16 +134,20 @@ def semblance(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
time, cfg, refshifts, nstats, bs_weights=bs_weights)
if cfg.Int('dimz') != 0:
return semblance_py_cube(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY,
mint, new_frequence, minSampleCount, latv_1,
lonv_1, traveltime_1, trace_1, calcStreamMap,
time, cfg, refshifts, bs_weights=bs_weights)
return semblance_py_cube(ncpus, nostat, nsamp, ntimes, nstep,
dimX, dimY, mint, new_frequence,
minSampleCount, latv_1, lonv_1,
traveltime_1, trace_1, calcStreamMap,
time, cfg, refshifts,
bs_weights=bs_weights)
if cfg.Bool('bp_music') is True:
return music(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY,
mint, new_frequence, minSampleCount, latv_1,
lonv_1, traveltime_1, trace_1, calcStreamMap,
time, cfg, refshifts, nstats, bs_weights=bs_weights)
return music_wrapper(ncpus, nostat, nsamp, ntimes, nstep, dimX,
dimY, mint, new_frequence, minSampleCount,
latv_1, lonv_1, traveltime_1, trace_1,
calcStreamMap, time, cfg, refshifts,
nstats,
bs_weights=bs_weights)
if flag_rpe is True:
return semblance_py_fixed(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY,
@ -251,47 +256,6 @@ def semblance_py_dynamic_cf(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY,
return abs(backSemb)
def music_doa(Y,n,d):
# Y <- the ULA data
# n <- the number of sources
# d <- sensor spacing in wavelengths
# doa -> the vector of DOA estimates
[m,N]= size(data)
# compute the sample covariance matrix
R=Y*Y.T/N
# do the eigendecomposition; use svd because it sorts eigenvalues
[U,D,V]=num.linalg.svd(R)
G=U[:,n+1:m]
C=G*G.T
a =[]
# find the coefficients of the polynomial in (4.5.16)
for kk in range(0,2*m-1):
a[kk]=sum(num.diag(C,kk-m))
ra=num.roots(a)
# find the n roots of the a polynomial that are nearest and inside the unit circle,
[dum,ind]=sort(abs(ra))
rb=ra(ind[1:m-1])
# pick the n roots that are closest to the unit circle
[dumm,I]=sort(abs(abs(rb)-1))
w=num.angle(rb(I[1:n]))
# compute the doas
doa= math.asin(w/d/math.pi/2.)*180/math.pi
return doa
def semblance_py(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
new_frequence, minSampleCount, latv_1, lonv_1, traveltime_1,
trace_1, calcStreamMap, time, cfg, refshifts, nstats,
@ -364,7 +328,7 @@ def semblance_py(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
# markers.append(m)
ibeg = max(0, t2ind_fast(tmin-tr.tmin, tr.deltat, snap[0]))
index_begins[str(j)+str(k)]= [ibeg, tmin]
index_begins[str(j)+str(k)] = [ibeg, tmin]
iend = min(
tr.data_len(),
@ -376,9 +340,9 @@ def semblance_py(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
index_steps.append(iend_step-iend)
index_window.append(iend-ibeg)
# for debug:
# trld.snuffle(trs_orgs, markers=markers)
'''
Basic.writeMatrix(trace_txt, trace, nostat, minSampleCount, '%e')
Basic.writeMatrix(travel_txt, traveltime, nostat, dimX * dimY, '%e')
@ -502,7 +466,7 @@ def semblance_py(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
def music(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
def music_wrapper(ncpus, nostat, nsamp, ntimes, nstep, dimX, dimY, mint,
new_frequence, minSampleCount, latv_1, lonv_1, traveltime_1,
trace_1, calcStreamMap, time, cfg, refshifts, nstats,
bs_weights=None):


+ 2
- 2
src/process/waveform.py View File

@ -167,8 +167,8 @@ def resampledummy(Waveform, end_frequence):
return Waveform
def processWaveforms(WaveformDict, Config, Folder, network, MetaDict, Event,
switch, Xcorr):
def processWaveforms_obspy(WaveformDict, Config, Folder, network, MetaDict,
Event, switch, Xcorr):
Logfile.red('Start Processing')


+ 50
- 0
src/tools/palantiri_init.py View File

@ -0,0 +1,50 @@
import os
import fnmatch
import sys
from optparse import OptionParser
import logging
import imp
import obspy.core
from palantiri.process import ProcessProgs
from palantiri.cluster import ClusterProgs
from palantiri.common import CommonProgs
from palantiri.tools import config as config
import palantiri
import shutil
logger = logging.getLogger(sys.argv[0])
logger.setLevel(logging.DEBUG)
def main():
if len(sys.argv) < 2:
print("workdir path/name missing")
quit()
foldername = sys.argv[1]
workfolder = os.path.join(os.getcwd(), './', foldername)
eventsfolder = os.path.join(os.getcwd(), './', foldername, 'events')
tttgridsfolder = os.path.join(os.getcwd(), './', foldername, 'tttgrid')
tmpfolder = os.path.join(os.getcwd(), './', foldername, 'tmpProcess')
if os.access(workfolder, os.F_OK) is False:
os.makedirs(workfolder)
os.makedirs(tmpfolder)
os.makedirs(tttgridsfolder)
os.makedirs(eventsfolder)
logger.info('\033[31m Working Super-FOLDER CREATED \033[0m \n')
else:
print("workdir already exists!")
quit()
dstfolder = foldername
dstfile = ('global.conf')
path = palantiri.__path__
src = os.path.join(path[0], 'skeleton', 'global.conf')
dst = os.path.join(dstfolder, dstfile)
logger.info('\033[31m Created work directory \
%s \033[0m \n' % (dstfolder.split('/')[-1]))
shutil.copy(src, dst)

Loading…
Cancel
Save