1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import, division 

6 

7import math 

8import copy 

9import logging 

10import numpy as num 

11 

12from pyrocko import orthodrome 

13from pyrocko.orthodrome import wrap 

14from pyrocko.guts import Object, Float, String, List, dump_all 

15 

16from .location import Location 

17 

18logger = logging.getLogger('pyrocko.model.station') 

19 

20guts_prefix = 'pf' 

21 

22d2r = num.pi / 180. 

23 

24 

25class ChannelsNotOrthogonal(Exception): 

26 pass 

27 

28 

29def guess_azimuth_from_name(channel_name): 

30 if channel_name.endswith('N'): 

31 return 0. 

32 elif channel_name.endswith('E'): 

33 return 90. 

34 elif channel_name.endswith('Z'): 

35 return 0. 

36 

37 return None 

38 

39 

40def guess_dip_from_name(channel_name): 

41 if channel_name.endswith('N'): 

42 return 0. 

43 elif channel_name.endswith('E'): 

44 return 0. 

45 elif channel_name.endswith('Z'): 

46 return -90. 

47 

48 return None 

49 

50 

51def guess_azimuth_dip_from_name(channel_name): 

52 return guess_azimuth_from_name(channel_name), \ 

53 guess_dip_from_name(channel_name) 

54 

55 

56def mkvec(x, y, z): 

57 return num.array([x, y, z], dtype=float) 

58 

59 

60def are_orthogonal(enus, eps=0.05): 

61 return all(abs(x) < eps for x in [ 

62 num.dot(enus[0], enus[1]), 

63 num.dot(enus[1], enus[2]), 

64 num.dot(enus[2], enus[0])]) 

65 

66 

67def fill_orthogonal(enus): 

68 

69 nmiss = sum(x is None for x in enus) 

70 

71 if nmiss == 1: 

72 for ic in range(len(enus)): 

73 if enus[ic] is None: 

74 enus[ic] = num.cross(enus[(ic-2) % 3], enus[(ic-1) % 3]) 

75 

76 if nmiss == 2: 

77 for ic in range(len(enus)): 

78 if enus[ic] is not None: 

79 xenu = enus[ic] + mkvec(1, 1, 1) 

80 enus[(ic+1) % 3] = num.cross(enus[ic], xenu) 

81 enus[(ic+2) % 3] = num.cross(enus[ic], enus[(ic+1) % 3]) 

82 

83 if nmiss == 3: 

84 # holy camoly.. 

85 enus[0] = mkvec(1, 0, 0) 

86 enus[1] = mkvec(0, 1, 0) 

87 enus[2] = mkvec(0, 0, 1) 

88 

89 

90class Channel(Object): 

91 name = String.T() 

92 azimuth = Float.T(optional=True) 

93 dip = Float.T(optional=True) 

94 gain = Float.T(default=1.0) 

95 

96 def __init__(self, name, azimuth=None, dip=None, gain=1.0): 

97 if azimuth is None: 

98 azimuth = guess_azimuth_from_name(name) 

99 if dip is None: 

100 dip = guess_dip_from_name(name) 

101 

102 Object.__init__( 

103 self, 

104 name=name, 

105 azimuth=float_or_none(azimuth), 

106 dip=float_or_none(dip), 

107 gain=float(gain)) 

108 

109 @property 

110 def ned(self): 

111 if None in (self.azimuth, self.dip): 

112 return None 

113 

114 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r) 

115 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r) 

116 d = math.sin(self.dip*d2r) 

117 return mkvec(n, e, d) 

118 

119 @property 

120 def enu(self): 

121 if None in (self.azimuth, self.dip): 

122 return None 

123 

124 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r) 

125 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r) 

126 d = math.sin(self.dip*d2r) 

127 return mkvec(e, n, -d) 

128 

129 def __str__(self): 

130 return '%s %f %f %g' % (self.name, self.azimuth, self.dip, self.gain) 

131 

132 

133class Station(Location): 

134 network = String.T() 

135 station = String.T() 

136 location = String.T() 

137 name = String.T(default='') 

138 channels = List.T(Channel.T()) 

139 

140 def __init__(self, network='', station='', location='', 

141 lat=0.0, lon=0.0, 

142 elevation=0.0, depth=0.0, 

143 north_shift=0.0, east_shift=0.0, 

144 name='', channels=None): 

145 

146 Location.__init__( 

147 self, 

148 network=network, station=station, location=location, 

149 lat=float(lat), lon=float(lon), 

150 elevation=float(elevation), 

151 depth=float(depth), 

152 north_shift=float(north_shift), 

153 east_shift=float(east_shift), 

154 name=name or '', 

155 channels=channels or []) 

156 

157 self.dist_deg = None 

158 self.dist_m = None 

159 self.azimuth = None 

160 self.backazimuth = None 

161 

162 def copy(self): 

163 return copy.deepcopy(self) 

164 

165 def set_event_relative_data(self, event, distance_3d=False): 

166 surface_dist = self.distance_to(event) 

167 if distance_3d: 

168 if event.depth is None: 

169 logger.warn('No event depth given: using 0 m.') 

170 dd = 0.0 - self.depth 

171 else: 

172 dd = event.depth - self.depth 

173 

174 self.dist_m = math.sqrt(dd**2 + surface_dist**2) 

175 else: 

176 self.dist_m = surface_dist 

177 

178 self.dist_deg = surface_dist / orthodrome.earthradius_equator * \ 

179 orthodrome.r2d 

180 

181 self.azimuth, self.backazimuth = event.azibazi_to(self) 

182 

183 def set_channels_by_name(self, *args): 

184 self.set_channels([]) 

185 for name in args: 

186 self.add_channel(Channel(name)) 

187 

188 def set_channels(self, channels): 

189 self.channels = [] 

190 for ch in channels: 

191 self.add_channel(ch) 

192 

193 def get_channels(self): 

194 return list(self.channels) 

195 

196 def get_channel_names(self): 

197 return set(ch.name for ch in self.channels) 

198 

199 def remove_channel_by_name(self, name): 

200 todel = [ch for ch in self.channels if ch.name == name] 

201 for ch in todel: 

202 self.channels.remove(ch) 

203 

204 def add_channel(self, channel): 

205 self.remove_channel_by_name(channel.name) 

206 self.channels.append(channel) 

207 self.channels.sort(key=lambda ch: ch.name) 

208 

209 def get_channel(self, name): 

210 for ch in self.channels: 

211 if ch.name == name: 

212 return ch 

213 

214 return None 

215 

216 def rotation_ne_to_rt(self, in_channel_names, out_channel_names): 

217 

218 angle = wrap(self.backazimuth + 180., -180., 180.) 

219 in_channels = [self.get_channel(name) for name in in_channel_names] 

220 out_channels = [ 

221 Channel(out_channel_names[0], 

222 wrap(self.backazimuth+180., -180., 180.), 0., 1.), 

223 Channel(out_channel_names[1], 

224 wrap(self.backazimuth+270., -180., 180.), 0., 1.)] 

225 return angle, in_channels, out_channels 

226 

227 def _projection_to( 

228 self, to, in_channel_names, out_channel_names, use_gains=False): 

229 

230 in_channels = [self.get_channel(name) for name in in_channel_names] 

231 

232 # create orthogonal vectors for missing components, such that this 

233 # won't break projections when components are missing. 

234 

235 vecs = [] 

236 for ch in in_channels: 

237 if ch is None: 

238 vecs.append(None) 

239 else: 

240 vec = getattr(ch, to) 

241 if use_gains: 

242 vec /= ch.gain 

243 vecs.append(vec) 

244 

245 fill_orthogonal(vecs) 

246 if not are_orthogonal(vecs): 

247 raise ChannelsNotOrthogonal( 

248 'components are not orthogonal: station %s.%s.%s, ' 

249 'channels %s, %s, %s' 

250 % (self.nsl() + tuple(in_channel_names))) 

251 

252 m = num.hstack([vec2[:, num.newaxis] for vec2 in vecs]) 

253 

254 m = num.where(num.abs(m) < num.max(num.abs(m))*1e-16, 0., m) 

255 

256 if to == 'ned': 

257 out_channels = [ 

258 Channel(out_channel_names[0], 0., 0., 1.), 

259 Channel(out_channel_names[1], 90., 0., 1.), 

260 Channel(out_channel_names[2], 0., 90., 1.)] 

261 

262 elif to == 'enu': 

263 out_channels = [ 

264 Channel(out_channel_names[0], 90., 0., 1.), 

265 Channel(out_channel_names[1], 0., 0., 1.), 

266 Channel(out_channel_names[2], 0., -90., 1.)] 

267 

268 return m, in_channels, out_channels 

269 

270 def guess_channel_groups(self): 

271 cg = {} 

272 for channel in self.get_channels(): 

273 if len(channel.name) >= 1: 

274 kind = channel.name[:-1] 

275 if kind not in cg: 

276 cg[kind] = [] 

277 cg[kind].append(channel.name[-1]) 

278 

279 def allin(a, b): 

280 return all(x in b for x in a) 

281 

282 out_groups = [] 

283 for kind, components in cg.items(): 

284 for sys in ('ENZ', '12Z', 'XYZ', 'RTZ', '123'): 

285 if allin(sys, components): 

286 out_groups.append(tuple([kind+c for c in sys])) 

287 

288 return out_groups 

289 

290 def guess_projections_to_enu(self, out_channels=('E', 'N', 'U'), **kwargs): 

291 proj = [] 

292 for cg in self.guess_channel_groups(): 

293 try: 

294 proj.append(self.projection_to_enu( 

295 cg, out_channels=out_channels, **kwargs)) 

296 

297 except ChannelsNotOrthogonal as e: 

298 logger.warning(str(e)) 

299 

300 return proj 

301 

302 def guess_projections_to_rtu( 

303 self, out_channels=('R', 'T', 'U'), backazimuth=None, **kwargs): 

304 

305 if backazimuth is None: 

306 backazimuth = self.backazimuth 

307 out_channels_ = [ 

308 Channel( 

309 out_channels[0], wrap(backazimuth+180., -180., 180.), 0., 1.), 

310 Channel( 

311 out_channels[1], wrap(backazimuth+270., -180., 180.), 0., 1.), 

312 Channel( 

313 out_channels[2], 0., -90., 1.)] 

314 

315 proj = [] 

316 for (m, in_channels, _) in self.guess_projections_to_enu(**kwargs): 

317 phi = (backazimuth + 180.)*d2r 

318 r = num.array([[math.sin(phi), math.cos(phi), 0.0], 

319 [math.cos(phi), -math.sin(phi), 0.0], 

320 [0.0, 0.0, 1.0]]) 

321 proj.append((num.dot(r, m), in_channels, out_channels_)) 

322 

323 return proj 

324 

325 def projection_to_enu( 

326 self, 

327 in_channels, 

328 out_channels=('E', 'N', 'U'), 

329 **kwargs): 

330 

331 return self._projection_to('enu', in_channels, out_channels, **kwargs) 

332 

333 def projection_to_ned( 

334 self, 

335 in_channels, 

336 out_channels=('N', 'E', 'D'), 

337 **kwargs): 

338 

339 return self._projection_to('ned', in_channels, out_channels, **kwargs) 

340 

341 def projection_from_enu( 

342 self, 

343 in_channels=('E', 'N', 'U'), 

344 out_channels=('X', 'Y', 'Z'), 

345 **kwargs): 

346 

347 m, out_channels, in_channels = self._projection_to( 

348 'enu', out_channels, in_channels, **kwargs) 

349 

350 return num.linalg.inv(m), in_channels, out_channels 

351 

352 def projection_from_ned( 

353 self, 

354 in_channels=('N', 'E', 'D'), 

355 out_channels=('X', 'Y', 'Z'), 

356 **kwargs): 

357 

358 m, out_channels, in_channels = self._projection_to( 

359 'ned', out_channels, in_channels, **kwargs) 

360 

361 return num.linalg.inv(m), in_channels, out_channels 

362 

363 def nsl_string(self): 

364 return '.'.join((self.network, self.station, self.location)) 

365 

366 def nsl(self): 

367 return self.network, self.station, self.location 

368 

369 def cannot_handle_offsets(self): 

370 if self.north_shift != 0.0 or self.east_shift != 0.0: 

371 logger.warn( 

372 'Station %s.%s.%s has a non-zero Cartesian offset. Such ' 

373 'offsets cannot be saved in the basic station file format. ' 

374 'Effective lat/lons are saved only. Please save the stations ' 

375 'in YAML format to preserve the reference-and-offset ' 

376 'coordinates.' % self.nsl()) 

377 

378 def oldstr(self): 

379 self.cannot_handle_offsets() 

380 nsl = '%s.%s.%s' % (self.network, self.station, self.location) 

381 s = '%-15s %14.5f %14.5f %14.1f %14.1f %s' % ( 

382 nsl, self.effective_lat, self.effective_lon, self.elevation, 

383 self.depth, self.name) 

384 return s 

385 

386 

387def dump_stations(stations, filename): 

388 ''' 

389 Write stations file. 

390 

391 :param stations: list of :py:class:`Station` objects 

392 :param filename: filename as str 

393 ''' 

394 f = open(filename, 'w') 

395 for sta in stations: 

396 f.write(sta.oldstr()+'\n') 

397 for cha in sta.get_channels(): 

398 azimuth = 'NaN' 

399 if cha.azimuth is not None: 

400 azimuth = '%7g' % cha.azimuth 

401 

402 dip = 'NaN' 

403 if cha.dip is not None: 

404 dip = '%7g' % cha.dip 

405 

406 f.write('%5s %14s %14s %14g\n' % ( 

407 cha.name, azimuth, dip, cha.gain)) 

408 

409 f.close() 

410 

411 

412def dump_stations_yaml(stations, filename): 

413 ''' 

414 Write stations file in YAML format. 

415 

416 :param stations: list of :py:class:`Station` objects 

417 :param filename: filename as str 

418 ''' 

419 

420 dump_all(stations, filename=filename) 

421 

422 

423def float_or_none(s): 

424 if s is None: 

425 return None 

426 elif isinstance(s, str) and s.lower() == 'nan': 

427 return None 

428 else: 

429 return float(s) 

430 

431 

432def detect_format(filename): 

433 with open(filename, 'r') as f: 

434 for line in f: 

435 line = line.strip() 

436 if not line or line.startswith('#') or line.startswith('%'): 

437 continue 

438 if line.startswith('--- !pf.Station'): 

439 return 'yaml' 

440 else: 

441 return 'basic' 

442 

443 return 'basic' 

444 

445 

446def load_stations(filename, format='detect'): 

447 ''' 

448 Read stations file. 

449 

450 :param filename: filename 

451 :returns: list of :py:class:`Station` objects 

452 ''' 

453 

454 if format == 'detect': 

455 format = detect_format(filename) 

456 

457 if format == 'yaml': 

458 from pyrocko import guts 

459 stations = [ 

460 st for st in guts.load_all(filename=filename) 

461 if isinstance(st, Station)] 

462 

463 return stations 

464 

465 elif format == 'basic': 

466 stations = [] 

467 f = open(filename, 'r') 

468 station = None 

469 channel_names = [] 

470 for (iline, line) in enumerate(f): 

471 toks = line.split(None, 5) 

472 if line.strip().startswith('#') or line.strip() == '': 

473 continue 

474 

475 if len(toks) == 5 or len(toks) == 6: 

476 net, sta, loc = toks[0].split('.') 

477 lat, lon, elevation, depth = [float(x) for x in toks[1:5]] 

478 if len(toks) == 5: 

479 name = '' 

480 else: 

481 name = toks[5].rstrip() 

482 

483 station = Station( 

484 net, sta, loc, lat, lon, 

485 elevation=elevation, depth=depth, name=name) 

486 

487 stations.append(station) 

488 channel_names = [] 

489 

490 elif len(toks) == 4 and station is not None: 

491 name, azi, dip, gain = ( 

492 toks[0], 

493 float_or_none(toks[1]), 

494 float_or_none(toks[2]), 

495 float(toks[3])) 

496 if name in channel_names: 

497 logger.warning( 

498 'redefined channel! (line: %i, file: %s)' % 

499 (iline + 1, filename)) 

500 else: 

501 channel_names.append(name) 

502 

503 channel = Channel(name, azimuth=azi, dip=dip, gain=gain) 

504 station.add_channel(channel) 

505 

506 else: 

507 logger.warning('skipping invalid station/channel definition ' 

508 '(line: %i, file: %s' % (iline + 1, filename)) 

509 

510 f.close() 

511 return stations 

512 

513 else: 

514 from pyrocko.io.io_common import FileLoadError 

515 raise FileLoadError('unknown event file format: %s' % format) 

516 

517 

518def dump_kml(objects, filename): 

519 station_template = ''' 

520 <Placemark> 

521 <name>%(network)s.%(station)s.%(location)s</name> 

522 <description></description> 

523 <styleUrl>#msn_S</styleUrl> 

524 <Point> 

525 <coordinates>%(elon)f,%(elat)f,%(elevation)f</coordinates> 

526 </Point> 

527 </Placemark> 

528''' 

529 

530 f = open(filename, 'w') 

531 f.write('<?xml version="1.0" encoding="UTF-8"?>\n') 

532 f.write('<kml xmlns="http://www.opengis.net/kml/2.2">\n') 

533 f.write('<Document>\n') 

534 f.write(''' <Style id="sh_S"> 

535 <IconStyle> 

536 <scale>1.3</scale> 

537 <Icon> 

538 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href> 

539 </Icon> 

540 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/> 

541 </IconStyle> 

542 <ListStyle> 

543 <ItemIcon> 

544 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href> 

545 </ItemIcon> 

546 </ListStyle> 

547 </Style> 

548 <Style id="sn_S"> 

549 <IconStyle> 

550 <scale>1.1</scale> 

551 <Icon> 

552 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href> 

553 </Icon> 

554 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/> 

555 </IconStyle> 

556 <ListStyle> 

557 <ItemIcon> 

558 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href> 

559 </ItemIcon> 

560 </ListStyle> 

561 </Style> 

562 <StyleMap id="msn_S"> 

563 <Pair> 

564 <key>normal</key> 

565 <styleUrl>#sn_S</styleUrl> 

566 </Pair> 

567 <Pair> 

568 <key>highlight</key> 

569 <styleUrl>#sh_S</styleUrl> 

570 </Pair> 

571 </StyleMap> 

572''') 

573 for obj in objects: 

574 

575 if isinstance(obj, Station): 

576 d = obj.__dict__.copy() 

577 d['elat'], d['elon'] = obj.effective_latlon 

578 f.write(station_template % d) 

579 

580 f.write('</Document>') 

581 f.write('</kml>\n') 

582 f.close()