Source code for orientpy.utils

# Copyright 2019 Pascal Audet
#
# This file is part of OrientPy.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from geographiclib.geodesic import *
from obspy.signal.rotate import rotate_rt_ne, rotate_ne_rt
from obspy import UTCDateTime, read, Stream, Trace
from scipy.stats import circmean as cmean
from scipy.stats import circstd as cstd
from scipy.stats import hmean as hm
import scipy.signal as sig
import os
import math


[docs] def catclean(cat): """ This function looks for repeat events in a catalogue of earthquakes. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- cat : :class:`~obspy.core.event.Catalog` Catalogue XML object, container for Event objects. Returns ------- rep : :class:`~numpy.ndarray` 1D array of repeat events in the catalogue. """ def close(x1, x2, val): if np.abs(x1 - x2) < val: return True else: return False rep = np.array([], dtype=int) for k, event in enumerate(cat): for kk, ev in enumerate(cat): if (ev.origins[0].time - event.origins[0].time < 60.*60.) and \ close(ev.origins[0].latitude, event.origins[0].latitude, 0.8) and \ close(ev.origins[0].longitude, event.origins[0].longitude, 0.8) and \ (ev.magnitudes[0].mag < event.magnitudes[0].mag - 0.3): rep = np.append(rep, kk) return rep
[docs] def checklen(st, hrs): """ Function to check if there is enough downloaded data to run program ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- st : :class:`~obspy.core.Stream` Stream containing waveforms hrs : float Amount of time (hours) that should be available for analysis Returns ------- boolean Whether or not the check is successful """ L = len(st) for i in np.arange((L)): if (UTCDateTime(st[i].stats.endtime) - UTCDateTime(st[i].stats.starttime)) + 100. < hrs: return True if np.var(st[0].data) < 1 or np.var(st[1].data) < 1 or \ np.var(st[2].data) < 1: return True return False
[docs] def resiz(x1, x2, x3): """ Function to resize arrays to all identical shapes ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x* : :class:`~numpy.ndarray` Array for which to check length (here it's trace.data) Returns ------- x* : :class:`~numpy.ndarray` Array trimmed to shortest segment """ L = np.min(np.array([len(x1), len(x2), len(x3)])) return x1[0:L], x2[0:L], x3[0:L]
[docs] def getf(freq, A): """ Function to extract frequency in array. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- freq : float Frequency of interest (Hz) A : :class:`~numpy.ndarray` Array of frequencies (Hz) obtained from global dispersion maps Returns ------- v : float Nearest frequency in A (Hz) """ for i in np.arange((len(A))): if A[i, 0] == freq: v = A[i, 1] return v
[docs] def nv(x, v): """ Function to find nearest value in an array. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x : :class:`~numpy.ndarray` Input array v : float Value to compare within array Returns ------- x[idx] : float Nearest value to v in x """ return x[(np.abs(x - v)).argmin()]
[docs] def rms(x): """ Function to calculate root-mean-square of array. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x : :class:`~numpy.ndarray` Input array Returns ------- rms : float Root-Mean-Square value of `x` """ return np.sqrt(np.mean(np.abs(x)**2))
[docs] def mad(x): """ Function to calculate Median Absolute Deviation. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x : :class:`~numpy.ndarray` Input array Returns ------- mad : float Median Absolute deviation of `x` """ return np.median(np.abs(x - np.median(x)))
[docs] def outlier(x, lim=5.0): """ Function to remove outliers based on MAD threshold ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x : :class:`~numpy.ndarray` Input array Returns ------- x : :class:`~numpy.ndarray` Shortened array where robust standard units are lower than `lim` """ return x[np.abs(x - np.median(x))/mad(x) < lim]
[docs] def boot(x, bootnum): """ Function to calculate directional mean value from bootstrap. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- x : :class:`~numpy.ndarray` Input array bootnum : int Number of bootstrap estimates to produce Returns ------- m : :class:`~numpy.ndarray` Array of directional mean values """ m = np.zeros((bootnum)) L = len(x) for i in np.arange((bootnum)): a = np.random.choice(x, size=L, replace=True) m[i] = cmean(a, high=360) return m
[docs] def centerat(phi, m=0.): """ Function to re-center azimuth array to the mean value. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- phi : :class:`~numpy.ndarray` Input array m : float Value around which to center the azimuth Returns ------- phinew : :class:`~numpy.ndarray` Re-centered azimuth array """ phinew = np.copy(phi) if len(phi.shape) == 1: for i in np.arange((len(phi))): if phi[i] >= m+180.: phinew[i] -= 360. if phi[i] <= m-180.: phinew[i] += 360. else: for k in np.arange((phi.shape[1])): for i in np.arange((phi.shape[0])): if phi[i, k] >= m+180.: phinew[i, k] -= 360. if phi[i, k] <= m-180.: phinew[i, k] += 360. return phinew
[docs] def pathvels(lat1, lon1, lat2, lon2, map10, map15, map20, map25, map30, map35, map40): """ Overall function to get path-averaged group velocity. ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- lat1 : float Latitude of origin point (deg) lon1 : float Longitude of origin point (deg) lat2 : float Latitude of end point (deg) lon2 : float Longitude of end point (deg) map* : :class:`~numpy.ndarray` maps of Rayleigh-wave dispersion at various frequencies Returns ------- R1 : :class:`~numpy.ndarray` R1 velocity path R2 : :class:`~numpy.ndarray` R2 velocity path """ Rearth = 6371.25 circE = 2.*np.pi*Rearth # Get distance and azimuth p = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2) minor = float(p['s12']) / 1000. major = circE-minor l1 = Geodesic.WGS84.Line(lat1, lon1, p['azi1']) l2 = Geodesic.WGS84.Line(lat1, lon1, p['azi1'] - 180.) deg = 111.17 D1 = np.zeros([1, 2]) for i in np.arange((361)): b = l1.Position(deg*1000*i) p2 = np.array([b['lat2'], b['lon2']]) if i == 0: D1[0, 0] = p2[0] D1[0, 1] = p2[1] else: D1 = np.vstack((D1, p2)) bb = Geodesic.WGS84.Inverse(lat2, lon2, b['lat2'], b['lon2']) if bb['s12'] <= deg*1000.: break D2 = np.zeros([1, 2]) for i in np.arange((361)): b = l2.Position(deg*1000*i) p2 = np.array([b['lat2'], b['lon2']]) if i == 0: D2[0, 0] = p2[0] D2[0, 1] = p2[1] else: D2 = np.vstack((D2, p2)) bb = Geodesic.WGS84.Inverse(lat2, lon2, b['lat2'], b['lon2']) if bb['s12'] <= deg*1000.0: break # We now have lat and lon points along the major and minor great circles. # We calculate the group velocity of at each point, and then # find the average velocities. for k in np.arange((len(D1))): if D1[k, 1] < 0: D1[k, 1] += 360 for k in np.arange((len(D2))): if D2[k, 1] < 0: D2[k, 1] += 360 def Ray(D): U1 = np.zeros([len(D), 7]) for k in np.arange((len(D))): # Do latitude first # Get correct precision # designed to match results of Ma et al codes if np.abs(D[k, 1]) < 10: D[k, 1] = round(D[k, 1], 5) if np.abs(D[k, 1]) >= 10 and np.abs(D[k, 1]) < 100: D[k, 1] = round(D[k, 1], 4) if np.abs(D[k, 1]) > 100: D[k, 1] = round(D[k, 1], 3) if np.abs(D[k, 0]) < 10: D[k, 0] = round(D[k, 0], 5) if np.abs(D[k, 0]) >= 10: D[k, 0] = round(D[k, 0], 4) # find right index q = np.where(map10[:, 1] == nv(map10[:, 1], (D[k, 0])))[0] qq = np.where(map10[q, 0] == nv(map10[q, 0], (D[k, 1])))[0] idx = q[qq] # update path U1[k, 0] = map10[idx, 2] U1[k, 1] = map15[idx, 2] U1[k, 2] = map20[idx, 2] U1[k, 3] = map25[idx, 2] U1[k, 4] = map30[idx, 2] U1[k, 5] = map35[idx, 2] U1[k, 6] = map40[idx, 2] mhz = np.array([10, 15, 20, 25, 30, 35, 40]) return np.array((mhz, hm(U1, axis=0))).T return Ray(D1), Ray(D2)
[docs] def DLcalc(stream, Rf, LPF, HPF, epi, baz, A, winlen=10., ptype=0, zcomp='Z'): """ DORAN-LASKE calculation for one freq, one orbit of surface wave ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- stream : float Latitude of origin point (deg) lon1 : float Longitude of origin point (deg) lat2 : float Latitude of end point (deg) lon2 : float Longitude of end point (deg) map* : :class:`~numpy.ndarray` maps of Rayleigh-wave dispersion at various frequencies Returns ------- R1 : :class:`~numpy.ndarray` R1 velocity path R2 : :class:`~numpy.ndarray` R2 velocity path """ # Pre-process stream.taper(type='hann', max_percentage=0.05) stream.filter("lowpass", freq=LPF, corners=4, zerophase=True) stream.filter("highpass", freq=HPF, corners=4, zerophase=True) stream.detrend() # Window info Rvel = getf(Rf, A) # Group velocity at Rf R1window = (1.0/(Rf/1000.))*winlen arv = 1./Rvel * epi r1 = arv - R1window/2. r2 = arv + R1window/2. dt = stream[0].stats.starttime st = stream.slice(starttime=dt+r1, endtime=dt+r2) # Extract waveform data for each component try: tr1 = st.select(component='1')[0].data tr2 = st.select(component='2')[0].data except Exception: tr1 = st.select(component='N')[0].data tr2 = st.select(component='E')[0].data trZ = st.select(component=zcomp.upper())[0].data # Calculate Hilbert transform of vertical trace data trZ = np.imag(sig.hilbert(trZ)) # Ensure all data vectors are same length tr1, tr2, trZ = resiz(tr1, tr2, trZ) # Rotate through and find max normalized covariance dphi = 0.1 ang = np.arange(0., 360., dphi) cc1 = np.zeros(len(ang)) cc2 = np.zeros(len(ang)) for k, a in enumerate(ang): R, T = rotate_ne_rt(tr1, tr2, a) covmat = np.corrcoef(R, trZ) cc1[k] = covmat[0, 1] cstar = np.cov(trZ, R)/np.cov(trZ) cc2[k] = cstar[0, 1] # Get argument of maximum of cc2 ia = cc2.argmax() # Get azimuth and correct for angles above 360 phi = (baz - float(ia)*dphi) + 180. if phi < 0.: phi += 360. if phi >= 360.: phi -= 360. return phi, cc1[ia]
# Final calculation
[docs] def estimate(phi, ind): """ Function to estimate final azimuth from ADRIAN. K. DORAN and GABI LASKE, DLOPy VERSION 1.0, RELEASED APRIL 2017 Parameters ---------- phi : :class:`~numpy.ndarray` Input array of estimated azimuth ind : list List of index values that satisfy some QC condition for phi Returns ------- m : float Mean value of robust, bootstrapped estimates s : float 2 Standard deviations of robust, bootstrapped estimates """ # Get phi', where conditions are met phip = phi[ind] # Re-center at circular mean phip = centerat(phip, m=cmean(phi, high=360)) # Remove outliers phipp = outlier(phip, 5.) # bootstrap results for statistic m = boot(phipp, 5000) return cmean(m, high=360), 2*1.96*cstd(m, high=360)