[task] implementing new methods for pdf comparison.

This commit is contained in:
sebastianp 2016-06-29 15:29:58 +02:00
parent 8714616d1b
commit a8b7eff561
4 changed files with 103 additions and 4 deletions

View File

@ -1 +1 @@
7c5a-dirty
8714-dirty

View File

@ -180,7 +180,7 @@ def picksdict_from_picks(evt):
try:
onsets = picks[station]
except KeyError as e:
print(e)
#print(e)
onsets = {}
mpp = pick.time
spe = pick.time_errors.uncertainty

View File

@ -2,9 +2,9 @@
# -*- coding: utf-8 -*-
import copy
import os
import numpy as np
import glob
import matplotlib.pyplot as plt
from obspy import read_events
@ -336,3 +336,79 @@ class PDFDictionary(object):
plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False)
return fig
class PDFstatistics(object):
def __init__(self, directory):
self.directory = directory
self.stations = {}
self.p_std = {}
self.s_std = {}
self.makeDirlist()
self.getData()
self.getStatistics()
#self.showData()
def makeDirlist(self, fn_pattern='*'):
self.dirlist = glob.glob1(self.directory, fn_pattern)
#print self.dirlist
self.evtdict = {}
for rd in self.dirlist:
#print rd
self.evtdict[rd] = glob.glob1((os.path.join(self.directory, rd)), '*.xml')
#print rd, self.evtdict[rd]
def getData(self):
arraylen = 0
for dir in self.dirlist:
for evt in self.evtdict[dir]:
self.stations[evt] = []
self.p_std[evt] = []
self.s_std[evt] = []
self.getPDFDict(dir, evt)
for station, pdfs in self.pdfdict.pdf_data.items():
#print station, pdfs
try:
p_std = pdfs['P'].standard_deviation()
except KeyError:
p_std = 0
try:
s_std = pdfs['S'].standard_deviation()
except KeyError:
s_std = 0
self.stations[evt].append(station)
self.p_std[evt].append(p_std)
self.s_std[evt].append(s_std)
arraylen += 1
self.makeArray(arraylen)
def makeArray(self, arraylen):
index = 0
self.p_stdarray = np.zeros(arraylen)
self.s_stdarray = np.zeros(arraylen)
for evt in self.p_std.keys():
for n in range(len(self.p_std[evt])):
self.p_stdarray[index] = self.p_std[evt][n]
self.s_stdarray[index] = self.s_std[evt][n]
index += 1
def showData(self):
#figure(p, s) = plt.pyplot.subplots(2)
#p.hist(self.p_stdarray, Bins=100, range=[min(self.p_stdarray),max(self.p_stdarray)/256])
#s.hist(self.s_stdarray, Bins=100, range=[min(self.s_stdarray),max(self.s_stdarray)/256])
#p.set_title('Histogramm der P-Wellen-Standartabweichung')
#s.set_title('Histogramm der S-Wellen-Standartabweichung')
# plt.show()
pass
def getPDFDict(self, month, evt):
self.pdfdict = PDFDictionary.from_quakeml(os.path.join(self.directory,month,evt))
def getStatistics(self):
self.p_mean = self.p_stdarray.mean()
self.p_std_std = self.p_stdarray.std()
self.p_median = np.median(self.p_stdarray)
self.s_mean = self.s_stdarray.mean()
self.s_std_std = self.s_stdarray.std()
self.s_median = np.median(self.s_stdarray)

View File

@ -294,6 +294,29 @@ class ProbabilityDensityFunction(object):
return None
return self.data[find_nearest(self.axis, value)] * self.incr
def quantile(self, prob_value, eps=0.01):
l = self.axis[0]
r = self.axis[-1]
m = (r - l) / 2
diff = prob_value - self.prob_lt_val(m)
while abs(diff) > eps:
if diff > 0:
l = m
else:
r = m
m = (r - l) / 2
diff = prob_value - self.prob_lt_val(m)
return m
pass
def quantile_distance(self, prob_value):
ql = self.quantile(prob_value)
qu = self.quantile(1 - prob_value)
return qu - ql
def plot(self, label=None):
import matplotlib.pyplot as plt