Kasper D. Fischer
af7c6642ff
placeholder is "VVVVV" use patchVersion.sh to replace version for new release. Do not use this script in develop branch.
406 lines
14 KiB
Python
Executable File
406 lines
14 KiB
Python
Executable File
#!/usr/bin/env python
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
Get waveform data from FDSN web service and create a fancy plot
|
|
This programme runs as a script or as a WSGI application.
|
|
|
|
:version VVVVV
|
|
|
|
:license
|
|
Copyright 2020 Kasper Fischer <kasper.fischer@ruhr-uni-bochum.de>
|
|
|
|
This program is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation; either version 2 of the License, or
|
|
(at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with this program; if not, write to the Free Software
|
|
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|
MA 02110-1301, USA.
|
|
"""
|
|
|
|
|
|
def utc2local(utctime, timezone='Europe/Berlin'):
|
|
"""
|
|
utc2local(utctime, timezone='Europe/Berlin')
|
|
converts the UTCDateTime object utctime into a
|
|
datetime object of the given timezone
|
|
"""
|
|
|
|
import pytz
|
|
|
|
utc = pytz.utc
|
|
local = pytz.timezone(timezone)
|
|
utctime = utc.localize(utctime.datetime)
|
|
return local.normalize(utctime.astimezone(local))
|
|
|
|
|
|
def fancy_plot(st, wsgi=False, img_format='png', color=True):
|
|
"""
|
|
creating fancy plot from ObsPy stream st
|
|
returns cStringIO object if wsgi == True
|
|
:type st: ObsPy Stream object
|
|
:type wsgi: bool
|
|
:type img_format: str
|
|
:type color: bool
|
|
"""
|
|
|
|
import matplotlib as mpl
|
|
from matplotlib.dates import date2num, AutoDateLocator, AutoDateFormatter
|
|
|
|
if wsgi:
|
|
try:
|
|
import cStringIO as StringIO
|
|
except ImportError:
|
|
import StringIO as StringIO
|
|
import matplotlib.pyplot as plt
|
|
from numpy import arange
|
|
|
|
# setup figure parameters
|
|
if wsgi:
|
|
mpl.rcParams['figure.figsize'] = [12.0, 6.75] # 16:9 aspect ratio
|
|
else:
|
|
mpl.rcParams['figure.figsize'] = [16.0, 9.0] # 16:9 aspect ratio
|
|
|
|
# get starttime, endtime, localtime and timezone
|
|
tr = st[0]
|
|
stime = tr.stats.starttime
|
|
etime = tr.stats.endtime
|
|
localtime = utc2local(stime)
|
|
if localtime.tzname() == u'CET':
|
|
tzname = u'MEZ'
|
|
else:
|
|
tzname = u'MESZ'
|
|
tz = localtime.timetz().tzinfo
|
|
|
|
# create proper time vector
|
|
d = arange(date2num(stime), date2num(etime), (date2num(etime) - date2num(stime)) / tr.stats.npts)
|
|
|
|
# setup time axis decorator
|
|
locator = AutoDateLocator(interval_multiples=True, minticks=5, maxticks=8, tz=tz)
|
|
minor_locator = AutoDateLocator(interval_multiples=True, minticks=9, maxticks=70, tz=tz)
|
|
formatter = AutoDateFormatter(locator, tz=tz)
|
|
formatter.scaled[1. / (24. * 60.)] = '%H:%M:%S'
|
|
formatter.scaled[1. / (24. * 60. * 60. * 10.)] = '%H:%M:%S.%f'
|
|
|
|
# draw figure
|
|
if color:
|
|
trace_color = {'Z': 'r', 'N': 'b', 'E': 'g'}
|
|
else:
|
|
trace_color = {'Z': 'k', 'N': 'k', 'E': 'k'}
|
|
|
|
fig = plt.figure()
|
|
|
|
ax1 = fig.add_subplot(311)
|
|
ax1.xaxis.set_major_locator(locator)
|
|
ax1.xaxis.set_major_formatter(formatter)
|
|
ax1.plot_date(d, st.select(component='Z')[0].data * 1000., trace_color['Z'], label=u'Z', tz=tz)
|
|
plt.title(u'Station %s' % st[0].stats.station, fontsize=24)
|
|
plt.legend(loc='upper right')
|
|
ax1.grid(b=True)
|
|
|
|
ax2 = fig.add_subplot(312, sharex=ax1, sharey=ax1)
|
|
ax2.plot_date(d, st.select(component='N')[0].data * 1000., trace_color['N'], label=u'N', tz=tz)
|
|
ax2.grid(b=True)
|
|
plt.ylabel(u'Geschwindigkeit [mm/s]', fontsize=16)
|
|
plt.legend(loc='upper right')
|
|
|
|
ax3 = fig.add_subplot(313, sharex=ax1, sharey=ax1)
|
|
ax3.plot_date(d, st.select(component='E')[0].data * 1000., trace_color['E'], label=u'E', tz=tz)
|
|
ax3.grid(b=True)
|
|
plt.legend(loc='upper right')
|
|
|
|
plt.xlabel(u'%s (%s)' % (localtime.strftime('%d.%m.%Y'), tzname), fontsize=16)
|
|
ax1.minorticks_on()
|
|
ax1.xaxis.set_minor_locator(minor_locator)
|
|
|
|
if wsgi:
|
|
buf = StringIO.StringIO()
|
|
plt.savefig(buf, format=img_format, facecolor="lightgray")
|
|
return buf
|
|
else:
|
|
if cla['filename']:
|
|
plt.savefig(cla['filename'], dpi=300, transparent=True)
|
|
else:
|
|
plt.show()
|
|
|
|
|
|
def trace_dayplot(st, deltat = None,
|
|
ftype='none', fmin=1.0, fmax=7.0,
|
|
col=('b', 'r', 'g'), interval=20, outpattern='',
|
|
wsgi=False):
|
|
"""
|
|
|
|
:type st: object
|
|
:type ftype: str
|
|
:type fmin: float
|
|
:type fmax: float
|
|
:type col: tuple
|
|
:type interval: float
|
|
:type outpattern: str
|
|
:type wsgi: bool
|
|
"""
|
|
|
|
if wsgi:
|
|
try:
|
|
import cStringIO as StringIO
|
|
except ImportError:
|
|
import StringIO as StringIO
|
|
|
|
# filter
|
|
if (ftype == 'bandpass') or (ftype == 'bandstop'):
|
|
st.filter(ftype, freqmin=fmin, freqmax=fmax)
|
|
elif (ftype == 'lowpass') or (ftype == 'highpass'):
|
|
st.filter(ftype, freq=fmin)
|
|
st.merge()
|
|
|
|
stime = st[0].stats.starttime
|
|
if deltat:
|
|
etime = stime + deltat
|
|
else:
|
|
etime = st[0].stats.endtime
|
|
|
|
# plot
|
|
for i in range(0, len(st)):
|
|
if wsgi:
|
|
buf = StringIO.StringIO()
|
|
st[i].plot(outfile=buf, format=outpattern, type='dayplot', color=col, interval=interval,
|
|
starttime=stime, endtime=etime)
|
|
return buf
|
|
else:
|
|
if outpattern != '':
|
|
imagefile = outpattern.format(st[i].stats.network,
|
|
st[i].stats.station,
|
|
st[i].stats.channel,
|
|
st[i].stats.starttime.year,
|
|
st[i].stats.starttime.julday)
|
|
st[i].plot(type='dayplot', color=col, interval=interval,
|
|
outfile=imagefile)
|
|
else:
|
|
st[i].plot(type='dayplot', color=col, interval=interval)
|
|
|
|
|
|
def get_fdsn(bulk_request, output='VEL', base_url='https://ariadne.geophysik.ruhr-uni-bochum.de'):
|
|
"""
|
|
Fetches waveform data from FDSN web service and returns
|
|
instrument corrected seismogram of type given in parameter output.
|
|
Acceptable values for output are "VEL", "DISP" and "ACC" (see ObsPy documentation)
|
|
|
|
:rtype : object
|
|
:param bulk_request: list
|
|
:param output: str
|
|
:param base_url: str
|
|
:return: ObsPy Stream() object
|
|
"""
|
|
import warnings
|
|
from obspy.fdsn import Client
|
|
from obspy.fdsn.header import FDSNException
|
|
|
|
try:
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore")
|
|
client = Client(base_url=base_url, debug=False)
|
|
st = client.get_waveforms_bulk(bulk_request, attach_response=True)
|
|
|
|
st.merge()
|
|
stime = st[0].stats.starttime
|
|
etime = st[0].stats.endtime
|
|
for trace in st.traces:
|
|
stime = max(stime, trace.stats.starttime)
|
|
etime = min(etime, trace.stats.endtime)
|
|
st.trim(starttime=stime, endtime=etime)
|
|
|
|
# choose 1 s taper
|
|
taper_fraction = 1 / (etime - stime)
|
|
for trace in st.traces:
|
|
# filter, remove response
|
|
trace.filter('bandpass', freqmin=0.01, freqmax=25, corners=3, zerophase=False). \
|
|
remove_response(output=output, zero_mean=True, taper=True, taper_fraction=taper_fraction)
|
|
return st
|
|
except FDSNException:
|
|
return None
|
|
|
|
|
|
def main(backend=None, args=None, wsgi=False):
|
|
"""
|
|
Main function to create a waveform plot.
|
|
This functions calls get_fdsn and fancy_plot to create the figure.
|
|
If wsgi == True it returns an ioString object containing the figure.
|
|
|
|
:type backend: str
|
|
:type args: dict
|
|
:type wsgi: bool
|
|
:rtype : object
|
|
:param backend:
|
|
:param args:
|
|
:param wsgi:
|
|
:return: ioString object with created image
|
|
"""
|
|
|
|
import warnings
|
|
import matplotlib as mpl
|
|
|
|
if wsgi:
|
|
backend = 'Agg'
|
|
if backend:
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore")
|
|
mpl.use(backend)
|
|
from obspy import UTCDateTime
|
|
|
|
if args:
|
|
deltat = args['length']
|
|
if args['stime']:
|
|
otime = UTCDateTime(args['stime'])
|
|
else:
|
|
otime = UTCDateTime() - 3600. - deltat
|
|
network = args['station'].split('.')[0]
|
|
station = args['station'].split('.')[1]
|
|
if args['type'] == 'dayplot':
|
|
if network == 'Z3':
|
|
channel = 'HHZ'
|
|
else:
|
|
channel = 'BHZ'
|
|
else:
|
|
if args['length'] < 3600.:
|
|
channel = 'HH?'
|
|
else:
|
|
if network == 'Z3':
|
|
channel = 'HH?'
|
|
else:
|
|
channel = 'BH?'
|
|
else:
|
|
otime = UTCDateTime() - 3600.
|
|
deltat = 30
|
|
network = 'GR'
|
|
station = 'BUG'
|
|
channel = 'HH?'
|
|
|
|
st = get_fdsn([(network, station, '*', channel, otime, otime + deltat)], base_url=args['server'])
|
|
if st is not None:
|
|
stime = max(st[0].stats.starttime, otime)
|
|
etime = min(st[0].stats.endtime, otime + deltat)
|
|
st.trim(starttime=stime, endtime=etime)
|
|
|
|
if wsgi:
|
|
if args['type'] == 'dayplot':
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore")
|
|
return trace_dayplot(st, outpattern=args['format'], wsgi=wsgi, deltat=deltat)
|
|
else:
|
|
return fancy_plot(st, wsgi=wsgi, img_format=args['format'], color=args['color'])
|
|
else:
|
|
if args['type'] == 'dayplot':
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore")
|
|
trace_dayplot(st, wsgi=wsgi)
|
|
else:
|
|
fancy_plot(st, color=args['color'])
|
|
elif not wsgi:
|
|
warnings.warn('No data available!')
|
|
|
|
|
|
def application(environ, start_response):
|
|
"""
|
|
Function application - Wrapper to process wsgi request
|
|
:param environ: contains information on the wsgi environment
|
|
:type environ: dict
|
|
:param start_response: function to process response header by the wsgi server
|
|
:type start_response: function
|
|
:return: response to be sent to the client by the wsgi server
|
|
:rtype: list
|
|
"""
|
|
|
|
import os
|
|
from cgi import FieldStorage
|
|
|
|
# set HOME environment variable to a directory the httpd server can write to
|
|
# otherwise matplotlib won't work
|
|
os.environ['HOME'] = '/tmp/wsgi-kasper'
|
|
try:
|
|
os.mkdir('/tmp/wsgi-kasper')
|
|
except OSError:
|
|
pass
|
|
|
|
# fake command line arguments
|
|
# setting needed default values
|
|
wsgi_args = {'server': 'http://localhost/'}
|
|
|
|
# fill wsgi_args with environment variables
|
|
form = FieldStorage(fp=environ['wsgi.input'], environ=environ)
|
|
wsgi_args['station'] = form.getfirst('station', 'GR.BUG').upper()
|
|
wsgi_args['stime'] = form.getfirst('start_time', [])
|
|
wsgi_args['length'] = form.getfirst('length', '30')
|
|
wsgi_args['length'] = abs(int(wsgi_args['length']))
|
|
wsgi_args['format'] = form.getfirst('format', 'png').lower()
|
|
if wsgi_args['format'] != 'png' and wsgi_args['format'] != 'svg':
|
|
wsgi_args['format'] = 'png'
|
|
wsgi_args['color'] = form.getfirst('color', 'TRUE').upper()
|
|
if wsgi_args['color'] == 'FALSE':
|
|
wsgi_args['color'] = False
|
|
else:
|
|
wsgi_args['color'] = True
|
|
wsgi_args['type'] = form.getfirst('type', '').lower()
|
|
if (wsgi_args['type'] != 'dayplot' and wsgi_args['type'] != 'normal') or wsgi_args['type'] == '':
|
|
if wsgi_args['length'] < 43200:
|
|
wsgi_args['type'] = 'normal'
|
|
else:
|
|
wsgi_args['type'] = 'dayplot'
|
|
|
|
# process the request
|
|
buf = main(args=wsgi_args, backend='Agg', wsgi=True)
|
|
if buf is not None:
|
|
data = buf.getvalue()
|
|
buf.close()
|
|
data_length = len(data)
|
|
|
|
if wsgi_args['format'] == 'svg':
|
|
wsgi_args['format'] = 'svg+xml'
|
|
start_response('200 OK', [('Content-Type', 'image/%s' % wsgi_args['format']), ('Content-Length', '%d' % data_length)])
|
|
return [data]
|
|
else:
|
|
start_response('400 Bad Request', [])
|
|
return []
|
|
|
|
# __main__
|
|
if __name__ == "__main__":
|
|
import os
|
|
import argparse
|
|
|
|
parser = argparse.ArgumentParser(
|
|
description=u'Get event waveform data of all RuhrNet stations.',
|
|
epilog=u'$Revision$ ($Date$, $Author$)'.replace(
|
|
"$", ""))
|
|
parser.add_argument(u'-v', u'-V', u'--version', action='version',
|
|
version=u'VVVVV')
|
|
parser.add_argument(u'-u', u'--url', action='store', dest='server',
|
|
default=u'https://ariadne.geophysik.ruhr-uni-bochum.de',
|
|
help=u'Base URL of the FDSN web service (https://ariadne.geophysik.ruhr-uni-bochum.de).')
|
|
parser.add_argument(u'-t', u'--type', action='store', dest='type', metavar='TYPE',
|
|
help=u'Type of plot: normal or dayplot.')
|
|
parser.add_argument(u'-f', u'--file', action='store', dest='filename', metavar='FILENAME',
|
|
help=u'Save plot to file FILENAME')
|
|
parser.add_argument(u'-c', u'--color', action='store_true', help=u'Create color plot.')
|
|
parser.add_argument(u'-l', u'--length', action='store', type=int, default=30,
|
|
help=u'Length of waveform window in seconds (30 s).')
|
|
group = parser.add_mutually_exclusive_group(required=True)
|
|
group.add_argument(u'-s', u'--start_time', dest='stime', action='store', metavar='START',
|
|
help=u'Start time of waveform window.')
|
|
group.add_argument(u'-e', u'--event', action='store', metavar='EVENTID',
|
|
help=u'Get starttime from event P-phase onset at station.')
|
|
parser.add_argument(u'station', action='store', metavar='NET.STATION',
|
|
help=u'Station to plot.')
|
|
|
|
cla = vars(parser.parse_args())
|
|
|
|
if os.getenv('DISPLAY'):
|
|
main(args=cla)
|
|
else:
|
|
main('Agg', cla)
|