management.py 10.5 KB
Newer Older
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
from ..global_config import get_setting
from ..core.graph import StuettNode

import os 

import dask
import logging
import obspy
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obsplus import obspy_to_array

 #TODO: revisit the following packages
import numpy as np
import pandas as pd
import datetime as dt 

import warnings
class DataSource(StuettNode):
    def __init__(self):
        pass

    def configure(self,requests=None):
        """ Default configure for DataSource nodes
            Same as configure from StuettNode but adds is_source flag

        Arguments:
            request {list} -- List of requests

        Returns:
            dict -- Original request or merged requests 
        """
        requests = super().configure(requests)
        requests['requires_request'] = True

        return requests

class Freezer(StuettNode):
    def __init__(self):
        pass

matthmey's avatar
matthmey committed
42
43
    def configure(self,requests):
        """ 
44
45
46
47
48

        Arguments:
            request {list} -- List of requests

        Returns:
matthmey's avatar
matthmey committed
49
            dict -- Original, updated or merged request(s) 
50
51
        """
        requests = super().configure(requests)
matthmey's avatar
matthmey committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

        #TODO: check if data is available for requested time period

        #TODO: check how we need to update the boundaries such that we get data that fits and that is available
        #TODO: with only one start/end_time it might be inefficient for the case where [unavailable,available,unavailable] since we need to load everything
        #      one option could be to duplicate the graph by returning multiple requests...

        #TODO: make a distinction between requested start_time and freeze_output_start_time

        #TODO: add node specific hash to freeze_output_start_time (there might be multiple in the graph) <- probably not necessart becaue we receive a copy of the request which is unique to this node
        #TODO: maybe the configuration method must add (and delete) the node name in the request?


        # we always require a request to crop out the right time period
        requests['requires_request'] = True
67
68
69

        return requests

matthmey's avatar
matthmey committed
70
71
72
73
74
75
76
77
78
79
80
81
82
    @dask.delayed
    def __call__(self,x=None, requests=None):

        #TODO: check requests if none
        if requests is None:
            #bypass
            #TODO: maybe a warning
            return x

        #TODO: check request start_time and load the data which is available, store the data which is not available


        #TODO: crop
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


class SeismicSource(DataSource):
    def __init__(self,config={},use_arclink=False):
        """ Seismic data source to get data from permasense
            The user can predefine the source's settings or provide them in a request
            Predefined setting should never be updated (to stay thread safe), but will be ignored
            if request contains settings
        Arguments:
            DataSource {[type]} -- [description]
            config {dict} -- configuration for seismic sources
        """
        self.config = config
        self.use_arclink = use_arclink

        if 'source' not in self.config:
            self.config['source'] = None

        # TODO: make fallback to permasense_vault
        if use_arclink:
            arclink_user = get_setting('arclink_user')
            arclink_password = get_setting('arclink_password')
            self.fdsn_client = Client(base_url='http://arclink.ethz.ch',user=arclink_user,password=arclink_password)

    @dask.delayed
    def __call__(self,request=None):
        # logging.info('Loading seismic with fdsn')
        print('loading with seismic')
        
        # TODO: check if request contains settings
        config = self.config.copy()
        if request is not None:
            config.update(request)

        print('config',config)

        if self.use_arclink:
            st = self.fdsn_client.get_waveforms(network='4D', station=config['station'], location='A',
                                        channel=config['channel'], starttime=UTCDateTime(config['start_time']), endtime=UTCDateTime(config['end_time']), attach_response=True)
        else:   #20180914 is last full day available in permasense_vault
            st = self.get_obspy_stream(config['start_time'],config['end_time'],config['station'],config['channel'])
            print(st)
        x = obspy_to_array(st)


matthmey's avatar
matthmey committed
128
        # print(x)
129
130
        starttime = obspy.UTCDateTime(x.starttime.values).datetime
        starttime =  pd.to_datetime(starttime,utc=True)
matthmey's avatar
matthmey committed
131
132
        # print('starttime',starttime)
        # print(x['time'].values)
133
        timedeltas = pd.to_timedelta(x['time'].values,unit='seconds')
matthmey's avatar
matthmey committed
134
135
        # print('timedeltas',timedeltas)
        # print(type(starttime),type(timedeltas))
136
        xt = starttime + timedeltas
matthmey's avatar
matthmey committed
137
        # print(type(xt),xt)
138
        x['time'] = pd.to_datetime(xt,utc=True) 
matthmey's avatar
matthmey committed
139
        # print('xtime',x['time'])
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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288

        return x

    def get_obspy_stream(self,
                        tbeg,
                        tend,
                        station,
                        channels,
                        pad=False,
                        verbose=False,
                        fill=0,
                        fill_sampling_rate=1000):
        """    
        Loads the microseismic data for the given timeframe into a miniseed file.

        Arguments:
            tbeg {datetime} -- start timestamp of the desired obspy stream
            tend {datetime} -- end timestamp of the desired obspy stream
        
        Keyword Arguments:
            pad {bool} -- If padding is true, the data will be zero padded if the data is not consistent
            fill {} -- If numpy.nan or fill value: error in the seismic stream will be filled with the value. If None no fill will be used
            verbose {bool} -- If info should be printed

        Returns:
            obspy stream -- obspy stream with up to three channels
                            the stream's channels will be sorted alphabetically
        """

        if not isinstance(channels,list):
            channels = [channels]
        datadir = os.path.join(get_setting('permasense_vault_dir'), '') + 'geophones/binaries/PS/%s/'%station

        if not os.path.isdir(datadir):
            # TODO: should this be an error or only a warning. In a period execution this could stop the whole script
            raise IOError('Cannot find the path {}. Please provide a correct path to the permasense_vault directory'.format(datadir))

        # We will get the full hours seismic data and trim it to the desired length afterwards
        tbeg_hours = pd.to_datetime(tbeg).replace(minute=0,second=0,microsecond=0)
        timerange = pd.date_range(start=tbeg_hours,end=tend, freq='H')

        non_existing_files_ts = []      # keep track of nonexisting files

        # drawback of memmap files is that we need to calculate the size beforehand
        stream = obspy.Stream()

        idx = 0
        ## loop through all hours
        for i in range(len(timerange)):
            # start = time.time()
            h = timerange[i]

            st_list  = obspy.Stream()

            datayear = timerange[i].strftime('%Y/')
            sn = 'MHDL' if station == 'MH36' else 'MHDT'    # TODO: do not hardcode it
            filenames = {}
            for channel in channels:
                filenames[channel] = datadir + datayear + '%s.D/PS.%s.A.%s.D.'%(channel,sn,channel)  + timerange[i].strftime('%Y%m%d_%H%M%S') + '.miniseed'
                print(filenames[channel])
                if (not os.path.isfile(filenames[channel])):
                    non_existing_files_ts.append(timerange[i])

                    warnings.warn(RuntimeWarning('Warning file not found for {} {}'.format(timerange[i],filenames[channel])))
                    # if verbose:
                    #     print('Warning files not found for',timerange[i])
                    #     print(filenames[channel])

                    if fill is not None:
                        # create empty stream of fill value
                        arr = np.ones((int(np.ceil(fill_sampling_rate*dt.timedelta(hours=1).total_seconds())),)) * fill
                        st = obspy.Stream([obspy.Trace(arr)])
                    else:
                        continue
                else:
                    st = obspy.read(filenames[channel])


                st_list += st

            stream_h    = st_list.merge(method=0, fill_value=fill)
            start_time  = tbeg if i == 0                else h
            end_time    = tend if i == len(timerange)-1 else h+dt.timedelta(hours=1)
            segment_h   = stream_h.trim(obspy.UTCDateTime(start_time),obspy.UTCDateTime(end_time),pad=pad,fill_value=fill)
            stream += segment_h

        stream = stream.merge(method=0, fill_value=fill)

        stream.sort(keys=['channel'])       # TODO: change this so that the order of the input channels list is maintained

        return stream


class LabelSource(DataSource):
    def __init__(self,config={}):
        """ Label source retrieves the available labels for a given request.
            The labels are stored in a specific format

        Arguments:
            DataSource {[type]} -- [description]
            config {dict} -- configuration for labels
        """
        self.config = config

        if 'source' not in self.config:
            self.config['source'] = None

        # TODO: make fallback to 
        arclink_user = get_setting('arclink_user')
        arclink_password = get_setting('arclink_password')
        self.fdsn_client = Client(base_url='http://arclink.ethz.ch',user=arclink_user,password=arclink_password)

    def __call__(self,request):
        if request is None:
            raise RuntimeError('No request provided, cannot provide data')
        
        # TODO: check if request contains settings
        config = self.config.copy()
        config.update(request)

        st = self.fdsn_client.get_waveforms(network='4D', station=config['station'], location='A',
                                    channel=config['channel'], starttime=UTCDateTime(config['start_time']), endtime=UTCDateTime(config['end_time']), attach_response=True)

        x = obspy_to_array(st)

        return x



class PytorchDataset(DataSource):
    def __init__(self,config={}):
        """ Creates a pytorch like dataset from a data source and a labelsource.
            
        Arguments:
            DataSource {[type]} -- [description]
            config {dict} -- configuration for labels
        """
        self.config = config

        if 'source' not in self.config:
            self.config['source'] = None


    def __call__(self,request):
        if request is None:
            raise RuntimeError('No request provided, cannot provide data')
        

        return x