Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/io/datacube.py: 57%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Reader for `DiGOS DATA-CUBE³ <https://digos.eu/the-seismic-data-recorder/>`_ 

8raw data. 

9''' 

10 

11 

12import os 

13import math 

14import logging 

15import warnings 

16 

17import numpy as num 

18 

19from pyrocko import trace, util, plot 

20from pyrocko.guts import Object, Int, String, Timestamp 

21from typing import NamedTuple 

22 

23from . import io_common 

24 

25logger = logging.getLogger(__name__) 

26 

27HOUR = 3600. 

28DAY = 24 * HOUR 

29WEEK = 7 * DAY 

30WNRO = 1024 * WEEK # GPS Week Number Rollover period 

31 

32N_GPS_TAGS_WANTED = 200 # must match definition in datacube_ext.c 

33 

34APPLY_SUBSAMPLE_SHIFT_CORRECTION = True 

35 

36 

37def color(c): 

38 c = plot.color(c) 

39 return tuple(x/255. for x in c) 

40 

41 

42class DataCubeError(io_common.FileLoadError): 

43 pass 

44 

45 

46class ControlPointError(Exception): 

47 pass 

48 

49 

50class GPSTags(NamedTuple): 

51 position: num.ndarray 

52 gps_time: num.ndarray 

53 fix: num.ndarray 

54 sat_count: num.ndarray 

55 latitudes: num.ndarray 

56 longitudes: num.ndarray 

57 elevations: num.ndarray 

58 temperatures: num.ndarray 

59 gps_utc_offsets: num.ndarray 

60 

61 @classmethod 

62 def from_tuple(cls, tup): 

63 return cls(*tup) 

64 

65 

66def check_WNRO(utc_gps_offsets: num.ndarray, leap_seconds: int) -> bool: 

67 """Checks for GPS week number rollover offset.""" 

68 # give some tolerance, we expect larger offsets from 1024 week rollovers 

69 max_utc_gps_offset = utc_gps_offsets.max() 

70 if max_utc_gps_offset == 0: 

71 return False 

72 return abs(leap_seconds + max_utc_gps_offset) > 1 

73 

74 

75def make_control_point(ipos_block, t_block, tref, deltat): 

76 

77 # reduce time (no drift would mean straight line) 

78 tred = (t_block - tref) - ipos_block*deltat 

79 

80 # first round, remove outliers 

81 q25, q75 = num.percentile(tred, (25., 75.)) 

82 iok = num.logical_and(q25 <= tred, tred <= q75) 

83 

84 # detrend 

85 slope, offset = num.polyfit(ipos_block[iok], tred[iok], 1) 

86 tred2 = tred - (offset + slope * ipos_block) 

87 

88 # second round, remove outliers based on detrended tred, refit 

89 q25, q75 = num.percentile(tred2, (25., 75.)) 

90 iok = num.logical_and(q25 <= tred2, tred2 <= q75) 

91 x = ipos_block[iok].copy() 

92 ipos0 = x[0] 

93 x -= ipos0 

94 y = tred[iok].copy() 

95 if x.size < 2: 

96 raise ControlPointError('Insufficient number control points after QC.') 

97 

98 elif x.size < 5: # needed for older numpy versions 

99 (slope, offset) = num.polyfit(x, y, 1) 

100 else: 

101 (slope, offset), cov = num.polyfit(x, y, 1, cov=True) 

102 

103 slope_err, offset_err = num.sqrt(num.diag(cov)) 

104 

105 slope_err_limit = 1.0e-10 

106 if ipos_block.size < N_GPS_TAGS_WANTED // 2: 

107 slope_err_limit *= (200. / ipos_block.size) 

108 

109 offset_err_limit = 5.0e-3 

110 

111 if slope_err > slope_err_limit: 

112 raise ControlPointError( 

113 'Slope error too large: %g (limit: %g' % ( 

114 slope_err, slope_err_limit)) 

115 

116 if offset_err > offset_err_limit: 

117 raise ControlPointError( 

118 'Offset error too large: %g (limit: %g' % ( 

119 offset_err, offset_err_limit)) 

120 

121 ic = ipos_block[ipos_block.size//2] 

122 tc = offset + slope * (ic - ipos0) 

123 

124 return ic, tc + ic * deltat + tref 

125 

126 

127def analyse_gps_tags(header, gps_tags, offset, nsamples): 

128 

129 ipos, times, fix, sat_count = gps_tags 

130 deltat = 1.0 / int(header['S_RATE']) 

131 

132 tquartz = offset + ipos * deltat 

133 

134 toff = times - tquartz 

135 toff_median = num.median(toff) 

136 

137 n = times.size 

138 

139 dtdt = (times[1:n] - times[0:n-1]) / (tquartz[1:n] - tquartz[0:n-1]) 

140 

141 ok = abs(toff_median - toff) < 10. 

142 

143 xok = num.abs(dtdt - 1.0) < 0.00001 

144 

145 if ok.size >= 1: 

146 

147 ok[0] = False 

148 ok[1:n] &= xok 

149 ok[0:n-1] &= xok 

150 ok[n-1] = False 

151 

152 ipos = ipos[ok] 

153 times = times[ok] 

154 fix = fix[ok] 

155 sat_count = sat_count[ok] 

156 

157 blocksize = N_GPS_TAGS_WANTED // 2 

158 if ipos.size < blocksize: 

159 blocksize = max(10, ipos.size // 4) 

160 warnings.warn( 

161 'Small number of GPS tags found. ' 

162 'Reducing analysis block size to %i tags. ' 

163 'Time correction may be unreliable.' % blocksize) 

164 

165 try: 

166 if ipos.size < blocksize: 

167 raise ControlPointError( 

168 'Cannot determine GPS time correction: ' 

169 'Too few GPS tags found: %i' % ipos.size) 

170 

171 j = 0 

172 control_points = [] 

173 tref = num.median(times - ipos*deltat) 

174 while j < ipos.size - blocksize: 

175 ipos_block = ipos[j:j+blocksize] 

176 t_block = times[j:j+blocksize] 

177 try: 

178 ic, tc = make_control_point(ipos_block, t_block, tref, deltat) 

179 control_points.append((ic, tc)) 

180 except ControlPointError as e: 

181 logger.debug(str(e)) 

182 

183 j += blocksize 

184 

185 ipos_last = ipos[-blocksize:] 

186 t_last = times[-blocksize:] 

187 try: 

188 ic, tc = make_control_point(ipos_last, t_last, tref, deltat) 

189 control_points.append((ic, tc)) 

190 except ControlPointError as e: 

191 logger.debug(str(e)) 

192 

193 if len(control_points) < 2: 

194 raise ControlPointError( 

195 'Could not safely determine time corrections from GPS: ' 

196 'unable to construct two or more control points') 

197 

198 i0, t0 = control_points[0] 

199 i1, t1 = control_points[1] 

200 i2, t2 = control_points[-2] 

201 i3, t3 = control_points[-1] 

202 if len(control_points) == 2: 

203 tmin = t0 - i0 * deltat - offset * deltat 

204 tmax = t3 + (nsamples - i3 - 1) * deltat 

205 else: 

206 icontrol = num.array( 

207 [x[0] for x in control_points], dtype=num.int64) 

208 tcontrol = num.array( 

209 [x[1] for x in control_points], dtype=float) 

210 # robust against steps: 

211 slope = num.median( 

212 (tcontrol[1:] - tcontrol[:-1]) 

213 / (icontrol[1:] - icontrol[:-1])) 

214 

215 tmin = t0 + (offset - i0) * slope 

216 tmax = t2 + (offset + nsamples - 1 - i2) * slope 

217 

218 if offset < i0: 

219 control_points[0:0] = [(offset, tmin)] 

220 

221 if offset + nsamples - 1 > i3: 

222 control_points.append((offset + nsamples - 1, tmax)) 

223 

224 icontrol = num.array([x[0] for x in control_points], dtype=num.int64) 

225 tcontrol = num.array([x[1] for x in control_points], dtype=float) 

226 

227 # corrected 2021-10-26: This sub-sample time shift introduced by the 

228 # Cube's ADC was previously not recognized. 

229 if APPLY_SUBSAMPLE_SHIFT_CORRECTION: 

230 tcontrol -= 0.199 * deltat + 0.0003 

231 

232 return tmin, tmax, icontrol, tcontrol, ok 

233 

234 except ControlPointError as e: 

235 logger.error(str(e)) 

236 

237 tmin = util.str_to_time(header['S_DATE'] + header['S_TIME'], 

238 format='%y/%m/%d%H:%M:%S') 

239 

240 idat = int(header['DAT_NO']) 

241 if idat == 0: 

242 tmin = tmin + util.gps_utc_offset(tmin) 

243 else: 

244 tmin = util.day_start(tmin + idat * DAY) \ 

245 + util.gps_utc_offset(tmin) 

246 

247 tmax = tmin + (nsamples - 1) * deltat 

248 icontrol = num.array([offset, offset+nsamples - 1], dtype=num.int64) 

249 tcontrol = num.array([tmin, tmax]) 

250 return tmin, tmax, icontrol, tcontrol, ok 

251 

252 

253def plot_gnss_location_timeline(fns): 

254 from matplotlib import pyplot as plt 

255 from pyrocko.orthodrome import latlon_to_ne_numpy 

256 not_ = num.logical_not 

257 

258 fig = plt.figure() 

259 

260 axes = [] 

261 for i in range(4): 

262 axes.append( 

263 fig.add_subplot(4, 1, i+1, sharex=axes[-1] if axes else None)) 

264 

265 background_colors = [ 

266 color('aluminium1'), 

267 color('aluminium2')] 

268 

269 tref = None 

270 for ifn, fn in enumerate(fns): 

271 header, gps_tags, nsamples = get_time_infos(fn) 

272 _, t, fix, nsvs, lats, lons, elevations, _, gps_utc_offset = gps_tags 

273 leap_seconds = util.utc_gps_offset(num.median(t)) 

274 if check_WNRO(gps_utc_offset, leap_seconds): 

275 t += WNRO 

276 

277 fix = fix.astype(bool) 

278 

279 if t.size < 2: 

280 logger.warning('Need at least 2 gps tags for plotting: %s' % fn) 

281 

282 if tref is None: 

283 tref = util.day_start(t[0]) 

284 lat, lon, elevation = coordinates_from_gps(gps_tags) 

285 

286 norths, easts = latlon_to_ne_numpy(lat, lon, lats, lons) 

287 

288 for ax, data in zip(axes, (norths, easts, elevations, nsvs)): 

289 

290 tspan = t[num.array([0, -1])] 

291 

292 ax.axvspan(*((tspan - tref) / HOUR), 

293 color=background_colors[ifn % 2]) 

294 med = num.median(data) 

295 ax.plot( 

296 (tspan - tref) / HOUR, 

297 [med, med], 

298 ls='--', 

299 c='k', 

300 lw=3, 

301 alpha=0.5) 

302 

303 ax.plot( 

304 (t[not_(fix)] - tref) / HOUR, data[not_(fix)], 'o', 

305 ms=1.5, 

306 mew=0, 

307 color=color('scarletred2')) 

308 

309 ax.plot( 

310 (t[fix] - tref) / HOUR, data[fix], 'o', 

311 ms=1.5, 

312 mew=0, 

313 color=color('aluminium6')) 

314 

315 for ax in axes: 

316 ax.grid(alpha=.3) 

317 

318 ax_lat, ax_lon, ax_elev, ax_nsv = axes 

319 

320 ax_lat.set_ylabel('Northing [m]') 

321 ax_lon.set_ylabel('Easting [m]') 

322 ax_elev.set_ylabel('Elevation [m]') 

323 ax_nsv.set_ylabel('Number of Satellites') 

324 

325 ax_lat.get_xaxis().set_tick_params(labelbottom=False) 

326 ax_lon.get_xaxis().set_tick_params(labelbottom=False) 

327 ax_nsv.set_xlabel( 

328 'Hours after %s' % util.time_to_str(tref, format='%Y-%m-%d')) 

329 

330 fig.suptitle( 

331 u'Lat: %.5f\u00b0 Lon: %.5f\u00b0 Elevation: %g m' % ( 

332 lat, lon, elevation)) 

333 

334 plt.show() 

335 

336 

337def coordinates_from_gps(gps_tags): 

338 ipos, t, fix, nsvs, lats, lons, elevations, temps, gps_utc_offset = \ 

339 gps_tags 

340 return tuple(num.median(x) for x in (lats, lons, elevations)) 

341 

342 

343def extract_stations(fns): 

344 import io 

345 import sys 

346 from pyrocko.model import Station 

347 from pyrocko.guts import dump_all 

348 

349 stations = {} 

350 

351 for fn in fns: 

352 sta_name = os.path.splitext(fn)[1].lstrip('.') 

353 if sta_name in stations: 

354 logger.warning('Cube %s already in list!', sta_name) 

355 continue 

356 

357 header, gps_tags, nsamples = get_time_infos(fn) 

358 

359 lat, lon, elevation = coordinates_from_gps(gps_tags) 

360 

361 sta = Station( 

362 network='', 

363 station=sta_name, 

364 name=sta_name, 

365 location='', 

366 lat=lat, 

367 lon=lon, 

368 elevation=elevation) 

369 

370 stations[sta_name] = sta 

371 

372 f = io.BytesIO() 

373 dump_all(stations.values(), stream=f) 

374 sys.stdout.write(f.getvalue().decode()) 

375 

376 

377def plot_timeline(fns): 

378 from matplotlib import pyplot as plt 

379 

380 fig = plt.figure() 

381 axes = fig.gca() 

382 

383 if isinstance(fns, str): 

384 fn = fns 

385 if os.path.isdir(fn): 

386 fns = [ 

387 os.path.join(fn, entry) for entry in sorted(os.listdir(fn))] 

388 

389 ipos, t, fix, nsvs, header, offset, nsamples = \ 

390 get_timing_context(fns) 

391 

392 else: 

393 ipos, t, fix, nsvs, header, offset, nsamples = \ 

394 get_extended_timing_context(fn) 

395 

396 else: 

397 ipos, t, fix, nsvs, header, offset, nsamples = \ 

398 get_timing_context(fns) 

399 

400 deltat = 1.0 / int(header['S_RATE']) 

401 

402 tref = num.median(t - ipos * deltat) 

403 tref = round(tref / deltat) * deltat 

404 

405 if APPLY_SUBSAMPLE_SHIFT_CORRECTION: 

406 tcorr = 0.199 * deltat + 0.0003 

407 else: 

408 tcorr = 0.0 

409 

410 x = ipos*deltat 

411 y = (t - tref) - ipos*deltat - tcorr 

412 

413 bfix = fix != 0 

414 bnofix = fix == 0 

415 

416 tmin, tmax, icontrol, tcontrol, ok = analyse_gps_tags( 

417 header, (ipos, t, fix, nsvs), offset, nsamples) 

418 

419 la = num.logical_and 

420 nok = num.logical_not(ok) 

421 

422 axes.plot( 

423 x[la(bfix, ok)]/HOUR, y[la(bfix, ok)], '+', 

424 ms=5, color=color('chameleon3')) 

425 axes.plot( 

426 x[la(bfix, nok)]/HOUR, y[la(bfix, nok)], '+', 

427 ms=5, color=color('aluminium4')) 

428 

429 axes.plot( 

430 x[la(bnofix, ok)]/HOUR, y[la(bnofix, ok)], 'x', 

431 ms=5, color=color('chocolate3')) 

432 axes.plot( 

433 x[la(bnofix, nok)]/HOUR, y[la(bnofix, nok)], 'x', 

434 ms=5, color=color('aluminium4')) 

435 

436 tred = tcontrol - icontrol*deltat - tref 

437 axes.plot(icontrol*deltat/HOUR, tred, color=color('aluminium6')) 

438 axes.plot(icontrol*deltat/HOUR, tred, 'o', ms=5, color=color('aluminium6')) 

439 

440 ymin = (math.floor(tred.min() / deltat)-1) * deltat 

441 ymax = (math.ceil(tred.max() / deltat)+1) * deltat 

442 # ymin = min(ymin, num.min(y)) 

443 # ymax = max(ymax, num.max(y)) 

444 

445 if ymax - ymin < 1000 * deltat: 

446 ygrid = math.floor(tred.min() / deltat) * deltat 

447 while ygrid < ymax: 

448 axes.axhline(ygrid, color=color('aluminium4'), alpha=0.3) 

449 ygrid += deltat 

450 

451 xmin = icontrol[0]*deltat/HOUR 

452 xmax = icontrol[-1]*deltat/HOUR 

453 xsize = xmax - xmin 

454 xmin -= xsize * 0.1 

455 xmax += xsize * 0.1 

456 axes.set_xlim(xmin, xmax) 

457 

458 axes.set_ylim(ymin, ymax) 

459 

460 axes.set_xlabel('Uncorrected (quartz) time [h]') 

461 axes.set_ylabel('Relative time correction [s]') 

462 

463 plt.show() 

464 

465 

466g_dir_contexts = {} 

467 

468 

469class DirContextEntry(Object): 

470 path = String.T() 

471 tstart = Timestamp.T() 

472 ifile = Int.T() 

473 

474 

475class DirContext(Object): 

476 path = String.T() 

477 mtime = Timestamp.T() 

478 entries = DirContextEntry.T() 

479 

480 def get_entry(self, fn): 

481 path = os.path.abspath(fn) 

482 for entry in self.entries: 

483 if entry.path == path: 

484 return entry 

485 

486 raise Exception('entry not found') 

487 

488 def iter_entries(self, fn, step=1): 

489 current = self.get_entry(fn) 

490 by_ifile = dict( 

491 (entry.ifile, entry) for entry in self.entries 

492 if entry.tstart == current.tstart) 

493 

494 icurrent = current.ifile 

495 while True: 

496 icurrent += step 

497 try: 

498 yield by_ifile[icurrent] 

499 

500 except KeyError: 

501 break 

502 

503 

504def context(fn): 

505 from pyrocko import datacube_ext 

506 

507 dpath = os.path.dirname(os.path.abspath(fn)) 

508 mtimes = [os.stat(dpath)[8]] 

509 

510 dentries = sorted([os.path.join(dpath, f) for f in os.listdir(dpath) 

511 if os.path.isfile(os.path.join(dpath, f))]) 

512 for dentry in dentries: 

513 fn2 = os.path.join(dpath, dentry) 

514 mtimes.append(os.stat(fn2)[8]) 

515 

516 mtime = float(max(mtimes)) 

517 

518 if dpath in g_dir_contexts: 

519 dir_context = g_dir_contexts[dpath] 

520 if dir_context.mtime == mtime: 

521 return dir_context 

522 

523 del g_dir_contexts[dpath] 

524 

525 entries = [] 

526 for dentry in dentries: 

527 fn2 = os.path.join(dpath, dentry) 

528 if not os.path.isfile(fn2): 

529 continue 

530 

531 with open(fn2, 'rb') as f: 

532 first512 = f.read(512) 

533 if not detect(first512): 

534 continue 

535 

536 with open(fn2, 'rb') as f: 

537 try: 

538 header, data_arrays, gps_tags, nsamples, _ = \ 

539 datacube_ext.load(f.fileno(), 3, 0, -1, None) 

540 

541 except datacube_ext.DataCubeError as e: 

542 e = DataCubeError(str(e)) 

543 e.set_context('filename', fn) 

544 raise e 

545 

546 header = dict(header) 

547 entries.append(DirContextEntry( 

548 path=os.path.abspath(fn2), 

549 tstart=util.str_to_time( 

550 '20' + header['S_DATE'] + ' ' + header['S_TIME'], 

551 format='%Y/%m/%d %H:%M:%S'), 

552 ifile=int(header['DAT_NO']))) 

553 

554 dir_context = DirContext(mtime=mtime, path=dpath, entries=entries) 

555 

556 return dir_context 

557 

558 

559def get_time_infos(fn): 

560 from pyrocko import datacube_ext 

561 

562 with open(fn, 'rb') as f: 

563 try: 

564 header, _, gps_tags, nsamples, _ = datacube_ext.load( 

565 f.fileno(), 1, 0, -1, None) 

566 

567 except datacube_ext.DataCubeError as e: 

568 e = DataCubeError(str(e)) 

569 e.set_context('filename', fn) 

570 raise e 

571 

572 return dict(header), GPSTags.from_tuple(gps_tags), nsamples 

573 

574 

575def get_timing_context(fns): 

576 joined = [[], [], [], []] 

577 ioff = 0 

578 for fn in fns: 

579 header, gps_tags, nsamples = get_time_infos(fn) 

580 

581 gps_tags = GPSTags.from_tuple(gps_tags) 

582 

583 ipos = gps_tags.position 

584 ipos += ioff 

585 

586 for i in range(4): 

587 joined[i].append(gps_tags[i]) 

588 

589 ioff += nsamples 

590 

591 ipos, t, fix, nsvs = [num.concatenate(x) for x in joined] 

592 

593 nsamples = ioff 

594 return ipos, t, fix, nsvs, header, 0, nsamples 

595 

596 

597def get_extended_timing_context(fn): 

598 c = context(fn) 

599 

600 header, gps_tags, nsamples_base = get_time_infos(fn) 

601 gps_tags = GPSTags.from_tuple(gps_tags) 

602 

603 ioff = 0 

604 aggregated = [gps_tags] 

605 

606 nsamples_total = nsamples_base 

607 

608 if num.sum(gps_tags.fix) == 0: 

609 

610 ioff = nsamples_base 

611 for entry in c.iter_entries(fn, 1): 

612 

613 _, gps_tags, nsamples = get_time_infos(entry.path) 

614 

615 ipos = gps_tags.position 

616 ipos += ioff 

617 

618 aggregated.append(gps_tags) 

619 nsamples_total += nsamples 

620 

621 if num.sum(gps_tags.fix) > 0: 

622 break 

623 

624 ioff += nsamples 

625 

626 ioff = 0 

627 for entry in c.iter_entries(fn, -1): 

628 

629 _, gps_tags, nsamples = get_time_infos(entry.path) 

630 

631 ioff -= nsamples 

632 

633 ipos = gps_tags.position 

634 ipos += ioff 

635 

636 aggregated[0:0] = [gps_tags] 

637 

638 nsamples_total += nsamples 

639 

640 if num.sum(gps_tags.fix) > 0: 

641 break 

642 

643 concat = num.concatenate 

644 

645 ipos = concat([tag.position for tag in aggregated]) 

646 time = concat([tag.gps_time for tag in aggregated]) 

647 fix = concat([tag.fix for tag in aggregated]) 

648 sat_count = concat([tag.sat_count for tag in aggregated]) 

649 gps_utc_offsets = concat([tag.gps_utc_offsets for tag in aggregated]) 

650 

651 return ipos, time, fix, sat_count, gps_utc_offsets, header, \ 

652 0, nsamples_base 

653 

654 

655def iload(fn, load_data=True, interpolation='sinc', yield_gps_tags=False): 

656 from pyrocko import datacube_ext 

657 from pyrocko import signal_ext 

658 

659 if interpolation not in ('sinc', 'off'): 

660 raise NotImplementedError( 

661 'no such interpolation method: %s' % interpolation) 

662 

663 with open(fn, 'rb') as f: 

664 loadflag = 2 if load_data else 1 

665 

666 try: 

667 header, data_arrays, gps_tags, nsamples, _ = datacube_ext.load( 

668 f.fileno(), loadflag, 0, -1, None) 

669 

670 except datacube_ext.DataCubeError as e: 

671 e = DataCubeError(str(e)) 

672 e.set_context('filename', fn) 

673 raise e 

674 

675 header = dict(header) 

676 deltat = 1.0 / int(header['S_RATE']) 

677 nchannels = int(header['CH_NUM']) 

678 

679 ipos, times, fix, sat_count, utc_gps_offsets, \ 

680 header_, offset_, nsamples_ = get_extended_timing_context(fn) 

681 

682 tmin, tmax, icontrol, tcontrol, _ = analyse_gps_tags( 

683 header_, (ipos, times, fix, sat_count), offset_, nsamples_) 

684 

685 if icontrol is None: 

686 warnings.warn( 

687 'No usable GPS timestamps found. Using datacube header ' 

688 'information to guess time. (file: "%s")' % fn) 

689 

690 tmin_ip = round(tmin / deltat) * deltat 

691 if interpolation != 'off': 

692 tmax_ip = round(tmax / deltat) * deltat 

693 else: 

694 tmax_ip = tmin_ip + (nsamples-1) * deltat 

695 

696 nsamples_ip = int(round((tmax_ip - tmin_ip)/deltat)) + 1 

697 # to prevent problems with rounding errors: 

698 tmax_ip = tmin_ip + (nsamples_ip-1) * deltat 

699 

700 leaps = num.array( 

701 [x[0] + util.gps_utc_offset(x[0]) for x in util.read_leap_seconds2()], 

702 dtype=float) 

703 

704 if load_data and icontrol is not None: 

705 ncontrol_this = num.sum( 

706 num.logical_and(0 <= icontrol, icontrol < nsamples)) 

707 

708 if ncontrol_this <= 1: 

709 warnings.warn( 

710 'Extrapolating GPS time information from directory context ' 

711 '(insufficient number of GPS timestamps in file: "%s").' % fn) 

712 

713 for i in range(nchannels): 

714 if load_data: 

715 arr = data_arrays[i] 

716 assert arr.size == nsamples 

717 

718 if interpolation == 'sinc' and icontrol is not None: 

719 

720 ydata = num.empty(nsamples_ip, dtype=float) 

721 try: 

722 signal_ext.antidrift( 

723 icontrol, tcontrol, 

724 arr.astype(float), tmin_ip, deltat, ydata) 

725 

726 except signal_ext.Error as e: 

727 e = DataCubeError(str(e)) 

728 e.set_context('filename', fn) 

729 e.set_context('n_control_points', icontrol.size) 

730 e.set_context('n_samples_raw', arr.size) 

731 e.set_context('n_samples_ip', ydata.size) 

732 e.set_context('tmin_ip', util.time_to_str(tmin_ip)) 

733 raise e 

734 

735 ydata = num.round(ydata).astype(arr.dtype) 

736 else: 

737 ydata = arr 

738 

739 tr_tmin = tmin_ip 

740 tr_tmax = None 

741 else: 

742 ydata = None 

743 tr_tmin = tmin_ip 

744 tr_tmax = tmax_ip 

745 

746 tr = trace.Trace('', header['DEV_NO'], '', 'p%i' % i, deltat=deltat, 

747 ydata=ydata, tmin=tr_tmin, tmax=tr_tmax, meta=header) 

748 bleaps = num.logical_and(tmin_ip <= leaps, leaps < tmax_ip) 

749 

750 if num.any(bleaps):