Python Magnitude.origin_id方法代码示例

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
    def _parse_record_ae(self, line, event):
        """
        Parses the 'additional hypocenter error and magnitude record' AE
        """
        orig_time_stderr = self._float_unused(line[2:7])
        latitude_stderr = self._float_unused(line[8:14])
        longitude_stderr = self._float_unused(line[15:21])
        depth_stderr = self._float_unused(line[22:27])
        gap = self._float_unused(line[28:33])
        mag1 = self._float(line[33:36])
        mag1_type = line[36:38]
        mag2 = self._float(line[43:46])
        mag2_type = line[46:48]

        evid = event.resource_id.id.split('/')[-1]
        # this record is to be associated to the latest origin
        origin = event.origins[-1]
        self._store_uncertainty(origin.time_errors, orig_time_stderr)
        self._store_uncertainty(origin.latitude_errors,
                                self._lat_err_to_deg(latitude_stderr))
        self._store_uncertainty(origin.longitude_errors,
                                self._lon_err_to_deg(longitude_stderr,
                                                     origin.latitude))
        self._store_uncertainty(origin.depth_errors, depth_stderr, scale=1000)
        origin.quality.azimuthal_gap = gap
        if mag1 > 0:
            mag = Magnitude()
            mag1_id = mag1_type.lower()
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, mag1_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(
                agency_id=origin.creation_info.agency_id)
            mag.mag = mag1
            mag.magnitude_type = mag1_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if mag2 > 0:
            mag = Magnitude()
            mag2_id = mag2_type.lower()
            if mag2_id == mag1_id:
                mag2_id += '2'
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, mag2_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(
                agency_id=origin.creation_info.agency_id)
            mag.mag = mag2
            mag.magnitude_type = mag2_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
def __toMagnitude(parser, magnitude_el, origin):
    """
    Parses a given magnitude etree element.

    :type parser: :class:`~obspy.core.util.xmlwrapper.XMLParser`
    :param parser: Open XMLParser object.
    :type magnitude_el: etree.element
    :param magnitude_el: magnitude element to be parsed.
    :return: A ObsPy :class:`~obspy.core.event.Magnitude` object.
    """
    global CURRENT_TYPE
    mag = Magnitude()
    mag.resource_id = ResourceIdentifier(prefix="/".join([RESOURCE_ROOT, "magnitude"]))
    mag.origin_id = origin.resource_id
    mag.mag, mag.mag_errors = __toFloatQuantity(parser, magnitude_el, "mag")
    # obspyck used to write variance (instead of std) in magnitude error fields
    if CURRENT_TYPE == "obspyck":
        if mag.mag_errors.uncertainty is not None:
            mag.mag_errors.uncertainty = math.sqrt(mag.mag_errors.uncertainty)
    mag.magnitude_type = parser.xpath2obj("type", magnitude_el)
    mag.station_count = parser.xpath2obj("stationCount", magnitude_el, int)
    mag.method_id = "%s/magnitude_method/%s/1" % (RESOURCE_ROOT,
        parser.xpath2obj('program', magnitude_el))
    if str(mag.method_id).lower().endswith("none"):
        mag.method_id = None

    return mag

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
    def _on_file_save(self):
        """
        Creates a new obspy.core.event.Magnitude object and writes the moment
        magnitude to it.
        """
        # Get the save filename.
        filename = QtGui.QFileDialog.getSaveFileName(caption="Save as...")
        filename = os.path.abspath(str(filename))
        mag = Magnitude()
        mag.mag = self.final_result["moment_magnitude"]
        mag.magnitude_type = "Mw"
        mag.station_count = self.final_result["station_count"]
        mag.evaluation_mode = "manual"
        # Link to the used origin.
        mag.origin_id = self.current_state["event"].origins[0].resource_id
        mag.method_id = "Magnitude picker Krischer"
        # XXX: Potentially change once this program gets more stable.
        mag.evaluation_status = "preliminary"
        # Write the other results as Comments.
        mag.comments.append( \
            Comment("Seismic moment in Nm: %g" % \
            self.final_result["seismic_moment"]))
        mag.comments.append( \
            Comment("Circular source radius in m: %.2f" % \
            self.final_result["source_radius"]))
        mag.comments.append( \
            Comment("Stress drop in Pa: %.2f" % \
            self.final_result["stress_drop"]))
        mag.comments.append( \
                Comment("Very rough Q estimation: %.1f" % \
            self.final_result["quality_factor"]))

        event = copy.deepcopy(self.current_state["event"])
        event.magnitudes.append(mag)
        cat = Catalog()
        cat.events.append(event)
        cat.write(filename, format="quakeml")

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [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 Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
    def _parse_record_e(self, line, event):
        """
        Parses the 'error and magnitude' record E
        """
        orig_time_stderr = self._float(line[2:7])
        latitude_stderr = self._float(line[8:14])
        longitude_stderr = self._float(line[15:21])
        depth_stderr = self._float(line[22:27])
        mb_mag = self._float(line[28:31])
        mb_nsta = self._int(line[32:35])
        ms_mag = self._float(line[36:39])
        ms_nsta = self._int(line[39:42])
        mag1 = self._float(line[42:45])
        mag1_type = line[45:47]
        mag1_source_code = line[47:51].strip()
        mag2 = self._float(line[51:54])
        mag2_type = line[54:56]
        mag2_source_code = line[56:60].strip()

        evid = event.resource_id.id.split('/')[-1]
        origin = event.origins[0]
        self._store_uncertainty(origin.time_errors, orig_time_stderr)
        self._store_uncertainty(origin.latitude_errors,
                                self._lat_err_to_deg(latitude_stderr))
        self._store_uncertainty(origin.longitude_errors,
                                self._lon_err_to_deg(longitude_stderr,
                                                     origin.latitude))
        self._store_uncertainty(origin.depth_errors, depth_stderr, scale=1000)
        if mb_mag is not None:
            mag = Magnitude()
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, 'mb'))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id='USGS-NEIC')
            mag.mag = mb_mag
            mag.magnitude_type = 'Mb'
            mag.station_count = mb_nsta
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if ms_mag is not None:
            mag = Magnitude()
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, 'ms'))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id='USGS-NEIC')
            mag.mag = ms_mag
            mag.magnitude_type = 'Ms'
            mag.station_count = ms_nsta
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if mag1 is not None:
            mag = Magnitude()
            mag1_id = mag1_type.lower()
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, mag1_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id=mag1_source_code)
            mag.mag = mag1
            mag.magnitude_type = mag1_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if mag2 is not None:
            mag = Magnitude()
            mag2_id = mag2_type.lower()
            if mag2_id == mag1_id:
                mag2_id += '2'
            res_id = '/'.join((res_id_prefix, 'magnitude', evid, mag2_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id=mag2_source_code)
            mag.mag = mag2
            mag.magnitude_type = mag2_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)

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

#.........这里部分代码省略.........
            corner_freqs = []
            for trace in stream:
                # Get the index of the pick.
                pick_index = int(round((pick.time - trace.stats.starttime) / \
                    trace.stats.delta))
                # Choose date window 0.5 seconds before and 1 second after pick.
                data_window = trace.data[pick_index - \
                    int(TIME_BEFORE_PICK * trace.stats.sampling_rate): \
                    pick_index + int(TIME_AFTER_PICK * trace.stats.sampling_rate)]
                # Calculate the spectrum.
                spec, freq = mtspec.mtspec(data_window, trace.stats.delta, 2)
                try:
                    fit = fit_spectrum(spec, freq, pick.time - origin_time,
                            spec.max(), 10.0)
                except:
                    continue
                if fit is None:
                    continue
                Omega_0, f_c, err, _ = fit
                Omega_0 = np.sqrt(Omega_0)
                omegas.append(Omega_0)
                corner_freqs.append(f_c)
            M_0 = 4.0 * np.pi * DENSITY * velocity ** 3 * distance * \
                np.sqrt(omegas[0] ** 2 + omegas[1] ** 2 + omegas[2] ** 2) / \
                radiation_pattern
            r = 3 * k * V_S / sum(corner_freqs)
            moments.append(M_0)
            source_radii.append(r)
            corner_frequencies.extend(corner_freqs)
        if not len(moments):
            print "No moments could be calculated for event %s" % \
                event.resource_id.resource_id
            continue

        # Calculate the seismic moment via basic statistics.
        moments = np.array(moments)
        moment = moments.mean()
        moment_std = moments.std()

        corner_frequencies = np.array(corner_frequencies)
        corner_frequency = corner_frequencies.mean()
        corner_frequency_std = corner_frequencies.std()

        # Calculate the source radius.
        source_radii = np.array(source_radii)
        source_radius = source_radii.mean()
        source_radius_std = source_radii.std()

        # Calculate the stress drop of the event based on the average moment and
        # source radii.
        stress_drop = (7 * moment) / (16 * source_radius ** 3)
        stress_drop_std = np.sqrt((stress_drop ** 2) * \
            (((moment_std ** 2) / (moment ** 2)) + \
            (9 * source_radius * source_radius_std ** 2)))
        if source_radius > 0 and source_radius_std < source_radius:
            print "Source radius:", source_radius, " Std:", source_radius_std
            print "Stress drop:", stress_drop / 1E5, " Std:", stress_drop_std / 1E5

        Mw = 2.0 / 3.0 * (np.log10(moment) - 9.1)
        Mw_std = 2.0 / 3.0 * moment_std / (moment * np.log(10))
        Mws_std.append(Mw_std)
        Mws.append(Mw)
        Mls.append(local_magnitude)
        calc_diff = abs(Mw - local_magnitude)
        Mw = ("%.3f" % Mw).rjust(7)
        Ml = ("%.3f" % local_magnitude).rjust(7)
        diff = ("%.3e" % calc_diff).rjust(7)

        ret_string = colorama.Fore.GREEN + \
            "For event %s: Ml=%s | Mw=%s | " % (event.resource_id.resource_id,
            Ml, Mw)
        if calc_diff >= 1.0:
            ret_string += colorama.Fore.RED
        ret_string += "Diff=%s" % diff
        ret_string += colorama.Fore.GREEN
        ret_string += " | Determined at %i stations" % len(moments)
        ret_string += colorama.Style.RESET_ALL
        print ret_string

        mag = Magnitude()
        mag.mag = Mw
        mag.mag_errors.uncertainty = Mw_std
        mag.magnitude_type = "Mw"
        mag.origin_id = event.origins[0].resource_id
        mag.method_id = "smi:com.github/krischer/moment_magnitude_calculator/automatic/1"
        mag.station_count = len(moments)
        mag.evaluation_mode = "automatic"
        mag.evaluation_status = "preliminary"
        mag.comments.append(Comment( \
            "Seismic Moment=%e Nm; standard deviation=%e" % (moment,
            moment_std)))
        mag.comments.append(Comment("Custom fit to Boatwright spectrum"))
        if source_radius > 0 and source_radius_std < source_radius:
            mag.comments.append(Comment( \
                "Source radius=%.2fm; standard deviation=%.2f" % (source_radius,
                source_radius_std)))
        event.magnitudes.append(mag)

    print "Writing output file..."
    cat.write(output_file, format="quakeml")

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
def iris2quakeml(url, output_folder=None):
    if not "/spudservice/" in url:
        url = url.replace("/spud/", "/spudservice/")
        if url.endswith("/"):
            url += "quakeml"
        else:
            url += "/quakeml"
    print "Downloading %s..." % url
    r = requests.get(url)
    if r.status_code != 200:
        msg = "Error Downloading file!"
        raise Exception(msg)

    # For some reason the quakeml file is escaped HTML.
    h = HTMLParser.HTMLParser()

    data = h.unescape(r.content)

    # Replace some XML tags.
    data = data.replace("long-period body waves", "body waves")
    data = data.replace("intermediate-period surface waves", "surface waves")
    data = data.replace("long-period mantle waves", "mantle waves")

    data = data.replace("<html><body><pre>", "")
    data = data.replace("</pre></body></html>", "")

    # Change the resource identifiers. Colons are not allowed in QuakeML.
    pattern = r"(\d{4})-(\d{2})-(\d{2})T(\d{2}):(\d{2}):(\d{2})\.(\d{6})"
    data = re.sub(pattern, r"\1-\2-\3T\4-\5-\6.\7", data)

    data = StringIO(data)

    try:
        cat = readEvents(data)
    except:
        msg = "Could not read downloaded event data"
        raise ValueError(msg)

    # Parse the event, and use only one origin, magnitude and focal mechanism.
    # Only the first event is used. Should not be a problem for the chosen
    # global cmt application.
    ev = cat[0]

    if ev.preferred_origin():
        ev.origins = [ev.preferred_origin()]
    else:
        ev.origins = [ev.origins[0]]
    if ev.preferred_focal_mechanism():
        ev.focal_mechanisms = [ev.preferred_focal_mechanism()]
    else:
        ev.focal_mechanisms = [ev.focal_mechanisms[0]]

    try:
        mt = ev.focal_mechanisms[0].moment_tensor
    except:
        msg = "No moment tensor found in file."
        raise ValueError
    seismic_moment_in_dyn_cm = mt.scalar_moment
    if not seismic_moment_in_dyn_cm:
        msg = "No scalar moment found in file."
        raise ValueError(msg)

    # Create a new magnitude object with the moment magnitude calculated from
    # the given seismic moment.
    mag = Magnitude()
    mag.magnitude_type = "Mw"
    mag.origin_id = ev.origins[0].resource_id
    # This is the formula given on the GCMT homepage.
    mag.mag = (2.0 / 3.0) * (math.log10(seismic_moment_in_dyn_cm) - 16.1)
    mag.resource_id = ev.origins[0].resource_id.resource_id.replace("Origin",
        "Magnitude")
    ev.magnitudes = [mag]
    ev.preferred_magnitude_id = mag.resource_id

    # Convert the depth to meters.
    org = ev.origins[0]
    org.depth *= 1000.0
    if org.depth_errors.uncertainty:
        org.depth_errors.uncertainty *= 1000.0

    # Ugly asserts -- this is just a simple script.
    assert(len(ev.magnitudes) == 1)
    assert(len(ev.origins) == 1)
    assert(len(ev.focal_mechanisms) == 1)

    # All values given in the QuakeML file are given in dyne * cm. Convert them
    # to N * m.
    for key, value in mt.tensor.iteritems():
        if key.startswith("m_") and len(key) == 4:
            mt.tensor[key] /= 1E7
        if key.endswith("_errors") and hasattr(value, "uncertainty"):
            mt.tensor[key].uncertainty /= 1E7
    mt.scalar_moment /= 1E7
    if mt.scalar_moment_errors.uncertainty:
        mt.scalar_moment_errors.uncertainty /= 1E7
    p_axes = ev.focal_mechanisms[0].principal_axes
    for ax in [p_axes.t_axis, p_axes.p_axis, p_axes.n_axis]:
        if ax is None or not ax.length:
            continue
        ax.length /= 1E7
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.event import Magnitude [as 别名]
# 或者: from obspy.core.event.Magnitude import origin_id [as 别名]
    def _parseRecordE(self, line, event):
        """
        Parses the 'error and magnitude' record E
        """
        orig_time_stderr = self._float(line[2:7])
        latitude_stderr = self._float(line[8:14])
        longitude_stderr = self._float(line[15:21])
        depth_stderr = self._float(line[22:27])
        mb_mag = self._float(line[28:31])
        mb_nsta = self._int(line[32:35])
        Ms_mag = self._float(line[36:39])
        Ms_nsta = self._int(line[39:42])
        mag1 = self._float(line[42:45])
        mag1_type = line[45:47]
        mag1_source_code = line[47:51].strip()
        mag2 = self._float(line[51:54])
        mag2_type = line[54:56]
        mag2_source_code = line[56:60].strip()

        evid = event.resource_id.id.split("/")[-1]
        origin = event.origins[0]
        self._storeUncertainty(origin.time_errors, orig_time_stderr)
        self._storeUncertainty(origin.latitude_errors, self._latErrToDeg(latitude_stderr))
        self._storeUncertainty(origin.longitude_errors, self._lonErrToDeg(longitude_stderr, origin.latitude))
        self._storeUncertainty(origin.depth_errors, depth_stderr, scale=1000)
        if mb_mag is not None:
            mag = Magnitude()
            res_id = "/".join((res_id_prefix, "magnitude", evid, "mb"))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id="USGS-NEIC")
            mag.mag = mb_mag
            mag.magnitude_type = "Mb"
            mag.station_count = mb_nsta
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if Ms_mag is not None:
            mag = Magnitude()
            res_id = "/".join((res_id_prefix, "magnitude", evid, "ms"))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id="USGS-NEIC")
            mag.mag = Ms_mag
            mag.magnitude_type = "Ms"
            mag.station_count = Ms_nsta
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if mag1 is not None:
            mag = Magnitude()
            mag1_id = mag1_type.lower()
            res_id = "/".join((res_id_prefix, "magnitude", evid, mag1_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id=mag1_source_code)
            mag.mag = mag1
            mag.magnitude_type = mag1_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)
        if mag2 is not None:
            mag = Magnitude()
            mag2_id = mag2_type.lower()
            if mag2_id == mag1_id:
                mag2_id += "2"
            res_id = "/".join((res_id_prefix, "magnitude", evid, mag2_id))
            mag.resource_id = ResourceIdentifier(id=res_id)
            mag.creation_info = CreationInfo(agency_id=mag2_source_code)
            mag.mag = mag2
            mag.magnitude_type = mag2_type
            mag.origin_id = origin.resource_id
            event.magnitudes.append(mag)