-
Notifications
You must be signed in to change notification settings - Fork 0
/
create_ppsd_npz_v2.py
executable file
·198 lines (167 loc) · 6.31 KB
/
create_ppsd_npz_v2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#!/usr/bin/env python3
import datetime
import os
from glob import glob
import traceback
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.patheffects as pe
import numpy as np
from numpy import loadtxt
import pandas as pd
import multiprocessing
from multiprocessing import cpu_count
from obspy import UTCDateTime, read
from obspy.clients.fdsn import Client
from obspy.signal import PPSD
from timeit import default_timer as timer
from itertools import cycle
def create_npz(network,station,channel,start,end,path):
# create a list of days
datelist = pd.date_range(start.datetime, end.datetime, freq="D")
# print(datelist)
# add GeoNet FDSN clients
c = Client(data_provider)
nrt_client = Client("http://service-nrt.geonet.org.nz")
# format if mseed is saved
nslc = "{}.{}.{}".format(network, station, channel)
# make sure that wildcard characters are not in nslc
nslc = nslc.replace("*", "").replace("?", "")
for day in datelist:
datestr = day.strftime("%Y-%m-%d")
fn = "{}_{}.mseed".format(datestr, nslc)
# print(fn)
if day != datelist[-1] and os.path.isfile(path+station+'/'+fn):
continue
else:
try:
location = '10'
args = (network, station, location, channel,
UTCDateTime(day)-1801, UTCDateTime(day)+86400+1801,)
# print('Trying location code 10 for '+station)
# print(day)
# print(station)
st = c.get_waveforms(*args, attach_response=True)
st.merge(fill_value='interpolate')
# print(st)
ppsd(st=st,network=network,station=station,location=location,channel=channel,day=day,path=path,datestr=datestr)
except Exception as e:
try:
location = '11'
args = (network, station, location, channel,
UTCDateTime(day)-1801, UTCDateTime(day)+86400+1801,)
# print('Trying location code 11 for '+station)
st = c.get_waveforms(*args, attach_response=True)
st.merge(fill_value='interpolate')
# print(st)
ppsd(st=st,network=network,station=station,location=location,channel=channel,day=day,path=path,datestr=datestr)
except Exception as e:
# print(e)
print(day)
print('No broadband FDSN data available for '+station)
# # st.write(fn)
# try:
# resp = c.get_stations(UTCDateTime(day), network=network, station=station, location='10',
# channel=channel, level="response")
# except Exception as e:
# try:
# resp = c.get_stations(UTCDateTime(day), network=network, station=station, location='11',
# channel=channel, level="response")
# except Exception as e:
# print('No inventory can be found for '+station)
# # print(resp)
# for id in list(set([tr.id for tr in st])):
# fn_out = "{}_{}.npz".format(datestr, id)
# if os.path.isfile(path+station+'/'+fn_out):
# print("%s completed previously."%fn_out)
# continue
# st = st.select(id=id)
# st.attach_response(resp)
# # ppsd = PPSD(st[0].stats, metadata=resp,
# # ppsd_length=1800, overlap=0.5,
# # period_smoothing_width_octaves=0.025,
# # period_step_octaves=0.0125,
# # period_limits=(0.008, 50),
# # db_bins=(-200, 20, 0.25))
# ppsd = PPSD(st[0].stats, metadata=resp)
# ppsd.add(st)
# ppsd.save_npz(path+station+'/'+fn_out[:-4])
# print(st)
# del st, ppsd
def ppsd(st,network,station,location,channel,day,path,datestr):
c = Client(data_provider)
resp = c.get_stations(UTCDateTime(day), network=network, station=station, location=location,
channel=channel, level="response")
print(resp)
for id in list(set([tr.id for tr in st])):
fn_out = "{}_{}.npz".format(datestr, id)
#print(fn_out)
if os.path.isfile(path+station+'/'+fn_out):
print("%s completed previously."%fn_out)
continue
st = st.select(id=id)
# print(st)
st.attach_response(resp)
# (after McNamara and Buland 2004, obspy default)
ppsd = PPSD(st[0].stats, metadata=resp,
db_bins=(-200, 20, 1.0), ppsd_length=3600.0,
overlap=0.5, special_handling=None,
period_smoothing_width_octaves=1.0,
period_step_octaves=0.125)
# (testing after ThomasLecocq/SeismoRMS wherby the PSD is set to nervous rather than smooth "default PQLX" and half hour 50% overlaps,
# the binning period of 0.0125 octaves and averaging over 0.125 octaves at each central frequency is not required here and creates issues
# at longer periods)
# ppsd = PPSD(st[0].stats, metadata=resp,
# ppsd_length=1800, overlap=0.5,
# period_smoothing_width_octaves=0.025,
# period_step_octaves=0.0125,
# period_limits=(0.008, 50),
# db_bins=(-200, 20, 0.25))
ppsd.add(st)
ppsd.save_npz(path+station+'/'+fn_out[:-4])
print(st)
del st, ppsd
def main():
pool = multiprocessing.Pool(processes=processes, maxtasksperchild=1)
pool = multiprocessing.get_context('spawn').Pool(processes=processes)
args = []
for station in stations:
args += [(network, station, channel, start, end, path)]
print(args)
pool.starmap(create_npz, args)
pool.close()
pool.join()
processes = cpu_count()
network = "NZ"
start = UTCDateTime("2020-01-01")
end = UTCDateTime("2020-12-31")
stations = loadtxt("NI_NN_list.txt", dtype=str, unpack=False)
channel = "HHZ"
path = "/home/conradb/git/NZ_noise_models/NI_stations/"
data_provider = "GEONET"
# some test sites and times to try sites with varying LC's and data gaps
# start = UTCDateTime('2021-03-06T00:00:00')
# end = UTCDateTime('2021-03-08T00:00:00')
# stations = ['BFZ', 'RIZ', 'FOZ', 'URZ']
# main program
if __name__ == "__main__":
start_code = timer()
multiprocessing.set_start_method("spawn")
main()
end_code = timer()
print(end_code - start_code)
# non-multiprocessing option, if running comment out the main() function and the __name__ == "__main__" block above
# for station in stations:
# try:
# create_npz(network=network,station=station,location='10',channel=channel,start=start,end=end,path=path)
# except:
# try:
# create_npz(network=network,station=station,location='11',channel=channel,start=start,end=end,path=path)
# except:
# print('no FDSN data available for site ' + station)
# end_code = timer()
# print(end_code - start_code)
# for reference using all NI stations (ex RIZ/CTZ) over 3 days processing times are:
# with multiprocessing = 246.40s
# without multiprocessing = 420.81s