seisPlotter.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# seisPlotter.py
# 
# Reads a bunch of seed files based on the time, and then displays a 12 hour
# helicorder plot ("dayplot").
#
# TODO:
# * Parse through the files on the server to make the dayplot
#       Put file in /srv/www/out
# * Generate a small XHTML page that can be imported containing:
#       Dayplot
#       List of Earthquakes in plot
#           "MXX LOCATION (XX km depth), Dist from receiver: XX km"
#           link location to USGS page (in feature['properties']['url']
#           or "No Earthquakes currently displayed"
#
import sys
from time import strftime, time, localtime

start_execute = time()

sys.stdout.write('loading dependencies... ')
sys.stdout.flush()
import matplotlib
matplotlib.use('Agg')
import datetime, urllib, json, subprocess
from obspy import read, UTCDateTime#, #Stream, Trace
from geopy.distance import vincenty

sys.stdout.write('done.\n')
sys.stdout.flush()

seismometer = (32.865468, -117.253436)
eventcoll = []
outstring = ""

# Get the hostname
p1=subprocess.check_output(['uname','-n']).rstrip()
p2=subprocess.check_output(['awk','/^domain/ {print $2}','/etc/resolv.conf']).rstrip()
hostname = p1 + "." + p2

datapath = "/srv/tc1/data/"

sys.stdout.write('querying USGS earthquake feed... ')
sys.stdout.flush()
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.geojson"
response = urllib.urlopen(url);
data = json.loads(response.read())
sys.stdout.write('done.\n')
sys.stdout.flush()

sys.stdout.write('filtering results by distance/magnitude... ')
sys.stdout.flush()

now   = datetime.datetime.now() + datetime.timedelta(hours=7, minutes=0, seconds=0) # convert time to UTC
past  = now - datetime.timedelta(hours=12, minutes=now.time().minute, seconds=now.time().second, microseconds=now.time().microsecond)
start = UTCDateTime(past) # put in UTCDateTime object
end   = UTCDateTime(now)

for feature in data['features']:
    place           = feature['properties']['place']
    mag             = feature['properties']['mag']
    url             = feature['properties']['url']
    (lon,lat,depth) = feature['geometry']['coordinates']
    eventtime       = UTCDateTime(feature['properties']['time']/1000)

    if (end - eventtime) <= 46800:
        # Decide which earthquakes to show
        show_quake = False
        dist = vincenty(seismometer, (lat,lon)).kilometers
        if (dist >= 8000):
            if mag >= 6.0:
                show_quake = True
        elif (dist < 8000 and dist >= 5000):
            if mag >= 5.9:
                show_quake = True
        elif (dist < 5000 and dist >= 3000):
            if mag >= 5.8: 
                show_quake = True
        elif (dist < 3000 and dist >= 1500):
            if mag >= 5.6:
                show_quake = True
        elif (dist < 1500 and dist >= 1000):
            if mag >= 4.7:
                show_quake = True
        elif (dist < 1000 and dist >= 300):
            if mag >= 3.5:
                show_quake = True
        elif (dist < 300 and dist >= 150):
            if mag >= 3.2:
                show_quake = True
        elif (dist < 150):
            if mag >= 3:
                show_quake = True
        
        if show_quake:
            # add to event collection
            event = {"time": eventtime, "text": ("M%s" % mag) + " " + place}
            eventcoll.append(event)
            outstring = outstring + "<a href=\"%s\" target=\"_new\">M%.1f %s (%.2f km depth), Distance from receiver: %.2f km</a><br />" % (url,mag,place,depth,dist)
sys.stdout.write('done.\n')
sys.stdout.flush()

yr    = end.year
ejday = end.julday
sjday = start.julday

#print "DEBUG: ejday: %d, sjday: %d" % (ejday,sjday)

sys.stdout.write('reading files... ')
sys.stdout.flush()

st = read(datapath + "%d/%d/*.mseed" % (yr,ejday), format="MSEED", starttime=start, endtime=end)
if sjday != ejday:
    st1 = read(datapath + "%d/%d/*.mseed" % (yr,sjday), format="MSEED", starttime=start, endtime=end)
    st = st + st1
    
sys.stdout.write('done.\n')
sys.stdout.flush()

numtraces = 0
try:
    numtraces = len(st.traces)
except AttributeError:
    sys.stdout.write('Error: There are no traces to gather.')
    sys.stdout.write('\nExiting process...\n')
    sys.stdout.flush()
    sys.exit(1)
    
ftime = st.traces[0].stats.starttime
ltime = st.traces[numtraces-1].stats.endtime
(fyear,fmon,fday) = (ftime.year,ftime.month,ftime.day)
(lyear,lmon,lday) = (ltime.year,ltime.month,ltime.day)

filedates = "%d%02d%02d-%d%02d%02d" % (lyear,lmon,lday,fyear,fmon,fday)
outfile = "/srv/www/%s/public/images/dayplots/iri.soca.%s.bhz.png" % (hostname, filedates)

sys.stdout.write('filtering data... ')
sys.stdout.flush()
st.filter("lowpass", freq=0.1, corners=2)
sys.stdout.write('done.\n')
sys.stdout.flush()

sys.stdout.write('generating plot... ')
sys.stdout.flush()
st.plot(type="dayplot", interval=60, right_vertical_labels=False,
        #title="%s (%s)" % (st.traces[0].id,hostname),
        title="",
        vertical_scaling_range=15, one_tick_per_line=True,
        color=['k', 'r', 'b', 'g'], 
        show_y_UTC_label=True,
        events=eventcoll,
        outfile=outfile
        )
sys.stdout.write('done.\n')
sys.stdout.flush()

timestring = strftime("%a, %d %b %Y %H:%M:%S %Z", localtime())

if len(eventcoll) == 0:
    outstring = "<i>No events currently displayed</i>"

imgurl = "http://%s/images/dayplots/iri.soca.%s.bhz.png" % (hostname, filedates)
htmlout = """
<html>
<head>
<meta http-equiv="Refresh" content="300">
<style>
table, th, td {
    border: 0px solid black;
    border-collapse: collapse;
}
</style>
</head>
<body>
<table>
<tr><td><img src="%s" /></td></tr>
<tr><th>Earthquakes shown on this plot</th></tr>
<tr><td style="text-align:center">%s</td></tr>
<tr><td style="text-align:center"><i>Note: Unlabeled seismic events are usually noise sources within the building.</i></td></tr>
<tr><td style="text-align:center"><i>Generated: %s on %s</i></td></tr>
</table>
</body>
</html>
""" % (imgurl, outstring, timestring, hostname)
sys.stdout.write('writing include file... ')
sys.stdout.flush()
f = open('/srv/www/%s/public/out/tc1.html' % hostname, 'w')
f.write(htmlout)
f.close()
sys.stdout.write('done.\n')
sys.stdout.write('processing finished.\n')

end_execute = time()
elapsed = end_execute - start_execute
print "Elapsed run time: ", elapsed, "seconds."