From 23ab30907f335edae5643cdecdebb21aa313d0f8 Mon Sep 17 00:00:00 2001 From: "Kasper D. Fischer" Date: Tue, 10 Mar 2015 13:25:41 +0000 Subject: [PATCH] Branching utils/trunk/python scripts to www/trunk/wsgi to start plotting webapp. --- wsgi/plotFDSN.py | 322 +++++++++++++++++++++++++++++++++++++++++++ wsgi/traceDayplot.py | 114 +++++++++++++++ 2 files changed, 436 insertions(+) create mode 100755 wsgi/plotFDSN.py create mode 100755 wsgi/traceDayplot.py diff --git a/wsgi/plotFDSN.py b/wsgi/plotFDSN.py new file mode 100755 index 0000000..8ffd530 --- /dev/null +++ b/wsgi/plotFDSN.py @@ -0,0 +1,322 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" + Get waveformdata from FDSN web service and create a fancy plot + This programm runs as a script or as a WSGI application. + + Subversion information: + $Id$ + + :license + Copyright 2015 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 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, cla=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 cla: dict + :type wsgi: bool + :rtype : object + :param backend: + :param cla: + :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 cla: + deltat = cla['length'] + if cla['stime']: + otime = UTCDateTime(cla['stime']) + else: + otime = UTCDateTime() - 3600 - deltat + network = cla['station'].split('.')[0] + station = cla['station'].split('.')[1] + else: + otime = UTCDateTime('2014-11-15T11:35:25Z') + deltat = 30 + network = 'GR' + station = 'BUG' + + st = get_fdsn([(network, station, '*', 'HH?', otime, otime + deltat)], base_url=cla['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: + # noinspection PyUnresolvedReferences + return fancy_plot(st, wsgi=wsgi, img_format=cla['format'], color=cla['color']) + else: + # noinspection PyUnresolvedReferences + fancy_plot(st, color=cla['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 + cla = {'server': 'http://localhost/'} + + # fill cla with environment variables + form = FieldStorage(fp=environ['wsgi.input'], environ=environ) + cla['station'] = form.getfirst('station', 'GR.BUG').upper() + cla['stime'] = form.getfirst('start_time', []) + cla['length'] = form.getfirst('length', '30') + cla['length'] = abs(int(cla['length'])) + cla['format'] = form.getfirst('format', 'png').lower() + if cla['format'] != 'png' and cla['format'] != 'svg': + cla['format'] = 'png' + cla['color'] = form.getfirst('color', 'TRUE').upper() + if cla['color'] == 'FALSE': + cla['color'] = False + else: + cla['color'] = True + + # process the request + buf = main(cla=cla, backend='Agg', wsgi=True) + if buf is not None: + data = buf.getvalue() + # noinspection PyUnresolvedReferences + buf.close() + data_length = len(data) + + if cla['format'] == u'svg': + cla['format'] = u'svg+xml' + start_response('200 OK', [('Content-Type', 'image/%s' % cla['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'$Revision$ ($Date$, \ + $Author$)'.replace('$', '')) + 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'-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(cla=cla) + else: + main('Agg', cla) \ No newline at end of file diff --git a/wsgi/traceDayplot.py b/wsgi/traceDayplot.py new file mode 100755 index 0000000..dfb51ac --- /dev/null +++ b/wsgi/traceDayplot.py @@ -0,0 +1,114 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- +""" +Produce a dayplot + +Subversion information: +$Id$ +""" +""" +Copyright 2012 Seismological Observatory, Ruhr-University Bochum +http://www.gmg.ruhr-uni-bochum.de/geophysik/seisobs +Contributors: + Martina Rische + Kasper D. Fischer + Sebastian Wehling-Benatelli + +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 3 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, see http://www.gnu.org/licenses/. +""" + +def traceDayplot(datafile, ftype='bandpass', fmin=1.0, fmax=7.0, + col=('b','r','g'), interval=20.0, outpattern=''): + + ''' + ''' + + + from obspy.core import read + + st = read(datafile) + + #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() + + #plot + for i in range(0,len(st)): + 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) + + +# __main__ +if __name__ == "__main__": + + # parse arguments + import argparse + parser = argparse.ArgumentParser( + description='Produce filtered 24h-plot (dayplot).', + epilog='$Rev: 403 $ ($Date: 2012-04-13 12:16:22 +0200 (Fri, 13 Apr 2012) $, $Author: kasper $)') + parser.add_argument('-v', '-V','--version', action='version', + version='$Rev: 403 $ ($Date: 2012-04-13 12:16:22 +0200 (Fri, 13 Apr 2012) $, $Author: kasper $)') + parser.add_argument('file', action='store', metavar='FILE', nargs='+', + help='File(s) to use for the dayplot. One dayplot will be used for every file. \ + Use wildcards to use multiple file for one plot') + parser.add_argument('-f', '--filter', action='store', dest='ftype', + default='bandpass', + choices=['none', 'bandpass', 'bandstop', 'lowpass', 'highpass'], + help='Select filtertype to filter the data (default: bandpass)\ + Note: For low- and highpass only fmin is used.') + parser.add_argument('--fmin', action='store', type=float, dest='fmin', + default=1.0, + help='Lower frequency of the filter in Hz (default: 1.0)') + parser.add_argument('--fmax', action='store', type=float, dest='fmax', + default=7.0, + help='Upper frequency of the filter in Hz (default: 7.0)') + parser.add_argument('-c', '--color', '--colour', action='store', dest='color', + default=('b', 'r', 'g'), + help='Color selection to use in the dayplot (default: brg)') + parser.add_argument('-i', '--interval', action='store', type=int, dest='interval', + default=20, + help='Interval length to show in each line of the dayplot in minutes (default: 20)') + parser.add_argument('-o', '--output', action='store', dest='outpattern', + default='', + help="Output filename pattern for the plot. (default: unset, show plot on screen only), \ + Use {0} to substitute network code, \ + {1} to substitute station code, \ + {2} to subsitute station channel, \ + {3} to substitute year, \ + {4} to substitute doy number \ + .format (supported formats: emf, eps, pdf, png, ps, raw, rgba, svg, svgz)") + + cla = parser.parse_args() + if cla.fmin < 0: + cla.fmin = -cla.fmin + if cla.fmax < 0: + cla.fmax = -cla.fmax + + # call traceDayplot(...) for each file + for datafile in cla.file: + traceDayplot(datafile, cla.ftype, cla.fmin, cla.fmax, cla.color, + cla.interval, cla.outpattern) +