Python Origin.depth方法代码示例

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def _parse_record_ah(self, line, event):
        """
        Parses the 'additional hypocenter' record AH
        """
        date = line[2:10]
        time = line[11:20]
        # unused: hypocenter_quality = line[20]
        latitude = self._float(line[21:27])
        lat_type = line[27]
        longitude = self._float(line[29:36])
        lon_type = line[36]
        # unused: preliminary_flag = line[37]
        depth = self._float(line[38:43])
        # unused: depth_quality = line[43]
        standard_dev = self._float_unused(line[44:48])
        station_number = self._int_unused(line[48:51])
        phase_number = self._int_unused(line[51:55])
        source_code = line[56:60].strip()

        evid = event.resource_id.id.split('/')[-1]
        origin = Origin()
        res_id = '/'.join((res_id_prefix, 'origin', evid, source_code.lower()))
        origin.resource_id = ResourceIdentifier(id=res_id)
        origin.creation_info = CreationInfo(agency_id=source_code)
        origin.time = UTCDateTime(date + time)
        origin.latitude = latitude * self._coordinate_sign(lat_type)
        origin.longitude = longitude * self._coordinate_sign(lon_type)
        origin.depth = depth * 1000
        origin.depth_type = 'from location'
        origin.quality = OriginQuality()
        origin.quality.standard_error = standard_dev
        origin.quality.used_station_count = station_number
        origin.quality.used_phase_count = phase_number
        origin.origin_type = 'hypocenter'
        event.origins.append(origin)

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def _parse_record_hy(self, line):
        """
        Parses the 'hypocenter' record HY
        """
        date = line[2:10]
        time = line[11:20]
        # unused: location_quality = line[20]
        latitude = self._float(line[21:27])
        lat_type = line[27]
        longitude = self._float(line[29:36])
        lon_type = line[36]
        depth = self._float(line[38:43])
        # unused: depth_quality = line[43]
        standard_dev = self._float(line[44:48])
        station_number = self._int(line[48:51])
        # unused: version_flag = line[51]
        fe_region_number = line[52:55]
        fe_region_name = self._decode_fe_region_number(fe_region_number)
        source_code = line[55:60].strip()

        event = Event()
        # FIXME: a smarter way to define evid?
        evid = date + time
        res_id = '/'.join((res_id_prefix, 'event', evid))
        event.resource_id = ResourceIdentifier(id=res_id)
        description = EventDescription(
            type='region name',
            text=fe_region_name)
        event.event_descriptions.append(description)
        description = EventDescription(
            type='Flinn-Engdahl region',
            text=fe_region_number)
        event.event_descriptions.append(description)
        origin = Origin()
        res_id = '/'.join((res_id_prefix, 'origin', evid))
        origin.resource_id = ResourceIdentifier(id=res_id)
        origin.creation_info = CreationInfo()
        if source_code:
            origin.creation_info.agency_id = source_code
        else:
            origin.creation_info.agency_id = 'USGS-NEIC'
        res_id = '/'.join((res_id_prefix, 'earthmodel/ak135'))
        origin.earth_model_id = ResourceIdentifier(id=res_id)
        origin.time = UTCDateTime(date + time)
        origin.latitude = latitude * self._coordinate_sign(lat_type)
        origin.longitude = longitude * self._coordinate_sign(lon_type)
        origin.depth = depth * 1000
        origin.depth_type = 'from location'
        origin.quality = OriginQuality()
        origin.quality.associated_station_count = station_number
        origin.quality.standard_error = standard_dev
        # associated_phase_count can be incremented in records 'P ' and 'S '
        origin.quality.associated_phase_count = 0
        # depth_phase_count can be incremented in record 'S '
        origin.quality.depth_phase_count = 0
        origin.origin_type = 'hypocenter'
        origin.region = fe_region_name
        event.origins.append(origin)
        return event

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
 def _deserialize(self, zmap_str):
     catalog = Catalog()
     for row in zmap_str.split('\n'):
         if len(row) == 0:
             continue
         origin = Origin()
         event = Event(origins=[origin])
         event.preferred_origin_id = origin.resource_id.id
         # Begin value extraction
         columns = row.split('\t', 13)[:13]  # ignore extra columns
         values = dict(zip(_STD_ZMAP_COLUMNS + _EXT_ZMAP_COLUMNS, columns))
         # Extract origin
         origin.longitude = self._str2num(values.get('lon'))
         origin.latitude = self._str2num(values.get('lat'))
         depth = self._str2num(values.get('depth'))
         if depth is not None:
             origin.depth = depth * 1000.0
         z_err = self._str2num(values.get('z_err'))
         if z_err is not None:
             origin.depth_errors.uncertainty = z_err * 1000.0
         h_err = self._str2num(values.get('h_err'))
         if h_err is not None:
             ou = OriginUncertainty()
             ou.horizontal_uncertainty = h_err
             ou.preferred_description = 'horizontal uncertainty'
             origin.origin_uncertainty = ou
         year = self._str2num(values.get('year'))
         if year is not None:
             t_fields = ['year', 'month', 'day', 'hour', 'minute', 'second']
             comps = [self._str2num(values.get(f)) for f in t_fields]
             if year % 1 != 0:
                 origin.time = self._decyear2utc(year)
             elif any(v > 0 for v in comps[1:]):
                 # no seconds involved
                 if len(comps) < 6:
                     utc_args = [int(v) for v in comps if v is not None]
                 # we also have to handle seconds
                 else:
                     utc_args = [int(v) if v is not None else 0
                                 for v in comps[:-1]]
                     # just leave float seconds as is
                     utc_args.append(comps[-1])
                 origin.time = UTCDateTime(*utc_args)
         mag = self._str2num(values.get('mag'))
         # Extract magnitude
         if mag is not None:
             magnitude = Magnitude(mag=mag)
             m_err = self._str2num(values.get('m_err'))
             magnitude.mag_errors.uncertainty = m_err
             event.magnitudes.append(magnitude)
             event.preferred_magnitude_id = magnitude.resource_id.id
         event.scope_resource_ids()
         catalog.append(event)
     return catalog

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
 def _deserialize(self, zmap_str):
     catalog = Catalog()
     for row in zmap_str.split("\n"):
         if len(row) == 0:
             continue
         origin = Origin()
         event = Event(origins=[origin])
         event.preferred_origin_id = origin.resource_id.id
         # Begin value extraction
         columns = row.split("\t", 13)[:13]  # ignore extra columns
         values = dict(zip(_STD_ZMAP_COLUMNS + _EXT_ZMAP_COLUMNS, columns))
         # Extract origin
         origin.longitude = self._str2num(values.get("lon"))
         origin.latitude = self._str2num(values.get("lat"))
         depth = self._str2num(values.get("depth"))
         if depth is not None:
             origin.depth = depth * 1000.0
         z_err = self._str2num(values.get("z_err"))
         if z_err is not None:
             origin.depth_errors.uncertainty = z_err * 1000.0
         h_err = self._str2num(values.get("h_err"))
         if h_err is not None:
             ou = OriginUncertainty()
             ou.horizontal_uncertainty = h_err
             ou.preferred_description = "horizontal uncertainty"
             origin.origin_uncertainty = ou
         year = self._str2num(values.get("year"))
         if year is not None:
             t_fields = ["year", "month", "day", "hour", "minute", "second"]
             comps = [self._str2num(values.get(f)) for f in t_fields]
             if year % 1 != 0:
                 origin.time = self._decyear2utc(year)
             elif any(v > 0 for v in comps[1:]):
                 utc_args = [int(v) for v in comps if v is not None]
                 origin.time = UTCDateTime(*utc_args)
         mag = self._str2num(values.get("mag"))
         # Extract magnitude
         if mag is not None:
             magnitude = Magnitude(mag=mag)
             m_err = self._str2num(values.get("m_err"))
             magnitude.mag_errors.uncertainty = m_err
             event.magnitudes.append(magnitude)
             event.preferred_magnitude_id = magnitude.resource_id.id
         catalog.append(event)
     return catalog

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
def event_to_quakeml(event, filename):
    """
    Write one of those events to QuakeML.
    """
    # Create all objects.
    cat = Catalog()
    ev = Event()
    org = Origin()
    mag = Magnitude()
    fm = FocalMechanism()
    mt = MomentTensor()
    t = Tensor()
    # Link them together.
    cat.append(ev)
    ev.origins.append(org)
    ev.magnitudes.append(mag)
    ev.focal_mechanisms.append(fm)
    fm.moment_tensor = mt
    mt.tensor = t

    # Fill values
    ev.resource_id = "smi:inversion/%s" % str(event["identifier"])
    org.time = event["time"]
    org.longitude = event["longitude"]
    org.latitude = event["latitude"]
    org.depth = event["depth_in_km"] * 1000

    mag.mag = event["Mw"]
    mag.magnitude_type = "Mw"

    t.m_rr = event["Mrr"]
    t.m_tt = event["Mpp"]
    t.m_pp = event["Mtt"]
    t.m_rt = event["Mrt"]
    t.m_rp = event["Mrp"]
    t.m_tp = event["Mtp"]

    cat.write(filename, format="quakeml")

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
def __toOrigin(parser, origin_el):
    """
    Parses a given origin etree element.

    :type parser: :class:`~obspy.core.util.xmlwrapper.XMLParser`
    :param parser: Open XMLParser object.
    :type origin_el: etree.element
    :param origin_el: origin element to be parsed.
    :return: A ObsPy :class:`~obspy.core.event.Origin` object.
    """
    origin = Origin()

    # I guess setting the program used as the method id is fine.
    origin.method_id = parser.xpath2obj('program', origin_el)

    # Standard parameters.
    origin.time, origin.time_errors = \
        __toTimeQuantity(parser, origin_el, "time")
    origin.latitude, origin.latitude_errors = \
        __toFloatQuantity(parser, origin_el, "latitude")
    origin.longitude, origin.longitude_errors = \
        __toFloatQuantity(parser, origin_el, "longitude")
    origin.depth, origin.depth_errors = \
        __toFloatQuantity(parser, origin_el, "depth")

    # Figure out the depth type.
    depth_type = parser.xpath2obj("depth_type", origin_el, str)
    # Map Seishub specific depth type to the QuakeML depth type.
    if depth_type == "from location program":
        depth_type == "from location"
    origin.depth_type = "from location"

    # Earth model.
    origin.earth_model_id = parser.xpath2obj("earth_mod", origin_el, str)

    # Parse th origin uncertainty. Rather verbose but should cover all cases.
    pref_desc = parser.xpath2obj("originUncertainty/preferredDescription",
                                 origin_el, str)
    hor_uncert = parser.xpath2obj("originUncertainty/horizontalUncertainty",
                                  origin_el, float)
    min_hor_uncert = parser.xpath2obj(\
        "originUncertainty/minHorizontalUncertainty", origin_el, float)
    max_hor_uncert = parser.xpath2obj(\
        "originUncertainty/maxHorizontalUncertainty", origin_el, float)
    azi_max_hor_uncert = parser.xpath2obj(\
        "originUncertainty/azimuthMaxHorizontalUncertainty", origin_el, float)
    origin_uncert = {}
    if pref_desc:
        origin_uncert["preferred_description"] = pref_desc
    if hor_uncert:
        origin_uncert["horizontal_uncertainty"] = hor_uncert
    if min_hor_uncert:
        origin_uncert["min_horizontal_uncertainty"] = min_hor_uncert
    if max_hor_uncert:
        origin_uncert["max_horizontal_uncertainty"] = max_hor_uncert
    if azi_max_hor_uncert:
        origin_uncert["azimuth_max_horizontal_uncertainty"] = \
        azi_max_hor_uncert

    if origin_uncert:
        origin.origin_uncertainty = origin_uncert

    # Parse the OriginQuality if applicable.
    if not origin_el.xpath("originQuality"):
        return origin

    origin_quality_el = origin_el.xpath("originQuality")[0]
    origin.quality = OriginQuality()
    origin.quality.associated_phase_count = \
        parser.xpath2obj("associatedPhaseCount", origin_quality_el, int)
    # QuakeML does apparently not distinguish between P and S wave phase
    # count. Some Seishub event files do.
    p_phase_count = parser.xpath2obj("P_usedPhaseCount", origin_quality_el,
                                     int)
    s_phase_count = parser.xpath2obj("S_usedPhaseCount", origin_quality_el,
                                     int)
    # Use both in case they are set.
    if p_phase_count and s_phase_count:
        phase_count = p_phase_count + s_phase_count
        # Also add two Seishub element file specific elements.
        origin.quality.p_used_phase_count = p_phase_count
        origin.quality.s_used_phase_count = s_phase_count
    # Otherwise the total usedPhaseCount should be specified.
    else:
        phase_count = parser.xpath2obj("usedPhaseCount",
                                       origin_quality_el, int)
    origin.quality.used_phase_count = phase_count

    origin.quality.associated_station_count = \
        parser.xpath2obj("associatedStationCount", origin_quality_el, int)
    origin.quality.used_station_count = \
        parser.xpath2obj("usedStationCount", origin_quality_el, int)
    origin.quality.depth_phase_count = \
        parser.xpath2obj("depthPhaseCount", origin_quality_el, int)
    origin.quality.standard_error = \
        parser.xpath2obj("standardError", origin_quality_el, float)
    origin.quality.azimuthal_gap = \
        parser.xpath2obj("azimuthalGap", origin_quality_el, float)
    origin.quality.secondary_azimuthal_gap = \
        parser.xpath2obj("secondaryAzimuthalGap", origin_quality_el, float)
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
def read_nlloc_hyp(filename, coordinate_converter=None, picks=None, **kwargs):
    """
    Reads a NonLinLoc Hypocenter-Phase file to a
    :class:`~obspy.core.event.Catalog` object.

    .. note::

        Coordinate conversion from coordinate frame of NonLinLoc model files /
        location run to WGS84 has to be specified explicitly by the user if
        necessary.

    .. note::

        An example can be found on the :mod:`~obspy.nlloc` submodule front
        page in the documentation pages.

    :param filename: File or file-like object in text mode.
    :type coordinate_converter: func
    :param coordinate_converter: Function to convert (x, y, z)
        coordinates of NonLinLoc output to geographical coordinates and depth
        in meters (longitude, latitude, depth in kilometers).
        If left `None` NonLinLoc (x, y, z) output is left unchanged (e.g. if
        it is in geographical coordinates already like for NonLinLoc in
        global mode).
        The function should accept three arguments x, y, z and return a
        tuple of three values (lon, lat, depth in kilometers).
    :type picks: list of :class:`~obspy.core.event.Pick`
    :param picks: Original picks used to generate the NonLinLoc location.
        If provided, the output event will include the original picks and the
        arrivals in the output origin will link to them correctly (with their
        `pick_id` attribute). If not provided, the output event will include
        (the rather basic) pick information that can be reconstructed from the
        NonLinLoc hypocenter-phase file.
    :rtype: :class:`~obspy.core.event.Catalog`
    """
    if not hasattr(filename, "read"):
        # Check if it exists, otherwise assume its a string.
        try:
            with open(filename, "rt") as fh:
                data = fh.read()
        except:
            try:
                data = filename.decode()
            except:
                data = str(filename)
            data = data.strip()
    else:
        data = filename.read()
        if hasattr(data, "decode"):
            data = data.decode()

    lines = data.splitlines()

    # remember picks originally used in location, if provided
    original_picks = picks
    if original_picks is None:
        original_picks = []

    # determine indices of block start/end of the NLLOC output file
    indices_hyp = [None, None]
    indices_phases = [None, None]
    for i, line in enumerate(lines):
        if line.startswith("NLLOC "):
            indices_hyp[0] = i
        elif line.startswith("END_NLLOC"):
            indices_hyp[1] = i
        elif line.startswith("PHASE "):
            indices_phases[0] = i
        elif line.startswith("END_PHASE"):
            indices_phases[1] = i
    if any([i is None for i in indices_hyp]):
        msg = ("NLLOC HYP file seems corrupt,"
               " could not detect 'NLLOC' and 'END_NLLOC' lines.")
        raise RuntimeError(msg)
    # strip any other lines around NLLOC block
    lines = lines[indices_hyp[0]:indices_hyp[1]]

    # extract PHASES lines (if any)
    if any(indices_phases):
        if not all(indices_phases):
            msg = ("NLLOC HYP file seems corrupt, 'PHASE' block is corrupt.")
            raise RuntimeError(msg)
        i1, i2 = indices_phases
        lines, phases_lines = lines[:i1] + lines[i2 + 1:], lines[i1 + 1:i2]
    else:
        phases_lines = []

    lines = dict([line.split(None, 1) for line in lines])
    line = lines["SIGNATURE"]

    line = line.rstrip().split('"')[1]
    signature, version, date, time = line.rsplit(" ", 3)
    creation_time = UTCDateTime().strptime(date + time, str("%d%b%Y%Hh%Mm%S"))

    # maximum likelihood origin location info line
    line = lines["HYPOCENTER"]

    x, y, z = map(float, line.split()[1:7:2])

    if coordinate_converter:
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
def __toOrigin(parser, origin_el):
    """
    Parses a given origin etree element.

    :type parser: :class:`~obspy.core.util.xmlwrapper.XMLParser`
    :param parser: Open XMLParser object.
    :type origin_el: etree.element
    :param origin_el: origin element to be parsed.
    :return: A ObsPy :class:`~obspy.core.event.Origin` object.
    """
    global CURRENT_TYPE

    origin = Origin()
    origin.resource_id = ResourceIdentifier(prefix="/".join([RESOURCE_ROOT, "origin"]))

    # I guess setting the program used as the method id is fine.
    origin.method_id = "%s/location_method/%s/1" % (RESOURCE_ROOT,
        parser.xpath2obj('program', origin_el))
    if str(origin.method_id).lower().endswith("none"):
        origin.method_id = None

    # Standard parameters.
    origin.time, origin.time_errors = \
        __toTimeQuantity(parser, origin_el, "time")
    origin.latitude, origin_latitude_error = \
        __toFloatQuantity(parser, origin_el, "latitude")
    origin.longitude, origin_longitude_error = \
        __toFloatQuantity(parser, origin_el, "longitude")
    origin.depth, origin.depth_errors = \
        __toFloatQuantity(parser, origin_el, "depth")

    if origin_longitude_error:
        origin_longitude_error = origin_longitude_error["uncertainty"]
    if origin_latitude_error:
        origin_latitude_error = origin_latitude_error["uncertainty"]

    # Figure out the depth type.
    depth_type = parser.xpath2obj("depth_type", origin_el)

    # Map Seishub specific depth type to the QuakeML depth type.
    if depth_type == "from location program":
        depth_type = "from location"
    if depth_type is not None:
        origin.depth_type = depth_type

    # XXX: CHECK DEPTH ORIENTATION!!

    if CURRENT_TYPE == "seiscomp3":
        origin.depth *= 1000
        if origin.depth_errors.uncertainty:
            origin.depth_errors.uncertainty *= 1000
    else:
        # Convert to m.
        origin.depth *= -1000
        if origin.depth_errors.uncertainty:
            origin.depth_errors.uncertainty *= 1000

    # Earth model.
    earth_mod = parser.xpath2obj('earth_mod', origin_el, str)
    if earth_mod:
        earth_mod = earth_mod.split()
        earth_mod = ",".join(earth_mod)
        origin.earth_model_id = "%s/earth_model/%s/1" % (RESOURCE_ROOT,
            earth_mod)

    if (origin_latitude_error is None or origin_longitude_error is None) and \
        CURRENT_TYPE not in ["seiscomp3", "toni"]:
        print "AAAAAAAAAAAAA"
        raise Exception

    if origin_latitude_error and origin_latitude_error:
        if CURRENT_TYPE in ["baynet", "obspyck"]:
            uncert = OriginUncertainty()
            if origin_latitude_error > origin_longitude_error:
                uncert.azimuth_max_horizontal_uncertainty = 0
            else:
                uncert.azimuth_max_horizontal_uncertainty = 90
            uncert.min_horizontal_uncertainty, \
                uncert.max_horizontal_uncertainty = \
                sorted([origin_longitude_error, origin_latitude_error])
            uncert.min_horizontal_uncertainty *= 1000.0
            uncert.max_horizontal_uncertainty *= 1000.0
            uncert.preferred_description = "uncertainty ellipse"
            origin.origin_uncertainty = uncert
        elif CURRENT_TYPE == "earthworm":
            uncert = OriginUncertainty()
            uncert.horizontal_uncertainty = origin_latitude_error
            uncert.horizontal_uncertainty *= 1000.0
            uncert.preferred_description = "horizontal uncertainty"
            origin.origin_uncertainty = uncert
        elif CURRENT_TYPE in ["seiscomp3", "toni"]:
            pass
        else:
            raise Exception

    # Parse the OriginQuality if applicable.
    if not origin_el.xpath("originQuality"):
        return origin

    origin_quality_el = origin_el.xpath("originQuality")[0]
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def _parse_first_line_origin(self, line, event, magnitudes):
        """
        Parse the first line of origin data.

        :type line: str
        :param line: Line to parse.
        :type event: :class:`~obspy.core.event.event.Event`
        :param event: Event of the origin.
        :type magnitudes: list of
            :class:`~obspy.core.event.magnitude.Magnitude`
        :param magnitudes: Store magnitudes in a list to keep
            their positions.
        :rtype: :class:`~obspy.core.event.origin.Origin`,
            :class:`~obspy.core.event.resourceid.ResourceIdentifier`
        :returns: Parsed origin or None, resource identifier of the
            origin.
        """
        magnitude_types = []
        magnitude_values = []
        magnitude_station_counts = []

        fields = self.fields['line_1']

        time_origin = line[fields['time']].strip()
        time_fixed_flag = line[fields['time_fixf']].strip()
        latitude = line[fields['lat']].strip()
        longitude = line[fields['lon']].strip()
        epicenter_fixed_flag = line[fields['epicenter_fixf']].strip()
        depth = line[fields['depth']].strip()
        depth_fixed_flag = line[fields['depth_fixf']].strip()
        phase_count = line[fields['n_def']].strip()
        station_count = line[fields['n_sta']].strip()
        azimuthal_gap = line[fields['gap']].strip()
        magnitude_types.append(line[fields['mag_type_1']].strip())
        magnitude_values.append(line[fields['mag_1']].strip())
        magnitude_station_counts.append(line[fields['mag_n_sta_1']].strip())
        magnitude_types.append(line[fields['mag_type_2']].strip())
        magnitude_values.append(line[fields['mag_2']].strip())
        magnitude_station_counts.append(line[fields['mag_n_sta_2']].strip())
        magnitude_types.append(line[fields['mag_type_3']].strip())
        magnitude_values.append(line[fields['mag_3']].strip())
        magnitude_station_counts.append(line[fields['mag_n_sta_3']].strip())
        author = line[fields['author']].strip()
        origin_id = line[fields['id']].strip()

        origin = Origin()
        origin.quality = OriginQuality()

        try:
            origin.time = UTCDateTime(time_origin.replace('/', '-'))
            origin.latitude = float(latitude)
            origin.longitude = float(longitude)
        except (TypeError, ValueError):
            self._warn('Missing origin data, skipping event')
            return None, None

        origin.time_fixed = time_fixed_flag.lower() == 'f'
        origin.epicenter_fixed = epicenter_fixed_flag.lower() == 'f'

        try:
            # Convert value from km to m
            origin.depth = float(depth) * 1000
        except ValueError:
            pass
        try:
            origin.depth_type = DEPTH_TYPES[depth_fixed_flag]
        except KeyError:
            origin.depth_type = OriginDepthType('from location')
        try:
            origin.quality.used_phase_count = int(phase_count)
            origin.quality.associated_phase_count = int(phase_count)
        except ValueError:
            pass
        try:
            origin.quality.used_station_count = int(station_count)
            origin.quality.associated_station_count = int(station_count)
        except ValueError:
            pass
        try:
            origin.quality.azimuthal_gap = float(azimuthal_gap)
        except ValueError:
            pass

        self.author = author
        origin.creation_info = self._get_creation_info()

        public_id = "origin/%s" % origin_id
        origin_res_id = self._get_res_id(public_id)

        for i in range(3):
            try:
                magnitude = Magnitude()
                magnitude.creation_info = self._get_creation_info()
                magnitude.magnitude_type = magnitude_types[i]
                magnitude.mag = float(magnitude_values[i])
                magnitude.station_count = int(magnitude_station_counts[i])
                magnitude.origin_id = origin_res_id
                magnitudes.append(magnitude)
                event.magnitudes.append(magnitude)
            except ValueError:
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def _parse_record_dp(self, line, event):
        """
        Parses the 'source parameter data - primary' record Dp
        """
        source_contributor = line[2:6].strip()
        computation_type = line[6]
        exponent = self._int_zero(line[7])
        scale = math.pow(10, exponent)
        centroid_origin_time = line[8:14] + '.' + line[14]
        orig_time_stderr = line[15:17]
        if orig_time_stderr == 'FX':
            orig_time_stderr = 'Fixed'
        else:
            orig_time_stderr = \
                self._float_with_format(orig_time_stderr, '2.1', scale)
        centroid_latitude = self._float_with_format(line[17:21], '4.2')
        lat_type = line[21]
        if centroid_latitude is not None:
            centroid_latitude *= self._coordinate_sign(lat_type)
        lat_stderr = line[22:25]
        if lat_stderr == 'FX':
            lat_stderr = 'Fixed'
        else:
            lat_stderr = self._float_with_format(lat_stderr, '3.2', scale)
        centroid_longitude = self._float_with_format(line[25:30], '5.2')
        lon_type = line[30]
        if centroid_longitude is not None:
            centroid_longitude *= self._coordinate_sign(lon_type)
        lon_stderr = line[31:34]
        if lon_stderr == 'FX':
            lon_stderr = 'Fixed'
        else:
            lon_stderr = self._float_with_format(lon_stderr, '3.2', scale)
        centroid_depth = self._float_with_format(line[34:38], '4.1')
        depth_stderr = line[38:40]
        if depth_stderr == 'FX' or depth_stderr == 'BD':
            depth_stderr = 'Fixed'
        else:
            depth_stderr = self._float_with_format(depth_stderr, '2.1', scale)
        station_number = self._int_zero(line[40:43])
        component_number = self._int_zero(line[43:46])
        station_number2 = self._int_zero(line[46:48])
        component_number2 = self._int_zero(line[48:51])
        # unused: half_duration = self._float_with_format(line[51:54], '3.1')
        moment = self._float_with_format(line[54:56], '2.1')
        moment_stderr = self._float_with_format(line[56:58], '2.1')
        moment_exponent = self._int(line[58:60])
        if (moment is not None) and (moment_exponent is not None):
            moment *= math.pow(10, moment_exponent)
        if (moment_stderr is not None) and (moment_exponent is not None):
            moment_stderr *= math.pow(10, moment_exponent)

        evid = event.resource_id.id.split('/')[-1]
        # Create a new origin only if centroid time is defined:
        origin = None
        if centroid_origin_time.strip() != '.':
            origin = Origin()
            res_id = '/'.join((res_id_prefix, 'origin',
                               evid, source_contributor.lower(),
                               'mw' + computation_type.lower()))
            origin.resource_id = ResourceIdentifier(id=res_id)
            origin.creation_info = \
                CreationInfo(agency_id=source_contributor)
            date = event.origins[0].time.strftime('%Y%m%d')
            origin.time = UTCDateTime(date + centroid_origin_time)
            # Check if centroid time is on the next day:
            if origin.time < event.origins[0].time:
                origin.time += timedelta(days=1)
            self._store_uncertainty(origin.time_errors, orig_time_stderr)
            origin.latitude = centroid_latitude
            origin.longitude = centroid_longitude
            origin.depth = centroid_depth * 1000
            if lat_stderr == 'Fixed' and lon_stderr == 'Fixed':
                origin.epicenter_fixed = True
            else:
                self._store_uncertainty(origin.latitude_errors,
                                        self._lat_err_to_deg(lat_stderr))
                self._store_uncertainty(origin.longitude_errors,
                                        self._lon_err_to_deg(lon_stderr,
                                                             origin.latitude))
            if depth_stderr == 'Fixed':
                origin.depth_type = 'operator assigned'
            else:
                origin.depth_type = 'from location'
                self._store_uncertainty(origin.depth_errors,
                                        depth_stderr, scale=1000)
            quality = OriginQuality()
            quality.used_station_count = \
                station_number + station_number2
            quality.used_phase_count = \
                component_number + component_number2
            origin.quality = quality
            origin.origin_type = 'centroid'
            event.origins.append(origin)
        focal_mechanism = FocalMechanism()
        res_id = '/'.join((res_id_prefix, 'focalmechanism',
                           evid, source_contributor.lower(),
                           'mw' + computation_type.lower()))
        focal_mechanism.resource_id = ResourceIdentifier(id=res_id)
        focal_mechanism.creation_info = \
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [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 Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
def par2quakeml(Par_filename, QuakeML_filename, rotation_axis=[0.0, 1.0, 0.0],
                rotation_angle=-57.5, origin_time="2000-01-01 00:00:00.0",
                event_type="other event"):
    # initialise event
    ev = Event()

    # open and read Par file
    fid = open(Par_filename, 'r')

    fid.readline()
    fid.readline()
    fid.readline()
    fid.readline()

    lat_old = 90.0 - float(fid.readline().strip().split()[0])
    lon_old = float(fid.readline().strip().split()[0])
    depth = float(fid.readline().strip().split()[0])

    fid.readline()

    Mtt_old = float(fid.readline().strip().split()[0])
    Mpp_old = float(fid.readline().strip().split()[0])
    Mrr_old = float(fid.readline().strip().split()[0])
    Mtp_old = float(fid.readline().strip().split()[0])
    Mtr_old = float(fid.readline().strip().split()[0])
    Mpr_old = float(fid.readline().strip().split()[0])

    # rotate event into physical domain

    lat, lon = rot.rotate_lat_lon(lat_old, lon_old, rotation_axis,
                                  rotation_angle)
    Mrr, Mtt, Mpp, Mtr, Mpr, Mtp = rot.rotate_moment_tensor(
        Mrr_old, Mtt_old, Mpp_old, Mtr_old, Mpr_old, Mtp_old, lat_old, lon_old,
        rotation_axis, rotation_angle)

    # populate event origin data
    ev.event_type = event_type

    ev_origin = Origin()
    ev_origin.time = UTCDateTime(origin_time)
    ev_origin.latitude = lat
    ev_origin.longitude = lon
    ev_origin.depth = depth
    ev.origins.append(ev_origin)

    # populte event moment tensor

    ev_tensor = Tensor()
    ev_tensor.m_rr = Mrr
    ev_tensor.m_tt = Mtt
    ev_tensor.m_pp = Mpp
    ev_tensor.m_rt = Mtr
    ev_tensor.m_rp = Mpr
    ev_tensor.m_tp = Mtp

    ev_momenttensor = MomentTensor()
    ev_momenttensor.tensor = ev_tensor
    ev_momenttensor.scalar_moment = np.sqrt(Mrr ** 2 + Mtt ** 2 + Mpp ** 2 +
                                            Mtr ** 2 + Mpr ** 2 + Mtp ** 2)

    ev_focalmechanism = FocalMechanism()
    ev_focalmechanism.moment_tensor = ev_momenttensor
    ev_focalmechanism.nodal_planes = NodalPlanes().setdefault(0, 0)

    ev.focal_mechanisms.append(ev_focalmechanism)

    # populate event magnitude
    ev_magnitude = Magnitude()
    ev_magnitude.mag = 0.667 * (np.log10(ev_momenttensor.scalar_moment) - 9.1)
    ev_magnitude.magnitude_type = 'Mw'
    ev.magnitudes.append(ev_magnitude)

    # write QuakeML file
    cat = Catalog()
    cat.append(ev)
    cat.write(QuakeML_filename, format="quakeml")

    # clean up
    fid.close()

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def build(self):
        """
        Build an obspy moment tensor focal mech event

        This makes the tensor output into an Event containing:
        1) a FocalMechanism with a MomentTensor, NodalPlanes, and PrincipalAxes
        2) a Magnitude of the Mw from the Tensor

        Which is what we want for outputting QuakeML using
        the (slightly modified) obspy code.

        Input
        -----
        filehandle => open file OR str from filehandle.read()

        Output
        ------
        event => instance of Event() class as described above
        """
        p = self.parser
        event         = Event(event_type='earthquake')
        origin        = Origin()
        focal_mech    = FocalMechanism()
        nodal_planes  = NodalPlanes()
        moment_tensor = MomentTensor()
        principal_ax  = PrincipalAxes()
        magnitude     = Magnitude()
        data_used     = DataUsed()
        creation_info = CreationInfo(agency_id='NN')
        ev_mode = 'automatic'
        ev_stat = 'preliminary'
        evid = None
        orid = None
        # Parse the entire file line by line.
        for n,l in enumerate(p.line):
            if 'REVIEWED BY NSL STAFF' in l:
                ev_mode = 'manual'
                ev_stat = 'reviewed'
            if 'Event ID' in l:
                evid = p._id(n)
            if 'Origin ID' in l:
                orid = p._id(n)
            if 'Ichinose' in l:
                moment_tensor.category = 'regional'
            if re.match(r'^\d{4}\/\d{2}\/\d{2}', l):
                ev = p._event_info(n)
            if 'Depth' in l:
                derived_depth = p._depth(n)
            if 'Mw' in l:
                magnitude.mag = p._mw(n) 
                magnitude.magnitude_type = 'Mw'
            if 'Mo' in l and 'dyne' in l:
                moment_tensor.scalar_moment = p._mo(n)
            if 'Percent Double Couple' in l:
                moment_tensor.double_couple = p._percent(n)
            if 'Percent CLVD' in l:
                moment_tensor.clvd = p._percent(n)
            if 'Epsilon' in l:
                moment_tensor.variance = p._epsilon(n)
            if 'Percent Variance Reduction' in l:
                moment_tensor.variance_reduction = p._percent(n)
            if 'Major Double Couple' in l and 'strike' in p.line[n+1]:
                np = p._double_couple(n)
                nodal_planes.nodal_plane_1 = NodalPlane(*np[0])
                nodal_planes.nodal_plane_2 = NodalPlane(*np[1])
                nodal_planes.preferred_plane = 1
            if 'Spherical Coordinates' in l:
                mt = p._mt_sphere(n)
                moment_tensor.tensor = Tensor(
                    m_rr = mt['Mrr'],
                    m_tt = mt['Mtt'],
                    m_pp = mt['Mff'],
                    m_rt = mt['Mrt'],
                    m_rp = mt['Mrf'],
                    m_tp = mt['Mtf'],
                    )
            if 'Eigenvalues and eigenvectors of the Major Double Couple' in l:
                ax = p._vectors(n)
                principal_ax.t_axis = Axis(ax['T']['trend'], ax['T']['plunge'], ax['T']['ev'])
                principal_ax.p_axis = Axis(ax['P']['trend'], ax['P']['plunge'], ax['P']['ev'])
                principal_ax.n_axis = Axis(ax['N']['trend'], ax['N']['plunge'], ax['N']['ev'])
            if 'Number of Stations' in l:
                data_used.station_count = p._number_of_stations(n)
            if 'Maximum' in l and 'Gap' in l:
                focal_mech.azimuthal_gap = p._gap(n)
            if re.match(r'^Date', l):
                creation_info.creation_time = p._creation_time(n)
        # Creation Time
        creation_info.version = orid
        # Fill in magnitude values
        magnitude.evaluation_mode = ev_mode
        magnitude.evaluation_status = ev_stat
        magnitude.creation_info = creation_info.copy()
        magnitude.resource_id = self._rid(magnitude)
        # Stub origin
        origin.time = ev.get('time')
        origin.latitude = ev.get('lat')
        origin.longitude = ev.get('lon')
        origin.depth = derived_depth * 1000.
        origin.depth_type = "from moment tensor inversion"
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def get_results(self):
        cids = []
        clusters = []
        results_file = "{}/{}".format(self.hypoDD_control.control_directory,
                              self.hypoDD_control.relocated_hypocenters_output
                              )
        residuals_file = "{}/{}".format(self.hypoDD_control.control_directory,
                                        self.hypoDD_control.data_residual_output
                                        )
        with open(results_file, "r") as f:
            for line in f:
                num = line.split()
                evid = num[0]
                lat = float(num[1])
                lon = float(num[2])
                dep = 1000 * float(num[3])  # km to m
                errx = num[7]
                erry = num[8]
                errz = num[9]
                yr = int(num[10])
                mo = int(num[11])
                dy = int(num[12])
                hr = int(num[13])
                mi = int(num[14])
                sc = float(num[15])
                mag = num[16]
                nccp = num[17]
                nccs = num[18]
                nctp = num[19]
                ncts = num[20]
                rcc = num[21]
                rct = num[22]
                cid = num[23]
                if cid not in cids:
                    cids.append(cid)
                    clusters.append(Cluster())
                    clusters[-1].hypoDD_id=cid
                    clusters[-1].successful_relocation=True
                    clusters[-1].catalog=Catalog()
                    clusters[-1].event_ids=[]
                origin=Origin()
                isec = int ( math.floor( sc ))
                micsec = int ( ( sc - isec) * 1000000 )
                origin.time = UTCDateTime(yr, mo, dy, hr, mi, isec, micsec)
                origin.longitude = lon
                origin.latitude = lat
                origin.depth = dep
                origin.method_id = "hypoDD"
                # TODO (@ogalanis): Add time/location errors (when
                # appropriate. Add quality and origin_uncertainty. Add arrivals.
                event=Event()
                event.creation_info=CreationInfo()
                event.creation_info.author = __package__
                event.creation_info.version = info.__version__
                event.origins=[origin]
                event.magnitude=Magnitude()
                event.magnitude.mag=mag
                idx=cids.index(cid)
                clusters[idx].catalog.events.append(event)
                clusters[idx].event_ids.append(evid)

        if self.hypoDD_control.cid != 0 :
            my_list = []
            clusters[0].connectedness = Connectedness()
            with open(residuals_file, "r") as f:
                for line in f:
                    num = line.split()
                    evid_1 = num[2]
                    evid_2 = num[3]
                    obs_type = num[4]
                    if obs_type == "1":
                        my_list = clusters[0].connectedness.cross_corr_P
                    elif obs_type == "2":
                        my_list = clusters[0].connectedness.cross_corr_S
                    elif obs_type == "3":
                        my_list = clusters[0].connectedness.catalog_P
                    elif obs_type == "4":
                        my_list = clusters[0].connectedness.catalog_S
                    else:
                        continue
                    in_list = [x for x in my_list if (( x[0] == evid_1 and
                                                        x[1] == evid_2
                                                        ) or
                                                      ( x[0] == evid_2 and
                                                        x[1] == evid_1
                                                        ))]
                    if in_list:
                        for x in my_list:
                            if (( x[0] == evid_1 and
                                  x[1] == evid_2
                                  ) or
                                ( x[0] == evid_2 and
                                  x[1] == evid_1
                                  )):
                                x[2] += 1
                    else:
                        my_list.append([evid_1,evid_2,1])

        return clusters

# 需要导入模块: from obspy.core.event import Origin [as 别名]
# 或者: from obspy.core.event.Origin import depth [as 别名]
    def _map_join2origin(self, db):
        """
        Return an Origin instance from an dict of CSS key/values
        
        Inputs
        ======
        db : dict of key/values of CSS fields related to the origin (see Join)

        Returns
        =======
        obspy.core.event.Origin

        Notes
        =====
        Any object that supports the dict 'get' method can be passed as
        input, e.g. OrderedDict, custom classes, etc.
        
        Join
        ----
        origin <- origerr (outer)

        """ 
        #-- Basic location ------------------------------------------
        origin = Origin()
        origin.latitude = db.get('lat')
        origin.longitude = db.get('lon')
        origin.depth = _km2m(db.get('depth'))
        origin.time = _utc(db.get('time'))
        origin.extra = {}
        
        #-- Quality -------------------------------------------------
        quality = OriginQuality(
            associated_phase_count = db.get('nass'),
            used_phase_count = db.get('ndef'),
            standard_error = db.get('sdobs'),
            )
        origin.quality = quality

        #-- Solution Uncertainties ----------------------------------
        # in CSS the ellipse is projected onto the horizontal plane
        # using the covariance matrix
        uncertainty = OriginUncertainty()
        a = _km2m(db.get('smajax'))
        b = _km2m(db.get('sminax'))
        s = db.get('strike')
        dep_u = _km2m(db.get('sdepth'))
        time_u = db.get('stime')

        uncertainty.max_horizontal_uncertainty = a
        uncertainty.min_horizontal_uncertainty = b
        uncertainty.azimuth_max_horizontal_uncertainty = s
        uncertainty.horizontal_uncertainty = a
        uncertainty.preferred_description = "horizontal uncertainty"

        if db.get('conf') is not None:
            uncertainty.confidence_level = db.get('conf') * 100.  

        if uncertainty.horizontal_uncertainty is not None:
            origin.origin_uncertainty = uncertainty

        #-- Parameter Uncertainties ---------------------------------
        if all([a, b, s]):
            n, e = _get_NE_on_ellipse(a, b, s)
            lat_u = _m2deg_lat(n)
            lon_u = _m2deg_lon(e, lat=origin.latitude)
            origin.latitude_errors = {'uncertainty': lat_u} 
            origin.longitude_errors = {'uncertainty': lon_u}
        if dep_u:
            origin.depth_errors = {'uncertainty': dep_u}
        if time_u:
            origin.time_errors = {'uncertainty': time_u}

        #-- Analyst-determined Status -------------------------------
        posted_author = _str(db.get('auth'))
        mode, status = self.get_event_status(posted_author)
        origin.evaluation_mode = mode
        origin.evaluation_status = status
        
        # Save etype per origin due to schema differences...
        css_etype = _str(db.get('etype'))
        # Compatible with future patch rename "_namespace" -> "namespace"
        origin.extra['etype'] = {
            'value': css_etype, 
            'namespace': CSS_NAMESPACE
            }

        origin.creation_info = CreationInfo(
            creation_time = _utc(db.get('lddate')),
            agency_id = self.agency, 
            version = db.get('orid'),
            author = posted_author,
            )
        origin.resource_id = self._rid(origin)
        return origin