Wiki source code of FDSN Guide
Show last authors
author | version | line-number | content |
---|---|---|---|
1 | {{box cssClass="floatinginfobox" title="**Contents**"}} | ||
2 | {{toc/}} | ||
3 | {{/box}} | ||
4 | |||
5 | = [[How to Install ObsPy>>url:https://github.com/obspy/obspy/wiki#installation]] = | ||
6 | |||
7 | = [[Seed-Vault>>https://github.com/AuScope/seed-vault]] = | ||
8 | |||
9 | = Connecting to an FDSN Server = | ||
10 | |||
11 | == How to connect to AusPass with & without authenticated access == | ||
12 | |||
13 | |||
14 | {{code language="python"}} | ||
15 | import obspy | ||
16 | from obspy.clients.fdsn import Client | ||
17 | |||
18 | # Initialize FDSN client for AusPass | ||
19 | |||
20 | # For open access data, no username or password is required. | ||
21 | client = Client('AUSPASS') | ||
22 | |||
23 | # To access restricted data, supply your username and password | ||
24 | # Replace 'Z1' and '12345' with your actual credentials | ||
25 | client = Client('AUSPASS', user='Z1', password='12345') | ||
26 | {{/code}} | ||
27 | |||
28 | = Station Metadata = | ||
29 | |||
30 | Information such as site locations, sensor and data logger types, response information, etc are in the station metadata. This can be accessed directly(link) or via the obspy get_stations (link) tool. | ||
31 | |||
32 | |||
33 | == How to download event, station, instrument response == | ||
34 | |||
35 | |||
36 | {{code language="python"}} | ||
37 | import obspy | ||
38 | from obspy.clients.fdsn import Client | ||
39 | |||
40 | # Use AusPass client for station, waveform, and earthquake information | ||
41 | client = Client("AUSPASS") | ||
42 | |||
43 | |||
44 | # Download station information for AUMTC station in S1 network at the response level | ||
45 | inv = client.get_stations(network="S1", station="AUMTC", location="*", | ||
46 | channel="*", level="response") | ||
47 | |||
48 | # Inventory metadata is stored in a Inventory > Network > Station > Channel hierarchy | ||
49 | |||
50 | print(inv) #inventory level | ||
51 | |||
52 | print(inv[0]) # network level (the first network in the inventory) | ||
53 | |||
54 | print(inv[0][0]) # station level (the first station of the first network in the inventory) | ||
55 | |||
56 | print(inv[0][0][0]) # channel level (the first channel of the first station of the first network in the inventoy) | ||
57 | |||
58 | # you can also select items directly | ||
59 | |||
60 | print(inv.select(station='AUMTC',channel='HHZ')[0][0][0]) | ||
61 | |||
62 | # instrument response is attached to a channel object | ||
63 | |||
64 | response = inv.select(station='AUMTC',channel='HHZ')[0][0][0].response | ||
65 | {{/code}} | ||
66 | |||
67 | |||
68 | |||
69 | |||
70 | = Waveform Data = | ||
71 | |||
72 | Waveform data (e.g. the actual seismic data) can be accessed directly (link) or via obspy's get_waveforms (link) tool. It can also be accessed via various tools such as seed-vault, pyweed, etc (add links). | ||
73 | |||
74 | == Downloading and Storing data == | ||
75 | |||
76 | === How to download waveform data === | ||
77 | |||
78 | {{code language="python"}} | ||
79 | from obspy import UTCDateTime | ||
80 | from obspy.clients.fdsn import Client | ||
81 | |||
82 | # Initialize the FDSN client (you can also specify other data centers) | ||
83 | client = Client("AUSPASS") | ||
84 | |||
85 | # Event information | ||
86 | network = "S1" | ||
87 | station = "AUGRF" | ||
88 | starttime = UTCDateTime("2021-09-21T23:15:53") # The time of the earthquake | ||
89 | endtime = starttime + 360 # One hour of data after the earthquake | ||
90 | |||
91 | # Download the MiniSEED data | ||
92 | st = client.get_waveforms(network=network, station=station, location="*", channel="BHZ", | ||
93 | starttime=starttime, endtime=endtime) | ||
94 | # Save the stream to a MiniSEED file | ||
95 | st.write("Woodspoint_2021.mseed", format="MSEED") | ||
96 | print("Downloaded and saved the MiniSEED file.") | ||
97 | {{/code}} | ||
98 | |||
99 | === How to download a LOT of waveform data === | ||
100 | |||
101 | {{code language="python"}} | ||
102 | from obspy import UTCDateTime | ||
103 | from obspy.clients.fdsn import Client | ||
104 | import datetime | ||
105 | |||
106 | # Initialize client and set parameters | ||
107 | client = Client("IRIS") | ||
108 | start_date = UTCDateTime("2023-07-01") | ||
109 | end_date = UTCDateTime("2023-07-02") | ||
110 | network = "S1" | ||
111 | station = "AUMTS" | ||
112 | channel = "BHZ" | ||
113 | |||
114 | # Loop to download data one day at a time | ||
115 | while start_date <= end_date: | ||
116 | next_date = start_date + datetime.timedelta(days=1) | ||
117 | try: | ||
118 | st = client.get_waveforms(network, station, "*", channel, start_date, next_date) | ||
119 | st.write(f"{start_date.date}.mseed", format="MSEED") | ||
120 | except: | ||
121 | print(f"Failed for {start_date} to {next_date}") | ||
122 | start_date = next_date | ||
123 | {{/code}} | ||
124 | |||
125 | === How to store and archive waveform data in SeisComP Data Structure (SDS) === | ||
126 | |||
127 | {{code language="python"}} | ||
128 | import os | ||
129 | from obspy import UTCDateTime, read | ||
130 | from obspy.clients.fdsn import Client | ||
131 | |||
132 | # Initialize the client | ||
133 | client = Client("AUSPASS") # Replace with the correct client endpoint if different | ||
134 | |||
135 | # Define event time and time window | ||
136 | event_time = UTCDateTime("2021-09-21T23:15:53") | ||
137 | starttime = event_time - 600 # 10 minutes before the event | ||
138 | endtime = event_time + 1800 # 30 minutes after the event | ||
139 | |||
140 | # Download the waveform data | ||
141 | st = client.get_waveforms(network="S1", station="AUMTS", location="*", channel="*", starttime=starttime, endtime=endtime) | ||
142 | |||
143 | # Create SDS structure: ROOT/YEAR/NET/STA/CHAN.TYPE/NET.STA.LOC.CHAN.YEAR.DAY | ||
144 | sds_root = "." # Replace with your desired directory | ||
145 | |||
146 | for tr in st: | ||
147 | net = tr.stats.network | ||
148 | sta = tr.stats.station | ||
149 | loc = tr.stats.location | ||
150 | chan = tr.stats.channel | ||
151 | year = str(tr.stats.starttime.year) | ||
152 | jday = str(tr.stats.starttime.julday).zfill(3) | ||
153 | |||
154 | sds_path = os.path.join(sds_root, year, net, sta, f"{chan}.D", f"{net}.{sta}.{loc}.{chan}.{year}.{jday}") | ||
155 | |||
156 | # Create directories if they don't exist | ||
157 | os.makedirs(os.path.dirname(sds_path), exist_ok=True) | ||
158 | |||
159 | # Save the trace as a MiniSEED file | ||
160 | tr.write(sds_path, format="MSEED") | ||
161 | {{/code}} | ||
162 | |||
163 | == | ||
164 | Common Data Operations == | ||
165 | |||
166 | === How to remove instrument response === | ||
167 | |||
168 | {{code language="python"}} | ||
169 | from obspy import read | ||
170 | from obspy.core.util import AttribDict | ||
171 | |||
172 | # Load the MiniSEED file | ||
173 | st = read("Woodspoint_2021.mseed") | ||
174 | |||
175 | # Download the instrument response | ||
176 | inv = client.get_stations(network=network, station=station, location="*", | ||
177 | channel="*", starttime=starttime, endtime=endtime, | ||
178 | level="response") | ||
179 | |||
180 | # Remove the instrument response | ||
181 | output = 'VEL' # Output unit ('VEL' = velocity (default), 'DISP' = displacement, 'ACC' = acceleration) | ||
182 | |||
183 | for tr in st: | ||
184 | tr.remove_response(inventory=inv, output=output, plot=True) | ||
185 | |||
186 | # Save the corrected MiniSEED file | ||
187 | st.write("Woodspoint_2021_corrected.mseed", format="MSEED") | ||
188 | {{/code}} | ||
189 | |||
190 | === How to apply a bandpass filter === | ||
191 | |||
192 | {{code language="python"}} | ||
193 | from obspy import read | ||
194 | |||
195 | # Load the MiniSEED file | ||
196 | st = read("Woodspoint_2021.mseed") | ||
197 | |||
198 | # Define the frequency band | ||
199 | freq_min = 0.1 # Minimum frequency in Hz | ||
200 | freq_max = 1.0 # Maximum frequency in Hz | ||
201 | |||
202 | # Apply the bandpass filter | ||
203 | for tr in st: | ||
204 | tr.filter(type='bandpass', freqmin=freq_min, freqmax=freq_max) | ||
205 | |||
206 | # Save the filtered MiniSEED file | ||
207 | st.write("Woodspoint_2021_filtered.mseed", format="MSEED") | ||
208 | {{/code}} | ||
209 | |||
210 | === How to slice a waveform === | ||
211 | |||
212 | {{code language="python"}} | ||
213 | from obspy import read, UTCDateTime, Stream # Importing Stream here | ||
214 | |||
215 | # Load the filtered MiniSEED file | ||
216 | st = read("Woodspoint_2021_filtered.mseed") | ||
217 | |||
218 | # Define the time window for slicing | ||
219 | slice_start = UTCDateTime("2021-09-21T23:20:00") | ||
220 | slice_end = slice_start +10 | ||
221 | |||
222 | # Slice the waveform for each Trace in the Stream | ||
223 | sliced_st = Stream() # Now Stream is defined | ||
224 | for tr in st: | ||
225 | sliced_tr = tr.slice(starttime=slice_start, endtime=slice_end) | ||
226 | sliced_st.append(sliced_tr) | ||
227 | |||
228 | # Save the sliced MiniSEED file | ||
229 | sliced_st.write("Woodspoint_2021_filtered_sliced.mseed", format="MSEED") | ||
230 | {{/code}} | ||
231 | |||
232 | === How to save a waveform === | ||
233 | |||
234 | {{code language="python"}} | ||
235 | # Save the sliced file as MiniSEED | ||
236 | sliced_st.write("Woodspoint_2021_filtered_sliced.mseed", format="MSEED") | ||
237 | |||
238 | # Or, save the sliced SAC file | ||
239 | sliced_st.write("Woodspoint_2021_filtered_sliced.sac", format="SAC") | ||
240 | {{/code}} | ||
241 | |||
242 | === How to convert miniSEED to SAC === | ||
243 | |||
244 | {{code language="python"}} | ||
245 | from obspy import read | ||
246 | |||
247 | # Read the MiniSEED file | ||
248 | st = read("Woodspoint_2021.mseed") | ||
249 | |||
250 | # Take the first Trace from the Stream | ||
251 | tr = st[0] | ||
252 | |||
253 | # Save that Trace as a SAC file | ||
254 | tr.write("Woodspoint_2021.sac", format="SAC") | ||
255 | {{/code}} | ||
256 | |||
257 | |||
258 | = Earthquake Data = | ||
259 | |||
260 | Earthquake data can be accessed directly or via ObsPy's get_events code | ||
261 | |||
262 | == How to download an Earthquake Catalog == | ||
263 | |||
264 | {{code language="python"}} | ||
265 | from obspy.clients.fdsn import Client | ||
266 | from obspy import UTCDateTime | ||
267 | |||
268 | # Initialize the AusPass FDSN client | ||
269 | client = Client("AUSPASS") | ||
270 | |||
271 | # Define the time range for the earthquake catalog | ||
272 | start_time = UTCDateTime("2021-08-01") | ||
273 | end_time = UTCDateTime("2022-01-01") # End of year | ||
274 | |||
275 | # Define the geographic region (latitude and longitude for Woodspoint, Victoria, Australia) | ||
276 | latitude = -37.47 | ||
277 | longitude = 146.10 | ||
278 | max_radius = 5 # in degrees | ||
279 | |||
280 | # Download the earthquake catalog | ||
281 | catalog = client.get_events(starttime=start_time, endtime=end_time, | ||
282 | minmagnitude=2, latitude=latitude, longitude=longitude, | ||
283 | maxradius=max_radius) | ||
284 | |||
285 | # Save the catalog to a file (e.g., QuakeML format) | ||
286 | catalog.write("Woodspoint_earthquakes.xml", format="QUAKEML") | ||
287 | {{/code}} | ||
288 | |||
289 | == How to plot (Global) Earthquakes == | ||
290 | |||
291 | {{code language="python"}} | ||
292 | from obspy import UTCDateTime | ||
293 | from obspy.clients.fdsn import Client | ||
294 | |||
295 | # Initialize FDSN client to connect to the IRIS data center | ||
296 | client = Client("IRIS") | ||
297 | |||
298 | # Set the time range for fetching earthquake data | ||
299 | # Start time: January 1, 2023 | ||
300 | # End time: Current time | ||
301 | starttime = UTCDateTime("2023-01-01") | ||
302 | endtime = UTCDateTime() | ||
303 | |||
304 | # Fetch earthquake events with a minimum magnitude of 7 | ||
305 | catalog = client.get_events(starttime=starttime, endtime=endtime, minmagnitude=7) | ||
306 | #client.get_events(). This function returns a Catalog object that contains a list of Event objects. | ||
307 | #Each Event object, in turn, has an Origins attribute that contains the depth information | ||
308 | |||
309 | # Plot the fetched earthquake data using an orthographic projection | ||
310 | catalog.plot(projection="ortho", title="Global Earthquakes with Magnitude >= 7 since 2023") | ||
311 | #catalog.plot(), ObsPy automatically uses the depth information to color the events in the plot | ||
312 | {{/code}} | ||
313 | |||
314 | == How to plot (Local) Earthquakes == | ||
315 | |||
316 | {{code language="python"}} | ||
317 | from obspy import UTCDateTime | ||
318 | from obspy.clients.fdsn import Client | ||
319 | |||
320 | # Initialize FDSN client | ||
321 | client = Client("AUSPASS") | ||
322 | |||
323 | # Define time range | ||
324 | starttime = UTCDateTime("2023-01-01") | ||
325 | endtime = UTCDateTime() | ||
326 | |||
327 | # Latitude and longitude bounds for Australia | ||
328 | minlatitude = -44.0 | ||
329 | maxlatitude = -10.0 | ||
330 | minlongitude = 113.0 | ||
331 | maxlongitude = 154.0 | ||
332 | |||
333 | # Fetch event data for Australia with a minimum magnitude | ||
334 | catalog = client.get_events(starttime=starttime, endtime=endtime, minmagnitude=4, | ||
335 | minlatitude=minlatitude, maxlatitude=maxlatitude, | ||
336 | minlongitude=minlongitude, maxlongitude=maxlongitude) | ||
337 | |||
338 | # Plot the earthquakes | ||
339 | catalog.plot(projection="local", title="Australia Earthquakes", resolution="i") | ||
340 | {{/code}} | ||
341 | |||
342 |