Well, I’m back again…
I’m developing a section plot, but my waveforms for the plot are not detrending and I can’t work out why. I have a lot of debugging traps in the code to switch on and off as required to try to work out where the problem is, so the code looks a bit messy no doubt.
I think I’m doing everything the same as I normally do, and I’ve tried other methods such as detrending the final stream rather than each trace, but all the results are the same. I’ve also tried all options for the detrending method.
Ultimately I want to have the section plot as a subplot beside a map, but I’ve commented that code out for now (as it has it’s own problems) while I concentrate on the detrending.
Here’s what my section plot currently looks like:
And this is the stream plot I’m producing as a debugging trap just before the section plot, showing the detrending is not working there either:
My plan is to have the local shakes in a list so I can simply comment out those I do not want to show if they don’t have data - mainly for small local quakes.
This is my code ATM:
# -*- coding: utf-8 -*-
"""
Created on Sun May 21 12:22:53 2023
@author: al72
"""
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
#from obspy.signal import filter
import matplotlib.pyplot as plt
#from matplotlib.ticker import AutoMinorLocator
#import numpy as np
#from obspy.taup import TauPyModel
#import math
#import cartopy.crs as ccrs
#import cartopy.feature as cfeature
from matplotlib.transforms import blended_transform_factory
from obspy.geodetics import gps2dist_azimuth
rs = Client('https://data.raspberryshake.org/')
def stationList(sl):
sl += ['R21C0']
sl += ['R811A']
sl += ['R9AF3']
sl += ['R571C']
#sl += ['R6D2A']
#sl += ['RF35D']
#sl += ['R7AF5']
#sl += ['R3756']
#sl += ['R6707']
#sl += ['RB18D']
#print(sl)
return sl
def freqFilter(ft):
# bandpass filter - select to suit system noise and range of quake
#ft = "'bandpass', freqmin=0.1, freqmax=0.8"
#ft = "'bandpass', freqmin=0.3, freqmax=0.8"
#ft = "'bandpass', freqmin=0.5, freqmax=2"
#ft = "'bandpass', freqmin=0.7, freqmax=2" #distant quake
#ft = "'bandpass', freqmin=0.7, freqmax=3"
#ft = "'bandpass', freqmin=0.7, freqmax=4"
#ft = "'bandpass', freqmin=0.7, freqmax=6"
#ft = "'bandpass', freqmin=0.7, freqmax=8"
#ft = "'bandpass', freqmin=0.7, freqmax=10"
#ft = "'bandpass', freqmin=1, freqmax=20"
#ft = "'bandpass', freqmin=3, freqmax=20" #use for local quakes
#ft = [0.1, 0.1, 0.8, 0.9]
#ft = [0.3, 0.3, 0.8, 0.9]
#ft = [0.5, 0.5, 2, 2.1]
ft = [0.7, 0.7, 2, 2.1] #distant quake
#ft = [0.7, 0.7, 3, 3.1]
#ft = [0.7, 0.7, 4, 4.1]
#ft = [0.7, 0.7, 6, 6.1]
#ft = [0.7, 0.7, 8, 8.1]
#ft = [1, 1, 10, 10.1]
#ft = [1, 1, 20, 20.1]
#ft = [3, 3, 20, 20.1] #use for local quakes
print(ft)
return ft
def buildStream(strm):
n = len(slist)
print(n)
#stream = Stream()
for i in range(0, n):
inv = rs.get_stations(network='AM', station=slist[i], level='RESP') # get the instrument response
#print(inv[0])
#read each epoch to find the one that's active for the event
k=0
while True:
sta = inv[0][k] #station metadata
if sta.is_active(time=eventTime):
break
k += 1
latS = sta.latitude #active station latitude
lonS = sta.longitude #active station longitude
#eleS = sta.elevation #active station elevation
print(slist[i], latS, lonS)
print(sta)
trace = rs.get_waveforms('AM', station=slist[i], location = "00", channel = 'EHZ', starttime = start, endtime = end)
trace.merge(method=0, fill_value='latest') #fill in any gaps in the data to prevent a crash
trace.detrend(type='demean') #detrend the data
#trace.filter(filt)
tr1 = trace.remove_response(inventory=inv,zero_mean=True,pre_filt=filt,output='VEL',water_level=60, plot=False) # convert to Vel
#print(tr[0])
tr1[0].stats.distance = gps2dist_azimuth(latS, lonS, latE, lonE)[0]
strm += tr1
#strm.remove_response
#strm.detrend(type='demean')
strm.plot(method='full', equal_scale=False)
return strm
eventTime = UTCDateTime(2023, 5, 13, 5, 48, 00) # (YYYY, m, d, H, M, S) **** Enter data****
latE = -32.4 # quake latitude + N -S **** Enter data****
lonE = 150.0 # quake longitude + E - W **** Enter data****
depth = 1 # quake depth, km **** Enter data****
mag = 1 # quake magnitude **** Enter data****
eventID = 'unknown' # ID for the event **** Enter data****
locE = "Mudgeeish, NSW, Australia" # location name **** Enter data****
#print(eventID)
slist = []
#filt = ""
filt = []
stationList(slist)
freqFilter(filt)
#print(slist)
print(filt)
#set up the plot
delay = 0
duration = 100
start = eventTime + delay
end = start + duration
st = Stream()
buildStream(st)
print(st)
#for t in st:
#t.stats.distance = gps2dist_azimuth(t.stats.station.latitude, t.stats.station.longitude, latE, lonE)[0]
fig = plt.figure()
st.plot(type='section', plot_dx=100e3, recordstart=delay, recordlength=duration, equal_scale=False, time_down=True, linewidth=.25, grid_linewidth=.25, show=False, fig=fig)
# Plot customization: Add station labels to offset axis
ax = fig.axes[0]
transform = blended_transform_factory(ax.transData, ax.transAxes)
for tr in st:
ax.text(tr.stats.distance / 1e3, 1.0, tr.stats.station, rotation=270, va="bottom", ha="center", transform=transform, zorder=10)
print(tr.stats.distance)
"""
fig = plt.figure(figsize=(20,14), dpi=150) # set to page size in inches
ax1 = fig.add_subplot(1,2,1)
#ax2 = fig.add_subplot(1,2,2)
st.plot(type='section', plot_dx=20e3, recordlength=100, time_down=True, linewidth=.25, grid_linewidth=.25, show=False, fig=ax1)
# Plot customization: Add station labels to offset axis
ax = ax1.axes
transform = blended_transform_factory(ax.transData, ax.transAxes)
for t in st:
ax.text(t.stats.distance / 1e3, 1.0, t.stats.station, rotation=270,
va="bottom", ha="center", transform=transform, zorder=10)
"""
plt.show()
Whatever the problem is, I’m not seeing it ATM. Probably (hopefully) something very simple.
Thanks,
Al.
