Coverage for C: \ Users \ peaco \ OneDrive \ Documents \ GitHub \ mt_metadata \ mt_metadata \ timeseries \ filters \ plotting_helpers.py: 74%
101 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-10 00:11 -0800
1import matplotlib.pyplot as plt
2import numpy as np
3from matplotlib.gridspec import GridSpec
6def is_flat_amplitude(array):
7 """
8 Check of an amplitude response is basically flat. If so, it is best to tune the y-axis lims to
9 make numeric noise invisible
11 Parameters
12 ----------
13 array
15 Returns
16 -------
18 """
19 differences = np.diff(np.abs(array))
20 if np.isclose(differences, 0.0).all():
21 return True
22 else:
23 return False
26def cast_angular_frequency_to_period_or_hertz(w, units):
27 if units.lower() == "period":
28 x_axis = (2.0 * np.pi) / w
29 elif units.lower() == "frequency":
30 x_axis = w / (2.0 * np.pi)
31 return x_axis
34def plot_response(
35 frequencies,
36 complex_response,
37 poles=None,
38 zeros=None,
39 xlim=None,
40 title=None,
41 x_units="Period",
42 unwrap=True,
43 pass_band=None,
44 label=None,
45 normalization_frequency=None,
46):
47 """
48 This function was contributed by Ben Murphy at USGS
49 2021-03-17: there are some issues encountered when using this function to plot generic resposnes, looks
50 like the x-axis gets out of order when using frequency as the axis ...
52 Edited 2021-11-18 JP to be more instructive and general.
54 Parameters
55 ----------
56 w_obs
57 resp_obs
58 zpk_obs
59 zpk_pred
60 w_values
61 xlim
62 title
64 Returns
65 -------
67 """
69 fig = plt.figure(figsize=(14, 4))
71 if poles is not None or zeros is not None:
72 gs = GridSpec(2, 3)
73 ax_amp = fig.add_subplot(gs[0, :2])
74 ax_phs = fig.add_subplot(gs[1, :2], sharex=ax_amp)
75 ax_pz = fig.add_subplot(gs[:, 2], aspect="equal")
76 else:
77 gs = GridSpec(2, 2)
78 ax_amp = fig.add_subplot(gs[0, :2])
79 ax_phs = fig.add_subplot(gs[1, :2], sharex=ax_amp)
81 if not isinstance(complex_response, list):
82 complex_response = [complex_response]
84 if not isinstance(label, list):
85 label = [label]
87 assert len(label) == len(complex_response)
89 if x_units.lower() == "period":
90 frequencies = 1.0 / frequencies
92 amp_min = []
93 amp_max = []
94 phase_min = []
95 phase_max = []
96 lines = []
97 for name, cr in zip(label, complex_response):
98 response_amplitude = np.absolute(cr)
99 if is_flat_amplitude(cr):
100 response_amplitude[:] = response_amplitude[0]
102 if unwrap:
103 response_phase = np.rad2deg(np.unwrap(np.angle(cr, deg=False)))
104 else:
105 response_phase = np.angle(cr, deg=True)
107 amp_min.append(response_amplitude.min())
108 amp_max.append(response_amplitude.max())
109 phase_min.append(response_phase.min())
110 phase_max.append(response_phase.max())
112 # plot amplitude
113 l1 = ax_amp.plot(
114 frequencies,
115 response_amplitude,
116 linewidth=1.5,
117 linestyle="-",
118 label=name,
119 )
120 lines.append(l1)
122 # plot phase
123 ax_phs.plot(
124 frequencies,
125 response_phase,
126 linewidth=1.5,
127 linestyle="-",
128 )
130 # plot pass band
131 if pass_band is not None:
132 if x_units.lower() == "period":
133 pass_band = 1.0 / pass_band
135 ax_amp.fill_between(
136 pass_band,
137 [10e-20, 10e-20],
138 [10e10, 10e10],
139 color=(0.7, 0.7, 0.7),
140 alpha=0.7,
141 zorder=1,
142 )
144 ax_phs.fill_between(
145 pass_band,
146 [-100000, -100000],
147 [100000, 100000],
148 color=(0.7, 0.7, 0.7),
149 alpha=0.7,
150 zorder=1,
151 )
153 if poles is not None:
154 ax_pz.scatter(
155 np.real(poles),
156 np.imag(poles),
157 s=75,
158 marker="x",
159 ec="tab:blue",
160 fc="tab:blue",
161 label="Poles",
162 )
163 if zeros is not None:
164 ax_pz.scatter(
165 np.real(zeros),
166 np.imag(zeros),
167 s=75,
168 marker="o",
169 ec="tab:blue",
170 fc="w",
171 label="Zeros",
172 )
174 if normalization_frequency is not None:
175 if x_units.lower() == "period":
176 try:
177 normalization_frequency = 1.0 / normalization_frequency
178 except ZeroDivisionError:
179 pass
180 ax_amp.plot(
181 [normalization_frequency, normalization_frequency],
182 [10e-20, 10e10],
183 color="k",
184 lw=2,
185 )
186 ax_phs.plot(
187 [normalization_frequency, normalization_frequency],
188 [-10000, 10000],
189 color="k",
190 lw=2,
191 )
193 if xlim is not None:
194 ax_amp.set_xlim(xlim)
195 else:
196 if x_units.lower() == "period":
197 ax_amp.set_xlim([frequencies.max(), frequencies.min()])
198 else:
199 ax_amp.set_xlim([frequencies.min(), frequencies.max()])
201 ax_amp.set_xscale("log")
202 ax_amp.set_yscale("log")
203 ax_amp.set_ylabel("Amplitude Response")
204 ax_amp.grid()
205 ax_amp.set_ylim(
206 [
207 10 ** (np.floor(np.log10(min(amp_min))) - 1),
208 10 ** (np.ceil(np.log10(max(amp_max))) + 1),
209 ]
210 )
212 if not unwrap:
213 ax_phs.set_ylim([-200.0, 200.0])
215 else:
216 ax_phs.set_ylim([min(phase_min) - 10, max(phase_max) + 10])
218 ax_phs.set_xscale("log")
219 ax_phs.set_ylabel("Phase Response")
221 if x_units.lower() == "period":
222 x_label = "Period (s)"
223 elif x_units.lower() == "frequency":
224 x_label = "Frequency (Hz)"
226 ax_phs.set_xlabel(x_label)
227 ax_phs.grid()
229 if poles is not None or zeros is not None:
230 ax_pz.set_xlabel("Re(z)")
231 ax_pz.set_ylabel("Im(z)")
232 max_lim = max(
233 [
234 abs(ax_pz.get_ylim()[0]),
235 abs(ax_pz.get_ylim()[1]),
236 abs(ax_pz.get_xlim()[0]),
237 abs(ax_pz.get_xlim()[0]),
238 ]
239 )
240 ax_pz.set_ylim([-1.25 * max_lim, 1.25 * max_lim])
241 ax_pz.set_xlim([-1.25 * max_lim, 1.25 * max_lim])
242 ax_pz.grid()
243 ax_pz.legend()
245 if len(label) > 1:
246 fig.legend(
247 lines,
248 label,
249 ncol=len(label),
250 loc="upper center",
251 borderaxespad=-0.5,
252 )
253 else:
254 fig.suptitle(title)
255 plt.tight_layout()
257 plt.show()