|
| 1 | +import pathlib |
| 2 | + |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +import numpy as np |
| 5 | +import obspy |
| 6 | +from obspy.core.inventory.response import Response |
| 7 | + |
| 8 | +CACHE_PATH = pathlib.Path("./cache") |
| 9 | +DATA_PATH = pathlib.Path("./data") |
| 10 | + |
| 11 | +# Number of frequencies to test at. |
| 12 | +N_FREQUENCIES = 100 |
| 13 | + |
| 14 | +# TOLERANCES |
| 15 | +RTOL = 1e-4 |
| 16 | +ATOL_AS_FRAC_OF_ABS_MAX = 5e-4 |
| 17 | + |
| 18 | +# Collect phase responses here that don't compare favourable to evalresp |
| 19 | +# but which have been manually verified to be at least as correct in the |
| 20 | +# new scipy implementation. |
| 21 | +# This will have to be expanded in the course of this test. |
| 22 | +SKIP_VALIDATING_PHASE_REPONSE = ["IU.AFI..UHE", "IU.AFI..UHN", "IU.AFI..UHZ"] |
| 23 | + |
| 24 | + |
| 25 | +def compare_single_response(channel_id: str, response: Response): |
| 26 | + # detect sampling rate from response stages |
| 27 | + for stage in response.response_stages[::-1]: |
| 28 | + if ( |
| 29 | + stage.decimation_input_sample_rate is not None |
| 30 | + and stage.decimation_factor is not None |
| 31 | + ): |
| 32 | + sampling_rate = stage.decimation_input_sample_rate / stage.decimation_factor |
| 33 | + break |
| 34 | + else: |
| 35 | + # XXX: Something has to be done here. |
| 36 | + msg = "Failed to autodetect sampling rate of channel from " "response stages." |
| 37 | + raise Exception(msg) |
| 38 | + |
| 39 | + # Compute up to the Nyquist frequency - evalresp's phase usually goes crazy |
| 40 | + # afterwards. |
| 41 | + FREQUENCIES = np.logspace(-2, np.log10(0.5 * sampling_rate), N_FREQUENCIES) |
| 42 | + |
| 43 | + # Compute for evalresp as well as the scipy response. |
| 44 | + try: |
| 45 | + eval_resp = response.get_evalresp_response_for_frequencies( |
| 46 | + FREQUENCIES, output="VEL" |
| 47 | + ) |
| 48 | + except Exception as e: |
| 49 | + print( |
| 50 | + " evalresp failed to compute response. Thus no comparision can be" |
| 51 | + f"performed. Reason for evalresp failure: {str(e)}" |
| 52 | + ) |
| 53 | + return |
| 54 | + |
| 55 | + scipy_resp = response.get_response(FREQUENCIES, output="VEL") |
| 56 | + |
| 57 | + # Use amplitude and phase for the comparison just because it is more |
| 58 | + # intuitive. |
| 59 | + eval_resp_amplitude = np.abs(eval_resp) |
| 60 | + eval_resp_phase = np.angle(eval_resp) |
| 61 | + |
| 62 | + scipy_resp_amplitude = np.abs(scipy_resp) |
| 63 | + scipy_resp_phase = np.angle(scipy_resp) |
| 64 | + |
| 65 | + atol_amplitude = scipy_resp_amplitude.max() * ATOL_AS_FRAC_OF_ABS_MAX |
| 66 | + atol_phase = np.abs(scipy_resp_phase).max() * ATOL_AS_FRAC_OF_ABS_MAX |
| 67 | + |
| 68 | + try: |
| 69 | + np.testing.assert_allclose( |
| 70 | + eval_resp_amplitude, |
| 71 | + scipy_resp_amplitude, |
| 72 | + rtol=RTOL, |
| 73 | + atol=atol_amplitude, |
| 74 | + err_msg="amplitude mismatch", |
| 75 | + ) |
| 76 | + # Skip if manually verified. |
| 77 | + if channel_id not in SKIP_VALIDATING_PHASE_REPONSE: |
| 78 | + np.testing.assert_allclose( |
| 79 | + eval_resp_phase, |
| 80 | + scipy_resp_phase, |
| 81 | + rtol=RTOL, |
| 82 | + atol=atol_phase, |
| 83 | + err_msg="amplitude mismatch", |
| 84 | + ) |
| 85 | + except Exception as e: |
| 86 | + print(f"Failed comparison due to: {str(e)}") |
| 87 | + print("Will now produce a plot to help diagnose the issue.") |
| 88 | + |
| 89 | + plt.subplot(411) |
| 90 | + plt.title("Amplitude response") |
| 91 | + plt.loglog(FREQUENCIES, eval_resp_amplitude, label="Evalresp") |
| 92 | + plt.loglog(FREQUENCIES, scipy_resp_amplitude, label="scipy") |
| 93 | + plt.legend() |
| 94 | + |
| 95 | + plt.subplot(412) |
| 96 | + plt.title("Amplitude response difference") |
| 97 | + plt.loglog(FREQUENCIES, eval_resp_amplitude - scipy_resp_amplitude) |
| 98 | + |
| 99 | + plt.subplot(413) |
| 100 | + plt.title("Phase response") |
| 101 | + plt.semilogx(FREQUENCIES, eval_resp_phase, label="Evalresp") |
| 102 | + plt.semilogx(FREQUENCIES, scipy_resp_phase, label="scipy") |
| 103 | + plt.legend() |
| 104 | + |
| 105 | + plt.title("Phase response difference") |
| 106 | + plt.subplot(414) |
| 107 | + plt.semilogx(FREQUENCIES, eval_resp_phase - scipy_resp_phase) |
| 108 | + |
| 109 | + plt.show() |
| 110 | + raise e |
| 111 | + |
| 112 | + |
| 113 | +def test_single_stationxml_file(filename: pathlib.Path): |
| 114 | + def _p(msg, indent: int = 0): |
| 115 | + print(f"{' ' * indent}File '{filename}': {msg}") |
| 116 | + |
| 117 | + # Simplistic cache to be able to rerun this a lot and fix |
| 118 | + # bugs as they appear. |
| 119 | + cache_file = CACHE_PATH / filename.name |
| 120 | + if cache_file.exists(): |
| 121 | + _p("Already has been tested. Skipping ...", indent=2) |
| 122 | + return |
| 123 | + |
| 124 | + try: |
| 125 | + inv = obspy.read_inventory(str(filename)) |
| 126 | + except Exception as e: |
| 127 | + _p(f"Failed to parse due to: {str(e)}") |
| 128 | + raise e |
| 129 | + |
| 130 | + all_responses = [] |
| 131 | + for net in inv: |
| 132 | + for sta in net: |
| 133 | + for cha in sta: |
| 134 | + all_responses.append( |
| 135 | + [ |
| 136 | + net.code, |
| 137 | + sta.code, |
| 138 | + cha.location_code, |
| 139 | + cha.code, |
| 140 | + cha.start_date, |
| 141 | + cha.end_date, |
| 142 | + cha.response, |
| 143 | + ] |
| 144 | + ) |
| 145 | + for _i, c in enumerate(all_responses): |
| 146 | + _p(f"Comparing responses for {c[:-1]} ...", indent=2) |
| 147 | + compare_single_response(channel_id=".".join(c[:4]), response=c[-1]) |
| 148 | + |
| 149 | + # Finally just touch the cache file so it will be skipped the next run. |
| 150 | + cache_file.touch() |
| 151 | + |
| 152 | + |
| 153 | +def main(): |
| 154 | + CACHE_PATH.mkdir(exist_ok=True) |
| 155 | + all_files = list(DATA_PATH.glob("*.xml")) |
| 156 | + for _i, filename in enumerate(all_files): |
| 157 | + print(f"Reading StationXML file {_i + 1} of {len(all_files)}: {filename}") |
| 158 | + test_single_stationxml_file(filename) |
| 159 | + |
| 160 | + |
| 161 | +if __name__ == "__main__": |
| 162 | + main() |
0 commit comments