Python Pick.polarity方法代码示例

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]
def __toPick(parser, pick_el, evaluation_mode):
    """
    """
    pick = Pick()
    waveform = pick_el.xpath("waveform")[0]
    pick.waveform_id = WaveformStreamID(\
                                network_code=waveform.get("networkCode"),
                                station_code=waveform.get("stationCode"),
                                channel_code=waveform.get("channelCode"),
                                location_code=waveform.get("locationCode"))
    pick.time, pick.time_errors = __toTimeQuantity(parser, pick_el, "time")
    pick.phase_hint = parser.xpath2obj('phaseHint', pick_el)
    onset = parser.xpath2obj('onset', pick_el)
    if onset and onset.lower() in ["emergent", "impulsive", "questionable"]:
        pick.onset = onset.lower()
    # Evaluation mode of a pick is global in the SeisHub Event file format.
    pick.evaluation_mode = evaluation_mode
    # The polarity needs to be mapped.
    polarity = parser.xpath2obj('polarity', pick_el)
    pol_map_dict = {'up': 'positive', 'positive': 'positive',
                    'down': 'negative', 'negative': 'negative',
                    'undecidable': 'undecidable'}
    if polarity and polarity.lower() in pol_map_dict:
        pick.polarity = pol_map_dict[polarity.lower()]
    # Convert azimuth to backazmith
    azimuth = __toFloatQuantity(parser, pick_el, "azimuth")
    if len(azimuth) == 2 and azimuth[0] and azimuth[1]:
        # Convert to backazimuth
        pick.backazimuth = (azimuth[0] + 180.0) % 360.0
        pick.backzimuth_errors = azimuth[1]
    return pick

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]

#.........这里部分代码省略.........
            msg = ("Negative value in XX value of covariance matrix, not "
                   "setting longitude error (epicentral uncertainties will "
                   "still be set in origin uncertainty).")
            warnings.warn(msg)
        else:
            raise
    o.latitude = y
    try:
        o.latitude_errors.uncertainty = kilometer2degrees(sqrt(covariance_yy))
    except ValueError:
        if covariance_yy < 0:
            msg = ("Negative value in YY value of covariance matrix, not "
                   "setting longitude error (epicentral uncertainties will "
                   "still be set in origin uncertainty).")
            warnings.warn(msg)
        else:
            raise
    o.depth = z * 1e3  # meters!
    o.depth_errors.uncertainty = sqrt(covariance_zz) * 1e3  # meters!
    o.depth_errors.confidence_level = 68
    o.depth_type = str("from location")
    o.time = time

    ou.horizontal_uncertainty = hor_unc
    ou.min_horizontal_uncertainty = min_hor_unc
    ou.max_horizontal_uncertainty = max_hor_unc
    # values of -1 seem to be used for unset values, set to None
    for field in ("horizontal_uncertainty", "min_horizontal_uncertainty",
                  "max_horizontal_uncertainty"):
        if ou.get(field, -1) == -1:
            ou[field] = None
        else:
            ou[field] *= 1e3  # meters!
    ou.azimuth_max_horizontal_uncertainty = hor_unc_azim
    ou.preferred_description = str("uncertainty ellipse")
    ou.confidence_level = 68  # NonLinLoc in general uses 1-sigma (68%) level

    oq.standard_error = stderr
    oq.azimuthal_gap = az_gap
    oq.secondary_azimuthal_gap = sec_az_gap
    oq.used_phase_count = used_phase_count
    oq.used_station_count = used_station_count
    oq.associated_phase_count = assoc_phase_count
    oq.associated_station_count = assoc_station_count
    oq.depth_phase_count = depth_phase_count
    oq.ground_truth_level = gt_level
    oq.minimum_distance = kilometer2degrees(min_dist)
    oq.maximum_distance = kilometer2degrees(max_dist)
    oq.median_distance = kilometer2degrees(med_dist)

    # go through all phase info lines
    for line in phases_lines:
        line = line.split()
        arrival = Arrival()
        o.arrivals.append(arrival)
        station = str(line[0])
        phase = str(line[4])
        arrival.phase = phase
        arrival.distance = kilometer2degrees(float(line[21]))
        arrival.azimuth = float(line[23])
        arrival.takeoff_angle = float(line[24])
        arrival.time_residual = float(line[16])
        arrival.time_weight = float(line[17])
        pick = Pick()
        # network codes are not used by NonLinLoc, so they can not be known
        # when reading the .hyp file.. to conform with QuakeML standard set an
        # empty network code
        wid = WaveformStreamID(network_code="", station_code=station)
        # have to split this into ints for overflow to work correctly
        date, hourmin, sec = map(str, line[6:9])
        ymd = [int(date[:4]), int(date[4:6]), int(date[6:8])]
        hm = [int(hourmin[:2]), int(hourmin[2:4])]
        t = UTCDateTime(*(ymd + hm), strict=False) + float(sec)
        pick.waveform_id = wid
        pick.time = t
        pick.time_errors.uncertainty = float(line[10])
        pick.phase_hint = phase
        pick.onset = ONSETS.get(line[3].lower(), None)
        pick.polarity = POLARITIES.get(line[5].lower(), None)
        # try to determine original pick for each arrival
        for pick_ in original_picks:
            wid = pick_.waveform_id
            if station == wid.station_code and phase == pick_.phase_hint:
                pick = pick_
                break
        else:
            # warn if original picks were specified and we could not associate
            # the arrival correctly
            if original_picks:
                msg = ("Could not determine corresponding original pick for "
                       "arrival. "
                       "Falling back to pick information in NonLinLoc "
                       "hypocenter-phase file.")
                warnings.warn(msg)
        event.picks.append(pick)
        arrival.pick_id = pick.resource_id

    event.scope_resource_ids()

    return event

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]

#.........这里部分代码省略.........
    # negative values can appear on diagonal of covariance matrix due to a
    # precision problem in NLLoc implementation when location coordinates are
    # large compared to the covariances.
    o.longitude = x
    try:
        o.longitude_errors.uncertainty = kilometer2degrees(sqrt(covariance_XX))
    except ValueError:
        if covariance_XX < 0:
            msg = ("Negative value in XX value of covariance matrix, not "
                   "setting longitude error (epicentral uncertainties will "
                   "still be set in origin uncertainty).")
            warnings.warn(msg)
        else:
            raise
    o.latitude = y
    try:
        o.latitude_errors.uncertainty = kilometer2degrees(sqrt(covariance_YY))
    except ValueError:
        if covariance_YY < 0:
            msg = ("Negative value in YY value of covariance matrix, not "
                   "setting longitude error (epicentral uncertainties will "
                   "still be set in origin uncertainty).")
            warnings.warn(msg)
        else:
            raise
    o.depth = z * 1e3  # meters!
    o.depth_errors.uncertainty = sqrt(covariance_ZZ) * 1e3  # meters!
    o.depth_errors.confidence_level = 68
    o.depth_type = str("from location")
    o.time = time

    ou.horizontal_uncertainty = hor_unc
    ou.min_horizontal_uncertainty = min_hor_unc
    ou.max_horizontal_uncertainty = max_hor_unc
    # values of -1 seem to be used for unset values, set to None
    for field in ("horizontal_uncertainty", "min_horizontal_uncertainty",
                  "max_horizontal_uncertainty"):
        if ou.get(field, -1) == -1:
            ou[field] = None
        else:
            ou[field] *= 1e3  # meters!
    ou.azimuth_max_horizontal_uncertainty = hor_unc_azim
    ou.preferred_description = str("uncertainty ellipse")
    ou.confidence_level = 68  # NonLinLoc in general uses 1-sigma (68%) level

    oq.standard_error = stderr
    oq.azimuthal_gap = az_gap
    oq.secondary_azimuthal_gap = sec_az_gap
    oq.used_phase_count = used_phase_count
    oq.used_station_count = used_station_count
    oq.associated_phase_count = assoc_phase_count
    oq.associated_station_count = assoc_station_count
    oq.depth_phase_count = depth_phase_count
    oq.ground_truth_level = gt_level
    oq.minimum_distance = kilometer2degrees(min_dist)
    oq.maximum_distance = kilometer2degrees(max_dist)
    oq.median_distance = kilometer2degrees(med_dist)

    # go through all phase info lines
    for line in phases_lines:
        line = line.split()
        arrival = Arrival()
        o.arrivals.append(arrival)
        station = str(line[0])
        phase = str(line[4])
        arrival.phase = phase
        arrival.distance = kilometer2degrees(float(line[21]))
        arrival.azimuth = float(line[23])
        arrival.takeoff_angle = float(line[24])
        arrival.time_residual = float(line[16])
        arrival.time_weight = float(line[17])
        pick = Pick()
        wid = WaveformStreamID(station_code=station)
        date, hourmin, sec = map(str, line[6:9])
        t = UTCDateTime().strptime(date + hourmin, "%Y%m%d%H%M") + float(sec)
        pick.waveform_id = wid
        pick.time = t
        pick.time_errors.uncertainty = float(line[10])
        pick.phase_hint = phase
        pick.onset = ONSETS.get(line[3].lower(), None)
        pick.polarity = POLARITIES.get(line[5].lower(), None)
        # try to determine original pick for each arrival
        for pick_ in original_picks:
            wid = pick_.waveform_id
            if station == wid.station_code and phase == pick_.phase_hint:
                pick = pick_
                break
        else:
            # warn if original picks were specified and we could not associate
            # the arrival correctly
            if original_picks:
                msg = ("Could not determine corresponding original pick for "
                       "arrival. "
                       "Falling back to pick information in NonLinLoc "
                       "hypocenter-phase file.")
                warnings.warn(msg)
        event.picks.append(pick)
        arrival.pick_id = pick.resource_id

    return cat

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]
def __toPick(parser, pick_el, evaluation_mode):
    """
    """
    pick = Pick()
    pick.resource_id = ResourceIdentifier(prefix="/".join([RESOURCE_ROOT, "pick"]))

    # Raise a warnings if there is a phase delay
    phase_delay = parser.xpath2obj("phase_delay", pick_el, float)
    if phase_delay is not None:
        msg = "The pick has a phase_delay!"
        raise Exception(msg)

    waveform = pick_el.xpath("waveform")[0]
    network = waveform.get("networkCode")
    station = fix_station_name(waveform.get("stationCode"))
    # Map some station names.
    if station in STATION_DICT:
        station = STATION_DICT[station]
    if not network:
        network = NETWORK_DICT[station]

    location = waveform.get("locationCode") or ""
    channel = waveform.get("channelCode") or ""
    pick.waveform_id = WaveformStreamID(
                                network_code=network,
                                station_code=station,
                                channel_code=channel,
                                location_code=location)
    pick.time, pick.time_errors = __toTimeQuantity(parser, pick_el, "time")
    # Picks without time are not quakeml conform
    if pick.time is None:
        print "Pick has no time and is ignored: %s" % station
        return None
    pick.phase_hint = parser.xpath2obj('phaseHint', pick_el, str)
    onset = parser.xpath2obj('onset', pick_el)
    # Fixing bad and old typo ...
    if onset == "implusive":
        onset = "impulsive"
    if onset:
        pick.onset = onset.lower()
    # Evaluation mode of a pick is global in the SeisHub Event file format.
    #pick.evaluation_mode = evaluation_mode
    # The polarity needs to be mapped.
    polarity = parser.xpath2obj('polarity', pick_el)
    pol_map_dict = {'up': 'positive', 'positive': 'positive',
                    'forward': 'positive',
                    'forwards': 'positive',
                    'right': 'positive',
                    'backward': 'negative',
                    'backwards': 'negative',
                    'left': 'negative',
                    'down': 'negative', 'negative': 'negative',
                    'undecidable': 'undecidable',
                    'poorup': 'positive',
                    'poordown': 'negative'}
    if polarity:
        if polarity.lower() in pol_map_dict:
            pick.polarity = pol_map_dict[polarity.lower()]
        else:
            pick.polarity = polarity.lower()

    pick_weight = parser.xpath2obj('weight', pick_el, int)
    if pick_weight is not None:
        pick.extra = AttribDict()
        pick.extra.weight = {'value': pick_weight, 'namespace': NAMESPACE}
    return pick

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]
    def _parse_arrivals(self, event, origin, origin_res_id):
        # Skip header of arrivals
        next(self.lines)

        # Stop the loop after 2 empty lines (according to the standard).
        previous_line_empty = False

        for line in self.lines:
            line_empty = not line or line.isspace()

            if not self.event_point_separator:
                # Event are separated by two empty lines
                if line_empty and previous_line_empty:
                    break
            else:
                # Event are separated by '.'
                if line.startswith('.'):
                    break

            previous_line_empty = line_empty

            if line_empty:
                # Skip empty lines when the loop should be stopped by
                # point
                continue

            magnitude_types = []
            magnitude_values = []

            fields = self.fields['arrival']

            station = line[fields['sta']].strip()
            distance = line[fields['dist']].strip()
            event_azimuth = line[fields['ev_az']].strip()
            evaluation_mode = line[fields['picktype']].strip()
            direction = line[fields['direction']].strip()
            onset = line[fields['detchar']].strip()
            phase = line[fields['phase']].strip()
            time = line[fields['time']].strip().replace('/', '-')
            time_residual = line[fields['t_res']].strip()
            arrival_azimuth = line[fields['azim']].strip()
            azimuth_residual = line[fields['az_res']].strip()
            slowness = line[fields['slow']].strip()
            slowness_residual = line[fields['s_res']].strip()
            time_defining_flag = line[fields['t_def']].strip()
            azimuth_defining_flag = line[fields['a_def']].strip()
            slowness_defining_flag = line[fields['s_def']].strip()
            snr = line[fields['snr']].strip()
            amplitude_value = line[fields['amp']].strip()
            period = line[fields['per']].strip()
            magnitude_types.append(line[fields['mag_type_1']].strip())
            magnitude_values.append(line[fields['mag_1']].strip())
            magnitude_types.append(line[fields['mag_type_2']].strip())
            magnitude_values.append(line[fields['mag_2']].strip())
            line_id = line[fields['id']].strip()

            # Don't take pick and arrival with wrong time residual
            if '*' in time_residual:
                continue

            try:
                pick = Pick()
                pick.creation_info = self._get_creation_info()
                pick.waveform_id = WaveformStreamID()
                pick.waveform_id.station_code = station
                pick.time = UTCDateTime(time)

                network_code = self.default_network_code
                location_code = self.default_location_code
                channel_code = self.default_channel_code

                try:
                    network_code, channel = self._get_channel(station,
                                                              pick.time)
                    if channel:
                        channel_code = channel.code
                        location_code = channel.location_code
                except TypeError:
                    pass

                pick.waveform_id.network_code = network_code
                pick.waveform_id.channel_code = channel_code
                if location_code:
                    pick.waveform_id.location_code = location_code

                try:
                    ev_mode = EVALUATION_MODES[evaluation_mode]
                    pick.evaluation_mode = ev_mode
                except KeyError:
                    pass
                try:
                    pick.polarity = PICK_POLARITIES[direction]
                except KeyError:
                    pass
                try:
                    pick.onset = PICK_ONSETS[onset]
                except KeyError:
                    pass
                pick.phase_hint = phase
                try:
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]
def outputOBSPY(hp, event=None, only_fm_picks=False):
    """
    Make an Event which includes the current focal mechanism information from HASH
    
    Use the 'only_fm_picks' flag to only include the picks HASH used for the FocalMechanism.
    This flag will replace the 'picks' and 'arrivals' lists of existing events with new ones.
    
    Inputs
    -------
    hp    : hashpy.HashPype instance
    
    event : obspy.core.event.Event
    
    only_fm_picks : bool of whether to overwrite the picks/arrivals lists
    
    
    Returns
    -------
    obspy.core.event.Event
    
    Event will be new if no event was input, FocalMech added to existing event
    """
    # Returns new (or updates existing) Event with HASH solution
    n = hp.npol
    if event is None:
	event = Event(focal_mechanisms=[], picks=[], origins=[])
	origin = Origin(arrivals=[])
	origin.time = UTCDateTime(hp.tstamp)
	origin.latitude = hp.qlat
	origin.longitude = hp.qlon
	origin.depth = hp.qdep
	origin.creation_info = CreationInfo(version=hp.icusp)
	origin.resource_id = ResourceIdentifier('smi:hash/Origin/{0}'.format(hp.icusp))
	for _i in range(n):
	    p = Pick()
	    p.creation_info = CreationInfo(version=hp.arid[_i])
	    p.resource_id = ResourceIdentifier('smi:hash/Pick/{0}'.format(p.creation_info.version))
	    p.waveform_id = WaveformStreamID(network_code=hp.snet[_i], station_code=hp.sname[_i], channel_code=hp.scomp[_i])
	    if hp.p_pol[_i] > 0:
		p.polarity = 'positive'
	    else:
		p.polarity = 'negative'
	    a = Arrival()
	    a.creation_info = CreationInfo(version=hp.arid[_i])
	    a.resource_id = ResourceIdentifier('smi:hash/Arrival/{0}'.format(p.creation_info.version))
	    a.azimuth = hp.p_azi_mc[_i,0]
	    a.takeoff_angle = 180. - hp.p_the_mc[_i,0]
	    a.pick_id = p.resource_id
	    origin.arrivals.append(a)
	    event.picks.append(p)
	event.origins.append(origin)
	event.preferred_origin_id = str(origin.resource_id)
    else: # just update the changes
	origin = event.preferred_origin()
	picks = []
	arrivals = []
	for _i in range(n):
	    ind = hp.p_index[_i]
	    a = origin.arrivals[ind]
	    p = a.pick_id.getReferredObject()
	    a.takeoff_angle = hp.p_the_mc[_i,0]
	    picks.append(p)
	    arrivals.append(a)
	if only_fm_picks:
	    origin.arrivals = arrivals
	    event.picks = picks
    # Use me double couple calculator and populate planes/axes etc
    x = hp._best_quality_index
    # Put all the mechanisms into the 'focal_mechanisms' list, mark "best" as preferred
    for s in range(hp.nmult):
        dc = DoubleCouple([hp.str_avg[s], hp.dip_avg[s], hp.rak_avg[s]])
        ax = dc.axis
        focal_mech = FocalMechanism()
        focal_mech.creation_info = CreationInfo(creation_time=UTCDateTime(), author=hp.author)
        focal_mech.triggering_origin_id = origin.resource_id
        focal_mech.resource_id = ResourceIdentifier('smi:hash/FocalMechanism/{0}/{1}'.format(hp.icusp, s+1))
        focal_mech.method_id = ResourceIdentifier('HASH')
        focal_mech.nodal_planes = NodalPlanes()
        focal_mech.nodal_planes.nodal_plane_1 = NodalPlane(*dc.plane1)
        focal_mech.nodal_planes.nodal_plane_2 = NodalPlane(*dc.plane2)
        focal_mech.principal_axes = PrincipalAxes()
        focal_mech.principal_axes.t_axis = Axis(azimuth=ax['T']['azimuth'], plunge=ax['T']['dip'])
        focal_mech.principal_axes.p_axis = Axis(azimuth=ax['P']['azimuth'], plunge=ax['P']['dip'])
        focal_mech.station_polarity_count = n
        focal_mech.azimuthal_gap = hp.magap
        focal_mech.misfit = hp.mfrac[s]
        focal_mech.station_distribution_ratio = hp.stdr[s]
        focal_mech.comments.append(
            Comment(hp.qual[s], resource_id=ResourceIdentifier(str(focal_mech.resource_id) + '/comment/quality'))
            )
        #----------------------------------------
        event.focal_mechanisms.append(focal_mech)
        if s == x:
            event.preferred_focal_mechanism_id = str(focal_mech.resource_id)
    return event

# 需要导入模块: from obspy.core.event import Pick [as 别名]
# 或者: from obspy.core.event.Pick import polarity [as 别名]
    def _map_join2phase(self, db):
        """
        Return an obspy Arrival and Pick from an dict of CSS key/values
        corresponding to one record. See the 'Join' section for the implied
        database table join expected.
        
        Inputs
        ======
        db : dict of key/values of CSS fields related to the phases (see Join)

        Returns
        =======
        obspy.core.event.Pick, obspy.core.event.Arrival

        Notes
        =====
        Any object that supports the dict 'get' method can be passed as
        input, e.g. OrderedDict, custom classes, etc.
        
        Join
        ----
        assoc <- arrival <- affiliation (outer) <- schanloc [sta chan] (outer)
        
        """
        p = Pick()
        p.time = _utc(db.get('time'))
        def_net = self.agency[:2].upper()
        css_sta = db.get('sta')
        css_chan = db.get('chan')
        p.waveform_id = WaveformStreamID(
            station_code = db.get('fsta') or css_sta, 
            channel_code = db.get('fchan') or css_chan,
            network_code = db.get('snet') or def_net,
            location_code = db.get('loc'),
            )
        p.horizontal_slowness = db.get('slow')
        #p.horizontal_slowness_errors = self._create_dict(db, 'delslo')
        p.backazimuth = db.get('azimuth')
        #p.backazimuth_errors = self._create_dict(db, 'delaz')
        
        on_qual = _str(db.get('qual')).lower()
        if 'i' in on_qual:
            p.onset = "impulsive"
        elif 'e' in on_qual:
            p.onset = "emergent"
        elif 'w' in on_qual:
            p.onset = "questionable"
        else:
            p.onset =  None
        
        p.phase_hint = db.get('iphase')
        
        pol = _str(db.get('fm')).lower()
        if 'c' in pol or 'u' in pol:
            p.polarity = "positive"
        elif 'd' in pol or 'r' in pol:
            p.polarity = "negative"
        elif '.' in pol:
            p.polarity = "undecidable"
        else:
            p.polarity = None
        
        p.evaluation_mode = "automatic"
        if 'orbassoc' not in _str(db.get('auth')):
            p.evaluation_mode = "manual"
        
        p.evaluation_status = "preliminary"
        if p.evaluation_mode is "manual":
            p.evaluation_status = "reviewed"
        
        p.creation_info = CreationInfo(
            version = db.get('arid'), 
            creation_time = _utc(db.get('arrival.lddate')), 
            agency_id = self.agency, 
            author = db.get('auth'),
            )

        p.resource_id = self._rid(p)

        a = Arrival()
        a.pick_id = ResourceIdentifier(str(p.resource_id), referred_object=p)
        a.phase = db.get('phase')
        a.azimuth = db.get('esaz')
        a.distance = db.get('delta')
        a.takeoff_angle = db.get('ema')
        #a.takeoff_angle_errors = self._create_dict(db, 'emares')
        a.time_residual = db.get('timeres')
        a.horizontal_slowness_residual = db.get('slores')
        a.time_weight = db.get('wgt')
        a.earth_model_id = ResourceIdentifier(self._prefix+'/VelocityModel/'+_str(db.get('vmodel')))
        a.creation_info = CreationInfo(
            version = db.get('arid'), 
            creation_time = _utc(db.get('lddate')),
            agency_id = self.agency,
            )
        a.extra = {}
        a.extra['timedef'] = {
            'value': _str(db.get('timedef')), 
            'namespace': CSS_NAMESPACE
            }
#.........这里部分代码省略.........