1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5# python 2/3
6from __future__ import absolute_import
8import os
9import math
10import logging
12import numpy as num
14from pyrocko import trace, util, plot
15from pyrocko.guts import Object, Int, String, Timestamp
17from . import io_common
19logger = logging.getLogger(__name__)
21N_GPS_TAGS_WANTED = 200 # must match definition in datacube_ext.c
24def color(c):
25 c = plot.color(c)
26 return tuple(x/255. for x in c)
29class DataCubeError(io_common.FileLoadError):
30 pass
33class ControlPointError(Exception):
34 pass
37def make_control_point(ipos_block, t_block, tref, deltat):
39 # reduce time (no drift would mean straight line)
40 tred = (t_block - tref) - ipos_block*deltat
42 # first round, remove outliers
43 q25, q75 = num.percentile(tred, (25., 75.))
44 iok = num.logical_and(q25 <= tred, tred <= q75)
46 # detrend
47 slope, offset = num.polyfit(ipos_block[iok], tred[iok], 1)
48 tred2 = tred - (offset + slope * ipos_block)
50 # second round, remove outliers based on detrended tred, refit
51 q25, q75 = num.percentile(tred2, (25., 75.))
52 iok = num.logical_and(q25 <= tred2, tred2 <= q75)
53 x = ipos_block[iok].copy()
54 ipos0 = x[0]
55 x -= ipos0
56 y = tred[iok].copy()
57 (slope, offset), cov = num.polyfit(x, y, 1, cov=True)
59 slope_err, offset_err = num.sqrt(num.diag(cov))
60 slope_err_limit = 1.0e-10
61 offset_err_limit = 5.0e-3
63 if slope_err > slope_err_limit:
64 raise ControlPointError('slope error too large')
66 if offset_err > offset_err_limit:
67 raise ControlPointError('offset error too large')
69 ic = ipos_block[ipos_block.size//2]
70 tc = offset + slope * (ic - ipos0)
72 return ic, tc + ic * deltat + tref
75def analyse_gps_tags(header, gps_tags, offset, nsamples):
77 ipos, t, fix, nsvs = gps_tags
78 deltat = 1.0 / int(header['S_RATE'])
80 tquartz = offset + ipos * deltat
82 toff = t - tquartz
83 toff_median = num.median(toff)
85 n = t.size
87 dtdt = (t[1:n] - t[0:n-1]) / (tquartz[1:n] - tquartz[0:n-1])
89 ok = abs(toff_median - toff) < 10.
91 xok = num.abs(dtdt - 1.0) < 0.00001
93 if ok.size >= 1:
95 ok[0] = False
96 ok[1:n] &= xok
97 ok[0:n-1] &= xok
98 ok[n-1] = False
100 ipos = ipos[ok]
101 t = t[ok]
102 fix = fix[ok]
103 nsvs = nsvs[ok]
105 blocksize = N_GPS_TAGS_WANTED // 2
107 try:
108 if ipos.size < blocksize:
109 raise ControlPointError(
110 'could not safely determine time corrections from gps')
112 j = 0
113 control_points = []
114 tref = num.median(t - ipos*deltat)
115 while j < ipos.size - blocksize:
116 ipos_block = ipos[j:j+blocksize]
117 t_block = t[j:j+blocksize]
118 try:
119 ic, tc = make_control_point(ipos_block, t_block, tref, deltat)
120 control_points.append((ic, tc))
121 except ControlPointError:
122 pass
123 j += blocksize
125 ipos_last = ipos[-blocksize:]
126 t_last = t[-blocksize:]
127 try:
128 ic, tc = make_control_point(ipos_last, t_last, tref, deltat)
129 control_points.append((ic, tc))
130 except ControlPointError:
131 pass
133 if len(control_points) < 2:
134 raise ControlPointError(
135 'could not safely determine time corrections from gps')
137 i0, t0 = control_points[0]
138 i1, t1 = control_points[1]
139 i2, t2 = control_points[-2]
140 i3, t3 = control_points[-1]
141 if len(control_points) == 2:
142 tmin = t0 - i0 * deltat - offset * deltat
143 tmax = t3 + (nsamples - i3 - 1) * deltat
144 else:
145 icontrol = num.array(
146 [x[0] for x in control_points], dtype=num.int64)
147 tcontrol = num.array(
148 [x[1] for x in control_points], dtype=float)
149 # robust against steps:
150 slope = num.median(
151 (tcontrol[1:] - tcontrol[:-1])
152 / (icontrol[1:] - icontrol[:-1]))
154 tmin = t0 + (offset - i0) * slope
155 tmax = t2 + (offset + nsamples - 1 - i2) * slope
157 if offset < i0:
158 control_points[0:0] = [(offset, tmin)]
160 if offset + nsamples - 1 > i3:
161 control_points.append((offset + nsamples - 1, tmax))
163 icontrol = num.array([x[0] for x in control_points], dtype=num.int64)
164 tcontrol = num.array([x[1] for x in control_points], dtype=float)
166 return tmin, tmax, icontrol, tcontrol, ok
168 except ControlPointError:
170 tmin = util.str_to_time(header['S_DATE'] + header['S_TIME'],
171 format='%y/%m/%d%H:%M:%S')
173 idat = int(header['DAT_NO'])
174 if idat == 0:
175 tmin = tmin + util.gps_utc_offset(tmin)
176 else:
177 tmin = util.day_start(tmin + idat * 24.*3600.) \
178 + util.gps_utc_offset(tmin)
180 tmax = tmin + (nsamples - 1) * deltat
181 icontrol, tcontrol = None, None
182 return tmin, tmax, icontrol, tcontrol, None
185def plot_gnss_location_timeline(fns):
186 from matplotlib import pyplot as plt
187 from pyrocko.orthodrome import latlon_to_ne_numpy
188 not_ = num.logical_not
189 h = 3600.
191 fig = plt.figure()
193 axes = []
194 for i in range(4):
195 axes.append(
196 fig.add_subplot(4, 1, i+1, sharex=axes[-1] if axes else None))
198 background_colors = [
199 color('aluminium1'),
200 color('aluminium2')]
202 tref = None
203 for ifn, fn in enumerate(fns):
204 header, gps_tags, nsamples = get_time_infos(fn)
205 _, t, fix, nsvs, lats, lons, elevations, _ = gps_tags
207 fix = fix.astype(bool)
209 if t.size < 2:
210 logger.warning('Need at least 2 gps tags for plotting: %s' % fn)
212 if tref is None:
213 tref = util.day_start(t[0])
214 lat, lon, elevation = coordinates_from_gps(gps_tags)
216 norths, easts = latlon_to_ne_numpy(lat, lon, lats, lons)
218 for ax, data in zip(axes, (norths, easts, elevations, nsvs)):
220 tspan = t[num.array([0, -1])]
222 ax.axvspan(*((tspan - tref) / h), color=background_colors[ifn % 2])
223 med = num.median(data)
224 ax.plot(
225 (tspan - tref) / h,
226 [med, med],
227 ls='--',
228 c='k',
229 lw=3,
230 alpha=0.5)
232 ax.plot(
233 (t[not_(fix)] - tref) / h, data[not_(fix)], 'o',
234 ms=1.5,
235 mew=0,
236 color=color('scarletred2'))
238 ax.plot(
239 (t[fix] - tref) / h, data[fix], 'o',
240 ms=1.5,
241 mew=0,
242 color=color('aluminium6'))
244 for ax in axes:
245 ax.grid(alpha=.3)
247 ax_lat, ax_lon, ax_elev, ax_nsv = axes
249 ax_lat.set_ylabel('Northing [m]')
250 ax_lon.set_ylabel('Easting [m]')
251 ax_elev.set_ylabel('Elevation [m]')
252 ax_nsv.set_ylabel('Number of Satellites')
254 ax_lat.get_xaxis().set_tick_params(labelbottom=False)
255 ax_lon.get_xaxis().set_tick_params(labelbottom=False)
256 ax_nsv.set_xlabel(
257 'Hours after %s' % util.time_to_str(tref, format='%Y-%m-%d'))
259 fig.suptitle(
260 u'Lat: %.5f\u00b0 Lon: %.5f\u00b0 Elevation: %g m' % (
261 lat, lon, elevation))
263 plt.show()
266def coordinates_from_gps(gps_tags):
267 ipos, t, fix, nsvs, lats, lons, elevations, temps = gps_tags
268 return tuple(num.median(x) for x in (lats, lons, elevations))
271def extract_stations(fns):
272 import io
273 import sys
274 from pyrocko.model import Station
275 from pyrocko.guts import dump_all
277 stations = {}
279 for fn in fns:
280 sta_name = os.path.splitext(fn)[1].lstrip('.')
281 if sta_name in stations:
282 logger.warning('Cube %s already in list!', sta_name)
283 continue
285 header, gps_tags, nsamples = get_time_infos(fn)
287 lat, lon, elevation = coordinates_from_gps(gps_tags)
289 sta = Station(
290 network='',
291 station=sta_name,
292 name=sta_name,
293 location='',
294 lat=lat,
295 lon=lon,
296 elevation=elevation)
298 stations[sta_name] = sta
300 f = io.BytesIO()
301 dump_all(stations.values(), stream=f)
302 sys.stdout.write(f.getvalue().decode())
305def plot_timeline(fns):
306 from matplotlib import pyplot as plt
308 fig = plt.figure()
309 axes = fig.gca()
311 h = 3600.
313 if isinstance(fns, str):
314 fn = fns
315 if os.path.isdir(fn):
316 fns = [
317 os.path.join(fn, entry) for entry in sorted(os.listdir(fn))]
319 ipos, t, fix, nsvs, header, offset, nsamples = \
320 get_timing_context(fns)
322 else:
323 ipos, t, fix, nsvs, header, offset, nsamples = \
324 get_extended_timing_context(fn)
326 else:
327 ipos, t, fix, nsvs, header, offset, nsamples = \
328 get_timing_context(fns)
330 deltat = 1.0 / int(header['S_RATE'])
332 tref = num.median(t - ipos * deltat)
333 tref = round(tref / deltat) * deltat
335 x = ipos*deltat
336 y = (t - tref) - ipos*deltat
338 bfix = fix != 0
339 bnofix = fix == 0
341 tmin, tmax, icontrol, tcontrol, ok = analyse_gps_tags(
342 header, (ipos, t, fix, nsvs), offset, nsamples)
344 la = num.logical_and
345 nok = num.logical_not(ok)
347 axes.plot(
348 x[la(bfix, ok)]/h, y[la(bfix, ok)], '+',
349 ms=5, color=color('chameleon3'))
350 axes.plot(
351 x[la(bfix, nok)]/h, y[la(bfix, nok)], '+',
352 ms=5, color=color('aluminium4'))
354 axes.plot(
355 x[la(bnofix, ok)]/h, y[la(bnofix, ok)], 'x',
356 ms=5, color=color('chocolate3'))
357 axes.plot(
358 x[la(bnofix, nok)]/h, y[la(bnofix, nok)], 'x',
359 ms=5, color=color('aluminium4'))
361 tred = tcontrol - icontrol*deltat - tref
362 axes.plot(icontrol*deltat/h, tred, color=color('aluminium6'))
363 axes.plot(icontrol*deltat/h, tred, 'o', ms=5, color=color('aluminium6'))
365 ymin = (math.floor(tred.min() / deltat)-1) * deltat
366 ymax = (math.ceil(tred.max() / deltat)+1) * deltat
368 # axes.set_ylim(ymin, ymax)
369 if ymax - ymin < 1000 * deltat:
370 ygrid = math.floor(tred.min() / deltat) * deltat
371 while ygrid < ymax:
372 axes.axhline(ygrid, color=color('aluminium4'), alpha=0.3)
373 ygrid += deltat
375 xmin = icontrol[0]*deltat/h
376 xmax = icontrol[-1]*deltat/h
377 xsize = xmax - xmin
378 xmin -= xsize * 0.1
379 xmax += xsize * 0.1
380 axes.set_xlim(xmin, xmax)
382 axes.set_ylim(ymin, ymax)
384 axes.set_xlabel('Uncorrected (quartz) time [h]')
385 axes.set_ylabel('Relative time correction [s]')
387 plt.show()
390g_dir_contexts = {}
393class DirContextEntry(Object):
394 path = String.T()
395 tstart = Timestamp.T()
396 ifile = Int.T()
399class DirContext(Object):
400 path = String.T()
401 mtime = Timestamp.T()
402 entries = DirContextEntry.T()
404 def get_entry(self, fn):
405 path = os.path.abspath(fn)
406 for entry in self.entries:
407 if entry.path == path:
408 return entry
410 raise Exception('entry not found')
412 def iter_entries(self, fn, step=1):
413 current = self.get_entry(fn)
414 by_ifile = dict(
415 (entry.ifile, entry) for entry in self.entries
416 if entry.tstart == current.tstart)
418 icurrent = current.ifile
419 while True:
420 icurrent += step
421 try:
422 yield by_ifile[icurrent]
424 except KeyError:
425 break
428def context(fn):
429 from pyrocko import datacube_ext
431 dpath = os.path.dirname(os.path.abspath(fn))
432 mtimes = [os.stat(dpath)[8]]
434 dentries = sorted([os.path.join(dpath, f) for f in os.listdir(dpath)
435 if os.path.isfile(os.path.join(dpath, f))])
436 for dentry in dentries:
437 fn2 = os.path.join(dpath, dentry)
438 mtimes.append(os.stat(fn2)[8])
440 mtime = float(max(mtimes))
442 if dpath in g_dir_contexts:
443 dir_context = g_dir_contexts[dpath]
444 if dir_context.mtime == mtime:
445 return dir_context
447 del g_dir_contexts[dpath]
449 entries = []
450 for dentry in dentries:
451 fn2 = os.path.join(dpath, dentry)
452 if not os.path.isfile(fn2):
453 continue
455 with open(fn2, 'rb') as f:
456 first512 = f.read(512)
457 if not detect(first512):
458 continue
460 with open(fn2, 'rb') as f:
461 try:
462 header, data_arrays, gps_tags, nsamples, _ = \
463 datacube_ext.load(f.fileno(), 3, 0, -1, None)
465 except datacube_ext.DataCubeError as e:
466 e = DataCubeError(str(e))
467 e.set_context('filename', fn)
468 raise e
470 header = dict(header)
471 entries.append(DirContextEntry(
472 path=os.path.abspath(fn2),
473 tstart=util.str_to_time(
474 '20' + header['S_DATE'] + ' ' + header['S_TIME'],
475 format='%Y/%m/%d %H:%M:%S'),
476 ifile=int(header['DAT_NO'])))
478 dir_context = DirContext(mtime=mtime, path=dpath, entries=entries)
480 return dir_context
483def get_time_infos(fn):
484 from pyrocko import datacube_ext
486 with open(fn, 'rb') as f:
487 try:
488 header, _, gps_tags, nsamples, _ = datacube_ext.load(
489 f.fileno(), 1, 0, -1, None)
491 except datacube_ext.DataCubeError as e:
492 e = DataCubeError(str(e))
493 e.set_context('filename', fn)
494 raise e
496 return dict(header), gps_tags, nsamples
499def get_timing_context(fns):
500 joined = [[], [], [], []]
501 ioff = 0
502 for fn in fns:
503 header, gps_tags, nsamples = get_time_infos(fn)
505 ipos = gps_tags[0]
506 ipos += ioff
508 for i in range(4):
509 joined[i].append(gps_tags[i])
511 ioff += nsamples
513 ipos, t, fix, nsvs = [num.concatenate(x) for x in joined]
515 nsamples = ioff
516 return ipos, t, fix, nsvs, header, 0, nsamples
519def get_extended_timing_context(fn):
520 c = context(fn)
522 header, gps_tags, nsamples_base = get_time_infos(fn)
524 ioff = 0
525 aggregated = [gps_tags]
527 nsamples_total = nsamples_base
529 if num.sum(gps_tags[2]) == 0:
531 ioff = nsamples_base
532 for entry in c.iter_entries(fn, 1):
534 _, gps_tags, nsamples = get_time_infos(entry.path)
536 ipos = gps_tags[0]
537 ipos += ioff
539 aggregated.append(gps_tags)
540 nsamples_total += nsamples
542 if num.sum(gps_tags[2]) > 0:
543 break
545 ioff += nsamples
547 ioff = 0
548 for entry in c.iter_entries(fn, -1):
550 _, gps_tags, nsamples = get_time_infos(entry.path)
552 ioff -= nsamples
554 ipos = gps_tags[0]
555 ipos += ioff
557 aggregated[0:0] = [gps_tags]
559 nsamples_total += nsamples
561 if num.sum(gps_tags[2]) > 0:
562 break
564 ipos, t, fix, nsvs = [num.concatenate(x) for x in zip(*aggregated)][:4]
566# return ipos, t, fix, nsvs, header, ioff, nsamples_total
567 return ipos, t, fix, nsvs, header, 0, nsamples_base
570def iload(fn, load_data=True, interpolation='sinc'):
571 from pyrocko import datacube_ext
572 from pyrocko import signal_ext
574 if interpolation not in ('sinc', 'off'):
575 raise NotImplementedError(
576 'no such interpolation method: %s' % interpolation)
578 with open(fn, 'rb') as f:
579 if load_data:
580 loadflag = 2
581 else:
582 loadflag = 1
584 try:
585 header, data_arrays, gps_tags, nsamples, _ = datacube_ext.load(
586 f.fileno(), loadflag, 0, -1, None)
588 except datacube_ext.DataCubeError as e:
589 e = DataCubeError(str(e))
590 e.set_context('filename', fn)
591 raise e
593 header = dict(header)
594 deltat = 1.0 / int(header['S_RATE'])
595 nchannels = int(header['CH_NUM'])
597 ipos, t, fix, nsvs, header_, offset_, nsamples_ = \
598 get_extended_timing_context(fn)
600 tmin, tmax, icontrol, tcontrol, _ = analyse_gps_tags(
601 header_, (ipos, t, fix, nsvs), offset_, nsamples_)
603 if icontrol is None:
604 logger.warn(
605 'No usable GPS timestamps found. Using datacube header '
606 'information to guess time. (file: "%s")' % fn)
608 tmin_ip = round(tmin / deltat) * deltat
609 if interpolation != 'off':
610 tmax_ip = round(tmax / deltat) * deltat
611 else:
612 tmax_ip = tmin_ip + (nsamples-1) * deltat
614 nsamples_ip = int(round((tmax_ip - tmin_ip)/deltat)) + 1
615 # to prevent problems with rounding errors:
616 tmax_ip = tmin_ip + (nsamples_ip-1) * deltat
618 leaps = num.array(
619 [x[0] + util.gps_utc_offset(x[0]) for x in util.read_leap_seconds2()],
620 dtype=float)
622 if load_data and icontrol is not None:
623 ncontrol_this = num.sum(
624 num.logical_and(0 <= icontrol, icontrol < nsamples))
626 if ncontrol_this <= 1:
627 logger.warn(
628 'Extrapolating GPS time information from directory context '
629 '(insufficient number of GPS timestamps in file: "%s").' % fn)
631 for i in range(nchannels):
632 if load_data:
633 arr = data_arrays[i]
634 assert arr.size == nsamples
636 if interpolation == 'sinc' and icontrol is not None:
638 ydata = num.empty(nsamples_ip, dtype=float)
639 try:
640 signal_ext.antidrift(
641 icontrol, tcontrol,
642 arr.astype(float), tmin_ip, deltat, ydata)
644 except signal_ext.Error as e:
645 e = DataCubeError(str(e))
646 e.set_context('filename', fn)
647 e.set_context('n_control_points', icontrol.size)
648 e.set_context('n_samples_raw', arr.size)
649 e.set_context('n_samples_ip', ydata.size)
650 e.set_context('tmin_ip', util.time_to_str(tmin_ip))
651 raise e
653 ydata = num.round(ydata).astype(arr.dtype)
654 else:
655 ydata = arr
657 tr_tmin = tmin_ip
658 tr_tmax = None
659 else:
660 ydata = None
661 tr_tmin = tmin_ip
662 tr_tmax = tmax_ip
664 tr = trace.Trace('', header['DEV_NO'], '', 'p%i' % i, deltat=deltat,
665 ydata=ydata, tmin=tr_tmin, tmax=tr_tmax, meta=header)
667 bleaps = num.logical_and(tmin_ip <= leaps, leaps < tmax_ip)
669 if num.any(bleaps):
670 assert num.sum(bleaps) == 1
671 tcut = leaps[bleaps][0]
673 for tmin_cut, tmax_cut in [
674 (tr.tmin, tcut), (tcut, tr.tmax+tr.deltat)]:
676 try:
677 tr_cut = tr.chop(tmin_cut, tmax_cut, inplace=False)
678 tr_cut.shift(
679 util.utc_gps_offset(0.5*(tr_cut.tmin+tr_cut.tmax)))
680 yield tr_cut
682 except trace.NoData:
683 pass
685 else:
686 tr.shift(util.utc_gps_offset(0.5*(tr.tmin+tr.tmax)))
687 yield tr
690header_keys = {
691 str: b'GIPP_V DEV_NO E_NAME GPS_PO S_TIME S_DATE DAT_NO'.split(),
692 int: b'''P_AMPL CH_NUM S_RATE D_FILT C_MODE A_CHOP F_TIME GPS_TI GPS_OF
693 A_FILT A_PHAS GPS_ON ACQ_ON V_TCXO D_VOLT E_VOLT'''.split()}
695all_header_keys = header_keys[str] + header_keys[int]
698def detect(first512):
699 s = first512
701 if len(s) < 512:
702 return False
704 if ord(s[0:1]) >> 4 != 15:
705 return False
707 n = s.find(b'\x80')
708 if n == -1:
709 n = len(s)
711 s = s[1:n]
712 s = s.replace(b'\xf0', b'')
713 s = s.replace(b';', b' ')
714 s = s.replace(b'=', b' ')
715 kvs = s.split(b' ')
717 if len([k for k in all_header_keys if k in kvs]) == 0:
718 return False
719 return True
722if __name__ == '__main__':
723 import argparse
724 parser = argparse.ArgumentParser(description='Datacube reader')
726 parser.add_argument(
727 'action', choices=['timeline', 'gnss', 'stations'],
728 help='Action')
729 parser.add_argument(
730 'files', nargs='+')
732 args = parser.parse_args()
733 if args.action == 'timeline':
734 plot_timeline(args.files)
736 elif args.action == 'gnss':
737 plot_gnss_location_timeline(args.files)
739 elif args.action == 'stations':
740 extract_stations(args.files)