[adresses #195] read_data now working correctly on QuakeML data

This commit is contained in:
Sebastian Wehling-Benatelli 2016-03-31 08:50:09 +02:00
parent bd2bad7367
commit 1fecec1696
2 changed files with 20 additions and 8 deletions

View File

@ -11,7 +11,7 @@ __version__ = _getVersionString()
__author__ = 'sebastianw' __author__ = 'sebastianw'
def readData(fn): def read_data(fn):
""" """
Reads pick data from QuakeML files named FN and returns a dictionary Reads pick data from QuakeML files named FN and returns a dictionary
containing a ProbabilityDensityFunction object for each pick. containing a ProbabilityDensityFunction object for each pick.

View File

@ -152,7 +152,7 @@ class ProbabilityDensityFunction(object):
def __nonzero__(self): def __nonzero__(self):
prec = self.precision(self.incr) prec = self.precision(self.incr)
gtzero = np.any(self.data >= 0) gtzero = np.all(self.data >= 0)
probone = bool(np.round(self.prob_gt_val(self.axis[0]), prec) == 1.) probone = bool(np.round(self.prob_gt_val(self.axis[0]), prec) == 1.)
return (gtzero and probone) return (gtzero and probone)
@ -161,7 +161,8 @@ class ProbabilityDensityFunction(object):
@staticmethod @staticmethod
def precision(incr): def precision(incr):
return int(np.ceil(np.abs(np.log10(incr)))) prec = int(np.ceil(np.abs(np.log10(incr)))) - 2
return prec if prec >= 0 else 0
@property @property
def data(self): def data(self):
@ -180,7 +181,7 @@ class ProbabilityDensityFunction(object):
self._x = np.array(x) self._x = np.array(x)
@classmethod @classmethod
def fromPick(self, lbound, barycentre, rbound, incr=0.005, decfact=0.01, def fromPick(self, lbound, barycentre, rbound, incr=0.0001, decfact=0.01,
type='gauss'): type='gauss'):
''' '''
Initialize a new ProbabilityDensityFunction object. Initialize a new ProbabilityDensityFunction object.
@ -194,7 +195,7 @@ class ProbabilityDensityFunction(object):
''' '''
# derive adequate window of definition # derive adequate window of definition
margin = 1.5 * np.max([barycentre - lbound, rbound - barycentre]) margin = 2. * np.max([barycentre - lbound, rbound - barycentre])
# find midpoint accounting also for `~obspy.UTCDateTime` object usage # find midpoint accounting also for `~obspy.UTCDateTime` object usage
try: try:
@ -206,6 +207,10 @@ class ProbabilityDensityFunction(object):
midpoint = float(rbound + float(lbound)) / 2 midpoint = float(rbound + float(lbound)) / 2
# find x0 on a grid point and sufficient npts # find x0 on a grid point and sufficient npts
was_datetime = None
if isinstance(barycentre, UTCDateTime):
barycentre = float(barycentre)
was_datetime = True
n = int(np.ceil((barycentre - midpoint) / incr)) n = int(np.ceil((barycentre - midpoint) / incr))
m = int(np.ceil((margin / incr))) m = int(np.ceil((margin / incr)))
midpoint = barycentre - n * incr midpoint = barycentre - n * incr
@ -213,6 +218,9 @@ class ProbabilityDensityFunction(object):
x0 = midpoint - margin x0 = midpoint - margin
npts = 2 * m npts = 2 * m
if was_datetime:
barycentre = UTCDateTime(barycentre)
# calculate parameter for pdf representing function # calculate parameter for pdf representing function
params = parameter[type](lbound, barycentre, rbound, decfact) params = parameter[type](lbound, barycentre, rbound, decfact)
@ -275,7 +283,7 @@ class ProbabilityDensityFunction(object):
return None return None
return self.data[find_nearest(self.axis, value)] * self.incr return self.data[find_nearest(self.axis, value)] * self.incr
def plot(self): def plot(self, label=None):
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.plot(self.axis, self.data) plt.plot(self.axis, self.data)
@ -283,9 +291,13 @@ class ProbabilityDensityFunction(object):
plt.ylabel('f(x)') plt.ylabel('f(x)')
plt.autoscale(axis='x', tight=True) plt.autoscale(axis='x', tight=True)
if self: if self:
plt.title('Probability density function') title_str = 'Probability density function '
if label:
title_str += label
title_str.strip()
else: else:
plt.title('Function not suitable as probability density function') title_str = 'Function not suitable as probability density function'
plt.title(title_str)
plt.show() plt.show()
def commonlimits(self, incr, other, max_npts=1e5): def commonlimits(self, incr, other, max_npts=1e5):