Coverage for C:\Users\peaco\OneDrive\Documents\GitHub\mth5\mth5\utils\fdsn

1# -*- coding: utf-8 -*-

2"""

3FDSN standards tools.

5Tools for working with FDSN (Incorporated Research Institutions for Seismology)

6standards including SEED channel codes, period/measurement/orientation codes,

7and conversions between SEED and MT (magnetotelluric) channel formats.

9Notes

10-----

11Created on Wed Sep 30 11:47:01 2020

13Author: Jared Peacock

15License: MIT

17References

18----------

19FDSN Channel Codes: https://www.fdsn.org/seed_manual/SEEDManual_V2.4.pdf

20"""

22from __future__ import annotations

24from loguru import logger

27period_code_dict = {

28 "F": {"min": 1000, "max": 5000},

29 "G": {"min": 1000, "max": 5000},

30 "D": {"min": 250, "max": 1000},

31 "C": {"min": 250, "max": 1000},

32 "E": {"min": 80, "max": 250},

33 "S": {"min": 10, "max": 80},

34 "H": {"min": 80, "max": 250},

35 "B": {"min": 10, "max": 80},

36 "M": {"min": 1, "max": 10},

37 "L": {"min": 0.95, "max": 1.05},

38 "V": {"min": 0.095, "max": 0.105},

39 "U": {"min": 0.0095, "max": 0.0105},

40 "R": {"min": 0.0001, "max": 0.001},

41 "P": {"min": 0.00001, "max": 0.0001},

42 "T": {"min": 0.000001, "max": 0.00001},

43 "Q": {"min": 0, "max": 0.000001},

44}

46measurement_code_dict = {

47 "tilt": "A",

48 "creep": "B",

49 "calibration": "C",

50 "pressure": "D",

51 "magnetics": "F",

52 "gravity": "G",

53 "humidity": "I",

54 "temperature": "K",

55 "water_current": "O",

56 "electric": "Q",

57 "rain_fall": "R",

58 "linear_strain": "S",

59 "tide": "T",

60 "wind": "W",

61}

63measurement_code_dict_reverse = dict([(v, k) for k, v in measurement_code_dict.items()])

64# Add Y as fallback for unknown/auxiliary measurements

65measurement_code_dict_reverse["Y"] = "auxiliary"

67orientation_code_dict = {

68 "N": {"min": 0, "max": 15},

69 "E": {"min": 75, "max": 90},

70 "Z": {"min": 0, "max": 15},

71 "1": {"min": 15, "max": 45},

72 "2": {"min": 45, "max": 75},

73 "3": {"min": 15, "max": 75},

74}

76mt_code_dict = {"magnetics": "h", "electric": "e"}

79def get_location_code(channel_obj: object) -> str:

80 """

81 Generate FDSN location code from channel metadata.

83 Creates a 2-character location code from the channel's component

84 and channel number.

86 Parameters

87 ----------

88 channel_obj : object

89 Channel metadata object with `component` and `channel_number` attributes.

90 Expected to be of type `~mt_metadata.timeseries.Channel`.

92 Returns

93 -------

94 str

95 2-character location code formatted as first letter of component

96 and last digit of channel number (e.g., 'E1', 'H0').

98 Examples

99 --------

100 Generate location code for an electric component on channel 5::

101

102 >>> class MockChannel:

103 ... component = 'ex'

104 ... channel_number = 5

105 >>> get_location_code(MockChannel())

106 'E5'

107

108 Magnetic component on channel 12 (wraps to 2)::

109

110 >>> class MockChannel:

111 ... component = 'hx'

112 ... channel_number = 12

113 >>> get_location_code(MockChannel())

114 'H2'

115 """

116 # Type narrowing with duck typing for Mock compatibility

117 assert hasattr(channel_obj, "component") and hasattr(channel_obj, "channel_number")

118

119 location_code = "{0}{1}".format(

120 channel_obj.component[0].upper(), # type: ignore

121 channel_obj.channel_number % 10, # type: ignore

122 )

123

124 return location_code

125

126

127def get_period_code(sample_rate: float) -> str:

128 """

129 Get SEED sampling rate code from sample rate.

130

131 Determines the appropriate FDSN/SEED period code based on the sample

132 rate in samples per second. Codes range from 'Q' (highest frequency,

133 period < 1 μs) to 'F' (lowest frequency, period 1000-5000 s).

134

135 Parameters

136 ----------

137 sample_rate : float

138 Sample rate in samples per second.

139

140 Returns

141 -------

142 str

143 Single character SEED sampling rate code. Defaults to 'A' if no

144 code matches the sample rate.

145

146 Notes

147 -----

148 Code mapping (frequency/period ranges):

149 - 'F', 'G': 1-5 kHz

150 - 'D', 'C': 250-1000 Hz

151 - 'E', 'H': 80-250 Hz

152 - 'S', 'B': 10-80 Hz

153 - 'M': 1-10 Hz

154 - 'L': 0.95-1.05 Hz

155 - 'V': 0.095-0.105 Hz

156 - 'U': 0.0095-0.0105 Hz

157 - 'R': 0.0001-0.001 Hz

158 - 'P': 0.00001-0.0001 Hz

159 - 'T': 0.000001-0.00001 Hz

160 - 'Q': < 0.000001 Hz

161

162 Examples

163 --------

164 Get code for 100 Hz sample rate::

165

166 >>> get_period_code(100.0)

167 'B'

168

169 Get code for 1000 Hz::

170

171 >>> get_period_code(1000.0)

172 'D'

173

174 Get code for 10 Hz (default 'A')::

175

176 >>> get_period_code(10.0)

177 'M'

178 """

179 period_code = "A"

180 for key, v_dict in sorted(period_code_dict.items()):

181 if (sample_rate >= v_dict["min"]) and (sample_rate <= v_dict["max"]):

182 period_code = key

183 break

184 return period_code

185

186

187def get_measurement_code(measurement: str) -> str:

188 """

189 Get SEED sensor code from measurement type.

190

191 Maps measurement types to single-character SEED sensor codes.

192 Performs case-insensitive substring matching.

193

194 Parameters

195 ----------

196 measurement : str

197 Measurement type (e.g., 'electric', 'magnetics', 'temperature',

198 'tilt', 'pressure', 'humidity', 'gravity', 'calibration',

199 'rain_fall', 'water_current', 'wind', 'linear_strain', 'tide',

200 'creep').

201

202 Returns

203 -------

204 str

205 Single character SEED sensor code. Returns 'Y' if measurement

206 type not found in mapping dictionary.

207

208 Notes

209 -----

210 Measurement to code mapping:

211 - 'tilt' → 'A'

212 - 'creep' → 'B'

213 - 'calibration' → 'C'

214 - 'pressure' → 'D'

215 - 'magnetics' → 'F'

216 - 'gravity' → 'G'

217 - 'humidity' → 'I'

218 - 'temperature' → 'K'

219 - 'water_current' → 'O'

220 - 'electric' → 'Q'

221 - 'rain_fall' → 'R'

222 - 'linear_strain' → 'S'

223 - 'tide' → 'T'

224 - 'wind' → 'W'

225 - unknown/auxiliary → 'Y'

226

227 Examples

228 --------

229 Get code for electric measurement::

230

231 >>> get_measurement_code('electric')

232 'Q'

233

234 Get code for magnetic field::

235

236 >>> get_measurement_code('magnetics')

237 'F'

238

239 Unknown measurement returns 'Y'::

240

241 >>> get_measurement_code('unknown')

242 'Y'

243 """

244 sensor_code = "Y"

245 for key, code in measurement_code_dict.items():

246 if measurement.lower() in key:

247 sensor_code = code

248 return sensor_code

249

250

251def get_orientation_code(azimuth: float, orientation: str = "horizontal") -> str:

252 """

253 Get SEED orientation code from azimuth and orientation type.

254

255 Maps azimuth angle to SEED orientation code based on whether the

256 sensor is oriented horizontally or vertically.

257

258 Parameters

259 ----------

260 azimuth : float

261 Azimuth angle in degrees where 0 is north, 90 is east,

262 180 is south, 270 is west. For vertical orientation,

263 0 = vertical down.

264 orientation : {'horizontal', 'vertical'}, default 'horizontal'

265 Type of sensor orientation.

266

267 Returns

268 -------

269 str

270 Single character SEED orientation code.

271

272 Raises

273 ------

274 ValueError

275 If `orientation` is not 'horizontal' or 'vertical'.

276

277 Notes

278 -----

279 Horizontal orientation codes (azimuths):

280 - 'N': 0-15° (North)

281 - 'E': 75-90° (East)

282 - '1': 15-45° (NE quadrant)

283 - '2': 45-75° (SE quadrant)

284

285 Vertical orientation codes:

286 - 'Z': 0-15° (Primary vertical)

287 - '3': 15-75° (Alternate vertical)

288

289 Examples

290 --------

291 Get code for northerly azimuth::

292

293 >>> get_orientation_code(10.0, orientation='horizontal')

294 'N'

295

296 Get code for easterly azimuth::

297

298 >>> get_orientation_code(85.0, orientation='horizontal')

299 'E'

300

301 Get code for vertical sensor::

302

303 >>> get_orientation_code(0.0, orientation='vertical')

304 'Z'

305 """

306 orientation_code = "1"

307 horizontal_keys = ["N", "E", "1", "2"]

308 vertical_keys = ["Z", "3"]

309

310 azimuth = abs(azimuth) % 91

311 if orientation == "horizontal":

312 test_keys = horizontal_keys

313 elif orientation == "vertical":

314 test_keys = vertical_keys

315 else:

316 raise ValueError(

317 f"{orientation} not supported must be [ 'horizontal' | 'vertical' ]"

318 )

319 for key in test_keys:

320 angle_min = orientation_code_dict[key]["min"]

321 angle_max = orientation_code_dict[key]["max"]

322 if (azimuth <= angle_max) and (azimuth >= angle_min):

323 orientation_code = key

324 break

325 return orientation_code

326

327

328def make_channel_code(channel_obj: object) -> str:

329 """

330 Generate 3-character SEED channel code from channel metadata.

331

332 Combines period code, measurement code, and orientation code into

333 a standard 3-character FDSN channel code.

334

335 Parameters

336 ----------

337 channel_obj : object

338 Channel metadata object with attributes: `sample_rate`, `component`,

339 `type`, `measurement_azimuth`, `measurement_tilt`.

340 Expected to be of type `~mt_metadata.timeseries.Channel`.

341

342 Returns

343 -------

344 str

345 3-character SEED channel code (e.g., 'BHZ', 'HHE').

346

347 Notes

348 -----

349 The channel code format is: [Period Code][Measurement Code][Orientation Code]

350

351 - Period code: based on sample_rate

352 - Measurement code: derived from component, with fallback to type

353 - Orientation code: depends on whether component is vertical ('z')

354

355 Examples

356 --------

357 Create channel code for horizontal electric component::

358

359 >>> class MockChannel:

360 ... sample_rate = 100.0

361 ... component = 'ex'

362 ... type = 'electric'

363 ... measurement_azimuth = 0.0

364 ... measurement_tilt = 0.0

365 >>> make_channel_code(MockChannel())

366 'BQN'

367

368 Create channel code for vertical magnetic component::

369

370 >>> class MockChannel:

371 ... sample_rate = 100.0

372 ... component = 'hz'

373 ... type = 'magnetic'

374 ... measurement_azimuth = 0.0

375 ... measurement_tilt = 0.0

376 >>> make_channel_code(MockChannel())

377 'BFZ'

378 """

379 # Type narrowing with duck typing for Mock compatibility

380 assert (

381 hasattr(channel_obj, "sample_rate")

382 and hasattr(channel_obj, "component")

383 and hasattr(channel_obj, "type")

384 and hasattr(channel_obj, "measurement_azimuth")

385 and hasattr(channel_obj, "measurement_tilt")

386 )

387

388 period_code = get_period_code(channel_obj.sample_rate) # type: ignore

389 # Try to get measurement code from component first, then fallback to type

390 sensor_code = get_measurement_code(channel_obj.component) # type: ignore

391 if sensor_code == "Y": # If component didn't match, try type

392 sensor_code = get_measurement_code(channel_obj.type) # type: ignore

393 if "z" in channel_obj.component.lower(): # type: ignore

394 orientation_code = get_orientation_code(

395 channel_obj.measurement_tilt, orientation="vertical" # type: ignore

396 )

397 else:

398 orientation_code = get_orientation_code(channel_obj.measurement_azimuth) # type: ignore

399 channel_code = "{0}{1}{2}".format(period_code, sensor_code, orientation_code)

400

401 return channel_code

402

403

404def read_channel_code(channel_code: str) -> dict[str, dict[str, int] | str | bool]:

405 """

406 Parse FDSN channel code into components.

407

408 Decodes a 3-character SEED channel code into its constituent parts:

409 period range, component type, orientation range, and vertical flag.

410

411 Parameters

412 ----------

413 channel_code : str

414 3-character FDSN channel code (e.g., 'BHZ', 'HHE').

415

416 Returns

417 -------

418 dict

419 Dictionary with keys:

420 - 'period' (dict): Period range with 'min' and 'max' keys (Hz).

421 - 'component' (str): Component type (e.g., 'electric', 'magnetics').

422 - 'orientation' (dict): Angle range with 'min' and 'max' keys (degrees).

423 - 'vertical' (bool): True if component is vertical, False otherwise.

424

425 Raises

426 ------

427 ValueError

428 If channel code is not 3 characters, contains invalid period code,

429 or contains invalid orientation code.

430

431 Notes

432 -----

433 Vertical components are identified by orientation codes 'Z' or '3'.

434

435 Examples

436 --------

437 Decode a horizontal channel code::

438

439 >>> result = read_channel_code('BHE')

440 >>> result['component']

441 'magnetics'

442 >>> result['vertical']

443 False

444

445 Decode a vertical channel code::

446

447 >>> result = read_channel_code('BHZ')

448 >>> result['vertical']

449 True

450 """

451

452 if len(channel_code) != 3:

453 msg = "Input FDSN channel code is not proper format, should be 3 letters"

454 logger.error(msg)

455 raise ValueError(msg)

456 try:

457 period_range = period_code_dict[channel_code[0].upper()]

458 except KeyError:

459 msg = (

460 f"Could not find period range for {channel_code[0]}. ",

461 "Setting to 1",

462 )

463 period_range = {"min": 1, "max": 1}

464 try:

465 component = measurement_code_dict_reverse[channel_code[1].upper()]

466 except KeyError:

467 msg = f"Could not find component for {channel_code[1]}"

468 logger.error(msg)

469 raise ValueError(msg)

470 vertical = False

471 try:

472 orientation = orientation_code_dict[channel_code[2].upper()]

473 if channel_code[2].upper() in ["3", "Z"]:

474 vertical = True

475 except KeyError:

476 msg = (

477 f"Could not find orientation for {channel_code[2]}. ",

478 "Setting to 0.",

479 )

480 logger.error(msg)

481 raise ValueError(msg)

482 return {

483 "period": period_range,

484 "component": component,

485 "orientation": orientation,

486 "vertical": vertical,

487 }

488

489

490def make_mt_channel(

491 code_dict: dict[str, dict[str, int] | str | bool], angle_tol: int = 15

492) -> str:

493 """

494 Convert FDSN code dictionary to magnetotelluric (MT) channel code.

495

496 Maps FDSN codes to MT channel naming convention (e.g., 'ex', 'hy', 'hz').

497

498 Parameters

499 ----------

500 code_dict : dict

501 Dictionary with keys:

502 - 'component' (str): Measurement type (e.g., 'magnetics', 'electric').

503 - 'vertical' (bool): True for vertical, False for horizontal.

504 - 'orientation' (dict): Orientation range with 'min' and 'max'.

505 angle_tol : int, default 15

506 Angle tolerance in degrees for determining cardinal directions.

507

508 Returns

509 -------

510 str

511 2-character MT channel code (e.g., 'ex', 'hy', 'hz').

512 Format: [component code][direction code]

513 - Component: 'e' (electric) or 'h' (magnetic)

514 - Direction: 'x', 'y', 'z' for cardinal or '1', '2', '3' for intermediate

515

516 Notes

517 -----

518 Direction mapping for horizontal channels (0-90°):

519 - 'x': North direction (0-15°)

520 - '1': NE quadrant (15-45°)

521 - 'y': East direction (90-angle_tol to 90°)

522 - '2': SE quadrant (45-90-angle_tol°)

523

524 Vertical channels:

525 - 'z': Primary vertical (0-15°)

526 - '3': Alternate vertical (15-90°)

527

528 Examples

529 --------

530 Create north-oriented electric channel::

531

532 >>> code_dict = {

533 ... 'component': 'electric',

534 ... 'vertical': False,

535 ... 'orientation': {'min': 0, 'max': 15}

536 ... }

537 >>> make_mt_channel(code_dict)

538 'ex'

539

540 Create vertical magnetic channel::

541

542 >>> code_dict = {

543 ... 'component': 'magnetics',

544 ... 'vertical': True,

545 ... 'orientation': {'min': 0, 'max': 15}

546 ... }

547 >>> make_mt_channel(code_dict)

548 'hz'

549 """

550

551 # Type narrowing: extract and validate component and orientation

552 component = code_dict["component"]

553 assert isinstance(component, str), "Component must be a string"

554

555 orientation = code_dict["orientation"]

556 assert isinstance(orientation, dict), "Orientation must be a dict"

557

558 vertical = code_dict["vertical"]

559 assert isinstance(vertical, bool), "Vertical flag must be a bool"

560

561 try:

562 mt_comp = mt_code_dict[component]

563 except KeyError:

564 mt_comp = component

565

566 mt_dir: str = "z" # Default direction

567

568 if not vertical:

569 if orientation["min"] >= 0 and orientation["max"] <= angle_tol:

570 mt_dir = "x"

571 elif orientation["min"] >= angle_tol and orientation["max"] <= 45:

572 mt_dir = "1"

573 elif orientation["min"] >= (90 - angle_tol) and orientation["max"] <= 90:

574 mt_dir = "y"

575 elif orientation["min"] >= 45 and orientation["max"] <= (90 - angle_tol):

576 mt_dir = "2"

577 else:

578 if orientation["min"] >= 0 and orientation["max"] <= angle_tol:

579 mt_dir = "z"

580 elif orientation["min"] >= angle_tol and orientation["max"] <= 90:

581 mt_dir = "3"

582

583 mt_code = f"{mt_comp}{mt_dir}"

584

585 return mt_code

Coverage for C: \ Users \ peaco \ OneDrive \ Documents \ GitHub \ mth5 \ mth5 \ utils \ fdsn_tools.py: 82%

106 statements