#!/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 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)