Source code for obspy.taup.ray_paths

# -*- coding: utf-8 -*-
"""
Calculations for 3D ray paths.

:author:
    Matthias Meschede
:copyright:
    The ObsPy Development Team (devs@obspy.org)
:license:
    GNU Lesser General Public License, Version 3
    (https://www.gnu.org/copyleft/lesser.html)
"""
import warnings
import numpy as np

import obspy.geodetics.base as geodetics


[docs] def get_ray_paths(inventory, catalog, phase_list=['P'], coordinate_system='XYZ', taup_model='iasp91'): """ This function returns lat, lon, depth coordinates from an event location to all stations in the inventory object :param inventory: an obspy station inventory :param catalog: an obspy event catalog :param phase_list: a list of seismic phase names that is passed to taup :param coordinate_system: can be either 'XYZ' or 'RTP'. :param taup_model: the taup model for which the greatcircle paths are computed :returns: a list of tuples ``[(gcircle, phase_name, station_label, event_timestamp, event_magnitude, event_id, origin_id), ...]``. ``gcircle`` is an array of shape ``[3, npoints]`` with the path coordinates. ``phase_name`` is the name of the seismic phase, ``station_label`` is the name of the station and network that belongs to the path. ``event_timestamp``, ``event_magnitude``, ``event_id`` and ``origin_id`` describe the event that belongs to the path. """ # GEOGRAPHICLIB is mandatory for this function if not geodetics.HAS_GEOGRAPHICLIB: raise ImportError('Geographiclib not found but required by ray path ' 'routine') stlats = [] stlons = [] stlabels = [] for network in inventory: for station in network: label_ = ".".join((network.code, station.code)) if station.latitude is None or station.longitude is None: msg = ("Station '%s' does not have latitude/longitude " "information and will not be plotted." % label_) warnings.warn(msg) continue stlats.append(station.latitude) stlons.append(station.longitude) stlabels.append(label_) # make a big list of event coordinates and names # this part should be included as a subroutine of catalog that extracts # a list of event properties. E.g. catalog.extract(['latitiude', # 'longitude', 'depth', 'mag', 'focal_mechanism') The same should be done # for an inventory with stations evlats = [] evlons = [] evdepths = [] event_ids = [] origin_ids = [] magnitudes = [] times = [] for event in catalog: if not event.origins: msg = ("Event '%s' does not have an origin and will not be " "plotted." % str(event.resource_id)) warnings.warn(msg) continue if not event.magnitudes: msg = ("Event '%s' does not have a magnitude and will not be " "plotted." % str(event.resource_id)) warnings.warn(msg) continue origin = event.preferred_origin() or event.origins[0] evlats.append(origin.latitude) evlons.append(origin.longitude) if not origin.get('depth'): # XXX do we really want to include events without depth??? origin.depth = 0. evdepths.append(origin.get('depth') * 1e-3) magnitude = event.preferred_magnitude() or event.magnitudes[0] mag = magnitude.mag event_ids.append(str(event.resource_id)) origin_ids.append(str(origin.resource_id)) magnitudes.append(mag) times.append(origin.time.timestamp) # initialize taup model if it is not provided if isinstance(taup_model, str): from obspy.taup import TauPyModel model = TauPyModel(model=taup_model) else: model = taup_model # now loop through all stations and source combinations r_earth = model.model.radius_of_planet greatcircles = [] for stlat, stlon, stlabel in zip(stlats, stlons, stlabels): for evlat, evlon, evdepth_km, time, magnitude, event_id, origin_id \ in zip(evlats, evlons, evdepths, times, magnitudes, event_ids, origin_ids): arrivals = model.get_ray_paths_geo( evdepth_km, evlat, evlon, stlat, stlon, phase_list=phase_list, resample=True) if len(arrivals) == 0: continue for arr in arrivals: radii = (r_earth - arr.path['depth']) / r_earth thetas = np.radians(90. - arr.path['lat']) phis = np.radians(arr.path['lon']) if coordinate_system == 'RTP': gcircle = np.array([radii, thetas, phis]) if coordinate_system == 'XYZ': gcircle = np.array([radii * np.sin(thetas) * np.cos(phis), radii * np.sin(thetas) * np.sin(phis), radii * np.cos(thetas)]) greatcircles.append((gcircle, arr.name, stlabel, time, magnitude, event_id, origin_id)) return greatcircles