diff --git a/datalogger.props.cfg b/datalogger.props.cfg new file mode 100644 index 0000000..22f53f0 --- /dev/null +++ b/datalogger.props.cfg @@ -0,0 +1,29 @@ +[connection] +port = /dev/ttyACM0 +baudrate = 9600 +timeout = 0.1 + +[data] +interval = 300 + +[device] +model = TC1 +latitude = 32.865468 +longitude = -117.253436 +altitude = 11 +network = IRI +station = SOCA +channel = BHZ +location = 00 +samplerate = 18.78 +dataquality = D + +[file] +datapath = /srv/tc1/data/ + +[logging] +logpath = /srv/tc1/log/ + +[calibration] +samplelimit = 405648 +offset = 34432 diff --git a/seisCalibration.py b/seisCalibration.py new file mode 100644 index 0000000..69a8bd0 --- /dev/null +++ b/seisCalibration.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import serial, sys, time, os.path, datetime +from time import strftime +import numpy as np +import ConfigParser as cp + +def printStats(alldata, alladj, counts, etime): + numsam = len(alldata) + (numpos, numneg, numzer) = counts + pospct = (numpos/(numsam-numzer))*100 + negpct = (numneg/(numsam-numzer))*100 + zerpct = (numzer/numsam)*100 + print("\n-----------------------------------------------------------------------------") + print "Running time seconds: %d" % etime + print "# of samples: %d" % numsam + print " Max/Adjusted max: %d / %d" % (np.amax(alldata),np.amax(alladj)) + print " Min/Adjusted min: %d / %d" % (np.amin(alldata),np.amin(alladj)) + print " Avg/Adjusted avg: %d / %d" % (np.mean(alldata),np.mean(alladj)) + print " # of positive: %d (%.2f%%)" % (numpos,pospct) + print " # of negative: %d (%.2f%%)" % (numneg,negpct) + print " # of zeroes: %d (%.2f%%)" % (numzer,zerpct) + print "System check completed." + +def msg(txt): + sys.stdout.write(txt) + sys.stdout.flush() + +def main(): + # Check if the data logger is running. If so, stop. + # We can't risk locking out the serial port from the data logger. + + if os.path.exists('/srv/tc1/run/seis_data_logger.pid'): + print "Error: Data Logger is already running! Exiting..." + exit(0) + + alldata = [] + alladj = [] + numpos = 0 + numneg = 0 + numzer = 0 + + config = cp.RawConfigParser() + config.read("/srv/tc1/conf/datalogger.props.cfg") + + port = config.get ("connection", "port") + baud = config.getint ("connection", "baudrate") + tout = config.getfloat("connection", "timeout") + offset = config.getint ("calibration","offset") + limit = config.getint ("calibration","samplelimit") + model = config.get ("device", "model") + srate = config.getfloat("device", "samplerate") + + ttlsec = limit / srate + m, s = divmod(ttlsec, 60) + h, m = divmod(m, 60) + now = datetime.datetime.now() + ftime = now + datetime.timedelta(hours=h, minutes=m, seconds=s) + (fhr,fmn,fsc) = (ftime.time().hour,ftime.time().minute,ftime.time().second) + + print("SIO Seismic Calibration Tool, v1.0") + print("* Initialized on " + strftime('%x at %X %Z')) + print("* Device info: model %s on %s (%d Hz transmission rate)" % (model, port, baud)) + print("* Using offset: %d" % offset) + print("* Sample limit: %d (approx. %d:%02d:%02d, finish at %d:%02d:%02d)" % (limit,h,m,s,fhr,fmn,fsc)) + print("-----------------------------------------------------------------------------") + print(" Sample Count | Sample Rate | Original Value | Adjusted Value | Current Mean ") + + device = serial.Serial(port,baud,timeout=tout) + stime = time.time() + try: + sampcount = 0 + runcount = 0 + curtotl = 0 + while len(alldata) < limit: + sproctime = time.time() + data = device.readline()[:-2] + sampcount = sampcount + 1 + if data: + eproctime = time.time() + + try: + rawdata = int(data) + runcount = runcount + 1 + alldata.append(rawdata) + adjdata = rawdata - offset + alladj.append(adjdata) + curtotl = curtotl + rawdata + curmean = int(curtotl / runcount) + if adjdata > 0: + numpos = numpos + 1 + elif adjdata < 0: + numneg = numneg + 1 + else: + numzer = numzer + 1 + + samprate = sampcount / (eproctime - sproctime) + + m = "{3:^13d} | {2:^11.2f} | {1:^14d} | {0:^14d} | {4:^14d}".format(adjdata,rawdata,samprate,runcount,curmean) + msg(m + chr(13)) + except ValueError: + continue + sampcount = 0 + except KeyboardInterrupt: + etime = time.time() + elapsed = etime - stime + counts = (float(numpos),float(numneg),float(numzer)) + printStats(alldata, alladj, counts, elapsed) + exit() + etime = time.time() + elapsed = etime - stime + counts = (float(numpos),float(numneg),float(numzer)) + printStats(alldata, alladj, counts, elapsed) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/seisDataLogger.py b/seisDataLogger.py new file mode 100644 index 0000000..f99c36c --- /dev/null +++ b/seisDataLogger.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# seisDataLogger.py +# +# Connects to educational seismometer and gathers data. +# Make save to current trace every user-defined interval of seconds +# Reads a config file for parameters, but runs as a process. +# + +from obspy.core import Stream, Trace, UTCDateTime +from seisDataLoggerDaemon import Daemon +from threading import Thread +from time import strftime +import serial, sys, os +import numpy as np + + +class DataLoggerDaemon(Daemon): + def run(self): + baseTime = UTCDateTime() + + # Assign all the configuration variables + port = self.config.get ("connection","port") + baud = self.config.getint ("connection","baudrate") + tout = self.config.getfloat("connection","timeout") + interval = self.config.getint ("data","interval") + offset = self.config.getint ("calibration","offset") + network = self.config.get ("device","network") + station = self.config.get ("device","station") + location = self.config.get ("device","location") + channel = self.config.get ("device","channel") + samprate = self.config.getfloat("device","samplerate") + dataqual = self.config.get ("device","dataquality") + + sampleIdx = 0 + traceData = np.array([]) + + self.logger.debug("["+ strftime('%X') + "] connecting...") + + rawData = serial.Serial(port,baud,timeout=tout) + + self.logger.debug("["+ strftime('%X') + "] listening for incoming data...") + + while True: # While loop that loops forever + while (rawData.inWaiting()==0): #Wait here until there is data + pass #do nothing + + dataPointString = rawData.readline() + + try: + traceData = np.append(traceData, int(dataPointString)) + except ValueError: + offset = int(np.mean(traceData)) + self.logger.debug("["+ strftime('%X') + "] * Bad value received. Replacing with current mean...") + traceData = np.append(traceData, offset) + + sampleIdx = sampleIdx + 1 + + currentTime = UTCDateTime() + elapsedTime = (currentTime - baseTime) + + # Write the data after x seconds + if elapsedTime >= (interval + (baseTime.microsecond / 1e6)): + # Fill header attributes + stats = {'network': network, + 'station': station, + 'location': location, + 'channel': channel, + 'npts': len(traceData), + 'sampling_rate': samprate, + 'mseed': {'dataquality': dataqual}, + 'starttime': baseTime} + + # Save the file using a different thread. + worker = Thread(target=self._writeData, args=(traceData, stats, baseTime)) + worker.setDaemon(True) + worker.start() + + baseTime = currentTime + + sampleIdx = 0 + + traceData = np.array([]) + + def _writeData(self, traceData, stats, timeObj): + streamObj = Stream([Trace(data=traceData, header=stats)]) + + filename = self._prepareFilename(timeObj) + offset = int(np.mean(streamObj.traces[0].data)) + streamObj.traces[0].data = np.array([x - offset for x in streamObj.traces[0].data]) + + self.logger.debug("["+ strftime('%X') + "] Saving %d samples (corrected by %d) to %s..." % (len(traceData), offset, filename)) + streamObj.write(filename, format='MSEED') + + def _prepareFilename(self, timeObj): + datapath = self.config.get("file","datapath") + filepath = datapath +"%d/%d/" % (timeObj.year, timeObj.julday) + + try: + if not os.path.exists(filepath): + os.makedirs(filepath) + except OSError as exception: + self.logger.debug("["+ strftime('%X') + "] * Error preparing path: (%d) %s" % (exception.errno, exception.strerror)) + + network = self.config.get("device","network") + station = self.config.get("device","station") + channel = self.config.get("device","channel") + filename = network+"."+station+".%02d%02d%4d_%02d%02d%02d." % \ + (timeObj.day,timeObj.month,timeObj.year, + timeObj.hour,timeObj.minute,timeObj.second) +channel+".mseed" + return (filepath+filename) + + def normalize(v): + norm=np.linalg.norm(v) + if norm==0: + return v + return v/norm + +if __name__ == "__main__": + daemon = DataLoggerDaemon('/srv/tc1/run/seis_data_logger.pid') + if len(sys.argv) == 2: + if 'start' == sys.argv[1]: + daemon.start() + elif 'stop' == sys.argv[1]: + daemon.stop() + elif 'restart' == sys.argv[1]: + daemon.restart() + else: + print "Unknown command" + sys.exit(2) + sys.exit(0) + else: + print "usage: %s start|stop|restart" % sys.argv[0] + sys.exit(2) \ No newline at end of file diff --git a/seisDataLoggerDaemon.py b/seisDataLoggerDaemon.py new file mode 100644 index 0000000..79271f6 --- /dev/null +++ b/seisDataLoggerDaemon.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import sys, os, atexit, logging, time +from time import strftime +from signal import SIGTERM +import ConfigParser as cp + +class Daemon: + """ + A generic daemon class. + + Usage: subclass the Daemon class and override the run() method + """ + def __init__(self, pidfile, stdin='/dev/null', \ + stdout='/srv/tc1/log/daemon.out.log', \ + stderr='/srv/tc1/log/daemon.err.log'): + self.stdin = stdin + self.stdout = stdout + self.stderr = stderr + self.pidfile = pidfile + + # Read and parse the config file. + self.config = cp.RawConfigParser() + self.config.read("/srv/tc1/conf/datalogger.props.cfg") + + # Set up logger + logpath = self.config.get("logging","logpath") + self.logger = logging.getLogger() + self.logger.setLevel(logging.DEBUG) + self.logger.addHandler(logging.FileHandler(logpath+"seisDataLogger.log")) + + def daemonize(self): + """ + do the UNIX double-fork magic, see Stevens' "Advanced + Programming in the UNIX Environment" for details (ISBN 0201563177) + http://www.erlenstar.demon.co.uk/unix/faq_2.html#SEC16 + """ + try: + pid = os.fork() + if pid > 0: + # exit first parent + sys.exit(0) + except OSError, e: + sys.stderr.write("fork #1 failed: %d (%s)\n" % (e.errno, e.strerror)) + sys.exit(1) + + # decouple from parent environment + os.chdir("/") + os.setsid() + os.umask(0) + + # do second fork + try: + pid = os.fork() + if pid > 0: + # exit from second parent + sys.exit(0) + except OSError, e: + sys.stderr.write("fork #2 failed: %d (%s)\n" % (e.errno, e.strerror)) + sys.exit(1) + + # redirect standard file descriptors + sys.stdout.flush() + sys.stderr.flush() + si = file(self.stdin, 'r') + so = file(self.stdout, 'a+') + se = file(self.stderr, 'a+', 0) + os.dup2(si.fileno(), sys.stdin.fileno()) + os.dup2(so.fileno(), sys.stdout.fileno()) + os.dup2(se.fileno(), sys.stderr.fileno()) + + # write pidfile + atexit.register(self.delpid) + pid = str(os.getpid()) + file(self.pidfile,'w+').write("%s\n" % pid) + + def delpid(self): + os.remove(self.pidfile) + + def start(self): + """ + Start the daemon + """ + # Check for a pidfile to see if the daemon already runs + try: + pf = file(self.pidfile,'r') + pid = int(pf.read().strip()) + pf.close() + except IOError: + pid = None + + if pid: + message = "pidfile %s already exist. Daemon already running?\n" + sys.stderr.write(message % self.pidfile) + sys.exit(1) + + # Start the daemon + self.daemonize() + + # Write log header + self.logger.info("-----------------------------------------") + self.logger.info("SIO Seismometer Data Logger Utility, v1.0") + self.logger.info("Initialized on " + strftime('%x at %X %Z')) + self.logger.info("Device type: " + self.config.get("device","model")) + self.logger.info("Lat: %f" % self.config.getfloat("device","latitude") + "°") + self.logger.info("Lon: %f" % self.config.getfloat("device","longitude") + "°") + self.logger.info("Alt: %.1f" % self.config.getfloat("device","altitude")+" meters") + + self.run() + + def stop(self): + """ + Stop the daemon + """ + self.logger.info("["+ strftime('%X') + "] shutdown requested...") + + # Get the pid from the pidfile + try: + pf = file(self.pidfile,'r') + pid = int(pf.read().strip()) + pf.close() + except IOError: + pid = None + + if not pid: + message = "pidfile %s does not exist. Daemon not running?\n" + sys.stderr.write(message % self.pidfile) + return # not an error in a restart + + # Try killing the daemon process + try: + while 1: + os.kill(pid, SIGTERM) + time.sleep(0.1) + except OSError, err: + err = str(err) + if err.find("No such process") > 0: + if os.path.exists(self.pidfile): + os.remove(self.pidfile) + else: + print str(err) + sys.exit(1) + + def restart(self): + """ + Restart the daemon + """ + self.stop() + self.start() + + def run(self): + """ + You should override this method when you subclass Daemon. It will be called after the process has been + daemonized by start() or restart(). + """ diff --git a/seisPlotAvg.py b/seisPlotAvg.py new file mode 100644 index 0000000..361d27f --- /dev/null +++ b/seisPlotAvg.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import serial, argparse +import numpy as np +from collections import deque +import ConfigParser as cp +import matplotlib.pyplot as plt +import matplotlib.animation as animation + + +# plot class +class AnalogPlot: + # constr + def __init__(self, strPort, maxLen): + # open serial port + self.ser = serial.Serial(strPort, 9600, timeout=.1) + + self.ax = deque([0.0]*maxLen) + self.ay = deque([0.0]*maxLen) + self.maxLen = maxLen + + config = cp.RawConfigParser() + config.read("/srv/tc1/conf/datalogger.props.cfg") + self.offset = config.getint("calibration","offset") + + self.dataArr = np.array([]) + self.sec = 0 + + # add to buffer + def addToBuf(self, buf, val): + if len(buf) < self.maxLen: + buf.append(val) + else: + buf.pop() + buf.appendleft(val) + + # add data + def add(self, data): + assert(len(data) == 2) + self.addToBuf(self.ax, data[0]) + self.addToBuf(self.ay, data[1]) + + # update plot + def update(self, frameNum, a0, a1): + try: + line = self.ser.readline() + try: + datapt = float(line) + #datapt = datapt - 33487 + if (datapt > (self.offset * 2) or datapt < (self.offset / 2)): + return self.offset + self.dataArr = np.append(self.dataArr, datapt) + self.sec = self.sec + 1 + data = [self.sec, self.dataArr.mean()] + # print data + self.add(data) + print "x: %d y: %d val: %d" % (data[0], data[1], datapt) + a0.set_data(range(self.maxLen), self.ax) + a1.set_data(range(self.maxLen), self.ay) + except ValueError: + pass + except KeyboardInterrupt: + print('exiting') + + return a0, + + # clean up + def close(self): + # close serial + self.ser.flush() + self.ser.close() + +# main() function +def main(): + config = cp.RawConfigParser() + config.read("/srv/tc1/conf/datalogger.props.cfg") + offset = config.getint("calibration","offset") + + strPort = '/dev/ttyACM0' + + print('reading from serial port %s...' % strPort) + + # plot parameters + analogPlot = AnalogPlot(strPort, 100) + + print('plotting data...') + # set up animation + fig = plt.figure() + ax = plt.axes(xlim=(0, 100), ylim=(offset-10, offset+10)) + a0, = ax.plot([], []) + a1, = ax.plot([], []) + anim = animation.FuncAnimation(fig, analogPlot.update, + fargs=(a0, a1), + interval=10) + + # show plot + plt.show() + + # clean up + analogPlot.close() + + print('exiting.') + + +# call main +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/seisPlotter.py b/seisPlotter.py new file mode 100644 index 0000000..56941f1 --- /dev/null +++ b/seisPlotter.py @@ -0,0 +1,198 @@ +#!/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 + "M%.1f %s (%.2f km depth), Distance from receiver: %.2f km
" % (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 = "No events currently displayed" + +imgurl = "http://%s/images/dayplots/iri.soca.%s.bhz.png" % (hostname, filedates) +htmlout = """ + + + + + + + + + + + + +
Earthquakes shown on this plot
%s
Note: Unlabeled seismic events are usually noise sources within the building.
Generated: %s on %s
+ + +""" % (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."