Branching utils/trunk/python scripts to www/trunk/wsgi to start plotting webapp.

r819 www/trunk

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 <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 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)

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 <martina.rische@rub.de>
+    Kasper D. Fischer <kasper.fischer@rub.de>
+    Sebastian Wehling-Benatelli <sebastian.wehling@rub.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 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)
+