example.py

# Imports
from obspy.core import UTCDateTime
from obspy.clients.fdsn import Client
from waveform_collection import gather_waveforms
from tools import waveform_tools, NIP, plotting
#%% FILTER AND WINDOW SETTINGS-------------------------------------------------------
FREQ_MIN = 1  # [Hz]
FREQ_MAX = 18
# Window length [sec]
WINDOW_LENGTH = 10
# Fraction of window overlap [0.0, 1.0)
WINDOW_OVERLAP = 0.90
# NETWORK, STATION, LOCATION, CHANNEL INFORMATION-----------------------------------
SOURCE = 'IRIS'  # FDSN source
NETWORK = 'TA'
STATION = '344A'
LOCATION = 'EP,--'  # EP infrasound for TA
CHANNEL = 'BHE,BHN,BHZ,BDF'

# Source location information
SOURCE_LAT = 32.569459  # Source latitude (Camp Minden)
SOURCE_LON = -93.352561  # Source longitude (Camp Minden)

STARTTIME = UTCDateTime('2012-10-16T04:40:30')
DURATION = 1.5*60  # [sec]
ENDTIME = STARTTIME + DURATION
# ------------------------------------------------------------------------------
#%% Get waveforms, inventory, and compute RECEIVER->SOURCE backazimuth
client = Client(SOURCE)
st = gather_waveforms(source=SOURCE, network=NETWORK, station=STATION, location=LOCATION, channel=CHANNEL,
                      starttime=STARTTIME, endtime=ENDTIME)
inv = client.get_stations(network=NETWORK, station=STATION, location=LOCATION, channel=CHANNEL, starttime=STARTTIME, endtime=ENDTIME, level="response")

# Computes source-> receiver backazimuth, removes instrument response, interpolates, checks for horizontal flipped polarity
st, backazimuth = waveform_tools.process_st(st, inv, SOURCE_LAT, SOURCE_LON)

#%% Compute Normalized Inner-Product (NIP), Mag^2 coherence (between pressure-vertical, or radial-vertical)
NIP_matrix, Cxy2, data = NIP.NIP(st, window_length=WINDOW_LENGTH, window_overlap=WINDOW_OVERLAP, backazimuth=backazimuth)
data.add_info(NIP_matrix, Cxy2, backazimuth, FREQ_MIN, FREQ_MAX)

#%% Plot waveforms coherence, and particle motion (NIP).
fig, axs = plotting.NIP_plot(st, data)