Source code for inputf

import copy
import logging
import os
from glob import glob

import numpy as num
from pyrocko import io, model

from beat import heart, utility

logger = logging.getLogger("inputf")

km = 1000.0
m = 0.000000001

[docs]def setup_stations(lats, lons, names, networks, event, rotate=True): """ Setup station objects, based on station coordinates and reference event. Parameters ---------- lats : :class:`num.ndarray` of station location latitude lons : :class:`num.ndarray` of station location longitude names : list of strings of station names networks : list of strings of network names for each station event : :class:`pyrocko.model.Event` Results ------- stations : list of :class:`pyrocko.model.Station` """ stations = [] for lat, lon, name, network in zip(lats, lons, names, networks): s = model.Station(lat=lat, lon=lon, station=name, network=network) s.set_event_relative_data(event) s.set_channels_by_name("E", "N", "Z") if rotate: p = s.guess_projections_to_rtu(out_channels=("R", "T", "Z")) s.set_channels(p[0][2]) stations.append(s) return stations
def load_matfile(datapath, **kwargs): try: return, **kwargs) except IOError: logger.warn("File %s does not exist." % datapath) return None
[docs]def load_SAR_data(datadir, names): """ Load SAR data in given directory and filenames. Returns Diff_IFG objects. """ diffgs = [] tobeloaded_names = set(copy.deepcopy(names)) for k in names: # open matlab.mat files data = load_matfile( datadir + "quad_" + k + ".mat", squeeze_me=True, struct_as_record=False ) covs = load_matfile( datadir + "CovMatrix_" + k + ".mat", squeeze_me=True, struct_as_record=False ) if data is not None and covs is not None: utmx = data["cfoc"][:, 0] utmy = data["cfoc"][:, 1] lons, lats = utility.utm_to_lonlat(utmx, utmy, 36) Lv = data["lvQT"] covariance = heart.Covariance(data=covs["Cov"]) diffgs.append( heart.DiffIFG( name=k, displacement=data["sqval"], utme=utmx, utmn=utmy, lons=lons, lats=lats, covariance=covariance, incidence=Lv.inci, heading=Lv.head, odw=data["ODW_sub"], ) ) tobeloaded_names.discard(k) else:"File %s was no SAR data?!" % datadir) names = list(tobeloaded_names) return diffgs
[docs]def load_kite_scenes(datadir, names): """ Load SAR data from the kite format. """ try: from kite import Scene from kite.scene import UserIOWarning except ImportError: raise ImportError("kite not installed! please checkout!") diffgs = [] for k in names: filepath = os.path.join(datadir, k) try:"Loading scene: %s" % k) sc = Scene.load(filepath) diffgs.append(heart.DiffIFG.from_kite_scene(sc))"Successfully imported kite scene %s" % k) except (ImportError, UserIOWarning): logger.warning("File %s not conform with kite format!" % k) return diffgs
[docs]def load_ascii_gnss_globk(filedir, filename, components=["east", "north", "up"]): """ Load ascii file columns containing: station name, Lon, Lat, ve, vn, vu, sigma_ve, sigma_vn, sigma_vu location [decimal deg] measurement unit [mm/yr] Returns ------- :class:`pyrocko.model.gnss.GNSSCampaign` """ from pyrocko.model import gnss filepath = os.path.join(filedir, filename) skiprows = 3 if os.path.exists(filepath): names = num.loadtxt(filepath, skiprows=skiprows, usecols=[12], dtype="str") d = num.loadtxt(filepath, skiprows=skiprows, usecols=range(12), dtype="float") elif len(os.path.splitext(filepath)[1]) == 0:"File %s is not an ascii text file!" % filepath) return else: raise ImportError("Did not find data under: %s" % filepath) component_to_idx_map = {"east": (2, 6), "north": (3, 7), "up": (9, 11)} # velocity_idxs = [2, 3, 9] # std_idxs = [6, 7, 11] if names.size != d.shape[0]: raise Exception("Number of stations and available data differs!") data = gnss.GNSSCampaign(name=filename) for i, name in enumerate(names): # code = str(name.split('_')[0]) # may produce duplicate sites code = str(name) gnss_station = gnss.GNSSStation( code=code, lon=float(d[i, 0]), lat=float(d[i, 1]) ) for j, comp in enumerate(components): vel_idx, std_idx = component_to_idx_map[comp] setattr( gnss_station, comp, gnss.GNSSComponent( shift=float(d[i, vel_idx] / km), sigma=float(d[i, std_idx] / km) ), ) logger.debug("Loaded station %s" % gnss_station.code) data.add_station(gnss_station) return data
def load_repsonses_from_file(projectpath): network = "" location = "" response_filename = os.path.join(projectpath, "responses.txt")"Loading responses from: %s", response_filename) responses = {} for line in open(response_filename, "r"): t = line.split() if len(t) == 8: sta, cha, instrument, lat, lon, mag, damp, period = t # please see the file format below if damp == "No_damping": damp = 0.001 lat, lon, mag, damp, period = [ float(x) for x in (lat, lon, mag, damp, period) ] responses[(network, sta, location, cha)] = (mag, damp, period) logger.debug("%s %s %s %s %s" % (sta, cha, mag, damp, period)) return responses
[docs]def load_and_blacklist_gnss( datadir, filename, blacklist, campaign=False, components=["north", "east", "up"] ): """ Load ascii GNSS data from GLOBK, apply blacklist and initialise targets. Parameters ---------- datadir: string of path to the directory filename: string filename to load blacklist: list of strings with station names to blacklist campaign : boolean if True return gnss.GNSSCampaign otherwise list of heart.GNSSCompoundComponent components : tuple of strings ('north', 'east', 'up') for displacement components to return """ gnss_campaign = load_ascii_gnss_globk(datadir, filename, components) if gnss_campaign: for station_code in blacklist: logger.debug("Removing station %s for campaign." % station_code) gnss_campaign.remove_station(station_code) station_names = [station.code for station in gnss_campaign.stations]"Loaded data of %i GNSS stations" % len(station_names)) logger.debug("Loaded GNSS stations %s" % utility.list2string(station_names)) if not campaign: return heart.GNSSCompoundComponent.from_pyrocko_gnss_campaign( gnss_campaign, components=components ) else: return gnss_campaign
[docs]def load_and_blacklist_stations(datadir, blacklist): """ Load stations from autokiwi output and apply blacklist """ stations = model.load_stations(datadir + "stations.txt") return utility.apply_station_blacklist(stations, blacklist)
def load_autokiwi(datadir, stations): return load_data_traces( datadir=datadir, stations=stations, divider="-", load_channels=["u", "r", "a"], name_prefix="reference", convert=True, ) channel_mappings = {"u": "Z", "r": "T", "a": "R", "BHE": "E", "BHN": "N", "BHZ": "Z"}
[docs]def load_obspy_data(datadir): """ Load data from the directory through obspy and convert to pyrocko objects. Parameters ---------- datadir : string absolute path to the data directory Returns ------- data_traces, stations """ import obspy from pyrocko import obspy_compat obspy_compat.plant() filenames = set(glob(datadir + "/*")) remaining_f = copy.deepcopy(filenames) print(filenames) stations = [] for f in filenames: print(f) try: inv = obspy.read_inventory(f) stations.extend(inv.to_pyrocko_stations()) remaining_f.discard(f) except TypeError: logger.debug("File %s not an inventory." % f) filenames = copy.deepcopy(remaining_f) print(filenames) data_traces = [] for f in filenames: print(f) try: stream = pyrocko_traces = stream.to_pyrocko_traces() for tr in pyrocko_traces: data_traces.append(heart.SeismicDataset.from_pyrocko_trace(tr)) remaining_f.discard(f) except TypeError: logger.debug("File %s not waveforms" % f) print(remaining_f) if len(remaining_f) > 0: logger.warning( "Could not import these files %s" % utility.list2string(list(filenames)) ) "Imported %i data_traces and %i stations" % (len(stations), len(data_traces)) ) return stations, data_traces
[docs]def load_data_traces( datadir, stations, load_channels=[], name_prefix=None, name_suffix=None, data_format="mseed", divider="-", convert=False, no_network=False, ): """ Load data traces for the given stations from datadir. """ data_trcs = [] # (r)ight transverse, (a)way radial, vertical (u)p for station in stations: if not load_channels: channels = station.channels else: channels = [model.Channel(name=cha) for cha in load_channels] for channel in channels: if no_network: trace_name = divider.join( (station.station, station.location, ) else: trace_name = divider.join( (, station.station, station.location, ) if name_suffix: trace_name = divider.join((trace_name, name_suffix)) if name_prefix: trace_name = divider.join((name_prefix, trace_name)) tracepath = os.path.join(datadir, trace_name) try: with open(tracepath): dt = io.load(tracepath, format=data_format)[0] # [nm] convert to m if convert: dt.set_ydata(dt.ydata * m) dt.set_channel( dt.set_station(station.station) dt.set_network( dt.set_location("0") # convert to BEAT seismic Dataset data_trcs.append(heart.SeismicDataset.from_pyrocko_trace(dt)) except IOError: logger.warn("Unable to open file: " + trace_name) return data_trcs
supported_channels = list(channel_mappings.values())
[docs]def rotate_traces_and_stations(datatraces, stations, event): """ Rotate traces and stations into RTZ with respect to the event. Updates channels of stations in place! Parameters --------- datatraces: list of :class:`pyrocko.trace.Trace` stations: list of :class:`pyrocko.model.Station` event: :class:`pyrocko.model.Event` Returns ------- rotated traces to RTZ """ from pyrocko import trace station2traces = utility.gather(datatraces, lambda t: t.station) trs_projected = [] for station in stations: station.set_event_relative_data(event) projections = station.guess_projections_to_rtu(out_channels=("R", "T", "Z")) try: traces = station2traces[station.station] except (KeyError): logger.warning( 'Did not find data traces for station "%s"' % stations.station ) continue ntraces = len(traces) if ntraces < 3: logger.warn( "Only found %i component(s) for station %s" % (ntraces, station.station) ) for matrix, in_channels, out_channels in projections: proc = trace.project(traces, matrix, in_channels, out_channels) for tr in proc: logger.debug("Outtrace: \n %s" % tr.__str__()) for ch in out_channels: if == station.add_channel(ch) if proc: logger.debug("Updated station: \n %s" % station.__str__()) trs_projected.extend(proc) return trs_projected