Changeset 799 for SHX/trunk


Ignore:
Timestamp:
10/25/12 17:21:42 (8 years ago)
Author:
marcus
Message:

...

Location:
SHX/trunk/SeismicHandler
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/modules/stations.py

    r796 r799  
    1313from SeismicHandler.basics import Singleton 
    1414 
     15FUTURE_DATE = UTCDateTime(3e10) 
     16 
    1517 
    1618Base = declarative_base() 
    17  
    18  
    1919class ChannelMeta(Base): 
    2020    """ 
     
    6767 
    6868    def __init__(self, network, station, location, stream, component, \ 
    69                        arraycode, latitude, longitude, elevation, offsetx, \ 
    70                        offsety, gain, poles, zeros, start, end=None, \ 
    71                        description="", depth="0."): 
     69                       latitude, longitude, elevation, gain, poles, zeros, \ 
     70                       start, end=None, description="", depth="0.", \ 
     71                       arraycode='', offsetx='', offsety=''): 
    7272 
    7373        self.channel = ".".join([network, station, location, \ 
     
    9090        self.ondatems = start.microsecond 
    9191 
     92        if end == FUTURE_DATE: 
     93            end = None 
    9294        if end: 
    9395            self.offdate = (end - end.microsecond/1e6).datetime 
     
    9698    def __repr__(self): 
    9799        if self.offdate: 
    98             return "<Channel:%s %s-%s>" % ( 
     100            return "<Channel:%s %s to %s>" % ( 
    99101                self.channel, 
    100102                self.ondate.strftime("%Y-%m-%d"), 
     
    102104            ) 
    103105        else: 
    104             return "<Channel:%s %s-NOW>" % ( 
     106            return "<Channel:%s %s to NOW>" % ( 
    105107                self.channel, 
    106108                self.ondate.strftime("%Y-%m-%d") 
     
    291293 
    292294            # None == open end 
    293             if m.end is not None and m.end < date: 
     295            if m.end is not None and m.end <= date: 
    294296                continue 
    295297 
     
    302304        raise KeyError("no meta data of '%s' found for '%s'" % (code, date)) 
    303305 
    304     def add(self, station, update=False): 
     306    def add(self, station, replace=False): 
    305307        """ 
    306308        Add station to database. 
     
    313315        session = self.dbsessions[self.dbreadwrite] 
    314316 
    315         # check for update 
     317        # check for conflict 
    316318        try: 
    317319            update = self[station.channel, station.start] 
     320            if update == station: 
     321                return 
     322 
     323            if not replace: 
     324                raise ValueError("Concurrent data present! % " % station) 
     325 
     326            session.delete(update) 
    318327        except: 
    319328            pass 
  • SHX/trunk/SeismicHandler/tools/import_inventory.py

    r796 r799  
    1414import re 
    1515import sys 
     16import codecs 
    1617import fnmatch 
    1718import optparse 
    1819from SeismicHandler.core import Stations 
    1920from SeismicHandler.config import Settings 
    20 from SeismicHandler.modules.stations import resolveStations, ChannelMeta 
     21from SeismicHandler.modules.stations import resolveStations, ChannelMeta, FUTURE_DATE 
    2122from SeismicHandler.modules.filters import FilterFile 
    2223from obspy.core import UTCDateTime 
    2324from obspy.sh.core import toUTCDateTime 
     25 
     26 
     27class Timespans(object): 
     28    def __init__(self): 
     29        self.timespans = {} 
     30 
     31    def add(self, start, stop, value): 
     32        self.timespans[(start, stop)] = value 
     33 
     34    def get_value(self, date): 
     35        lowest = (UTCDateTime(), UTCDateTime()) 
     36        for _i in self.timespans: 
     37            if _i[0] < lowest[0]: 
     38                lowest = _i 
     39 
     40            if _i[0] <= date and _i[1] >= date: 
     41                return self.timespans[_i] 
     42 
     43        # by convention: use lowest start if date equals epoch time zero 
     44        if date == UTCDateTime(0): 
     45            return self.timespans[lowest] 
     46 
     47        raise IndexError("Datum not found!") 
     48 
     49    def spans(self): 
     50        return self.timespans.keys() 
     51 
     52    def __len__(self): 
     53        return len(self.timespans) 
     54 
     55    def __add__(self, other): 
     56        """ 
     57        Intersect time spans. Returns only list of timespans NOT Timespans 
     58        instance!! 
     59        """ 
     60        chunks = [] 
     61        all = self.timespans.keys() + other.timespans.keys() 
     62        # make unique 
     63        loop = [] 
     64        while len(all): 
     65            t = all.pop() 
     66            if t not in all: 
     67                loop.append(t) 
     68 
     69        if len(loop) == 1: 
     70            return loop 
     71 
     72        loop.sort() 
     73        for _i, t in enumerate(loop): 
     74            for n in loop[_i+1:]: 
     75                if n == t: 
     76                    chunks.append(n) 
     77                    continue 
     78                if t[0] == n[0]: 
     79                    chunks.append([t[0], min(t[1], n[1])]) 
     80                    break 
     81                if t[1] >= n[0]: 
     82                    chunks.append([t[0], n[0]]) 
     83                    break 
     84        # last one is always appended 
     85        if list(t) not in chunks: 
     86            chunks.append(t) 
     87 
     88        return chunks 
     89 
    2490 
    2591def parse_statinf(name, verbose=False): 
     
    38104 
    39105    data = {} 
    40     with open(name) as f: 
     106    with codecs.open(name, "rb", "ascii") as f: 
    41107        content = f.readlines() 
    42108 
     
    151217 
    152218        if ffile not in filters: 
    153             fname = "TF_VEL_%s.FLF" % ffile 
    154             if verbose: 
    155                 print "reading %s..." % fname 
    156             pz_vel = FilterFile(os.path.join(fdir, fname)) 
    157             fname = "TF_DSP_%s.FLF" % ffile 
    158             if verbose: 
    159                 print "reading %s..." % fname 
    160             pz_dsp = FilterFile(os.path.join(fdir, fname)) 
    161             filters[ffile] = [pz_vel, pz_dsp] 
     219            fnames = ["TF_VEL_S+%s.FLF" % ffile, "TF_VEL_%s.FLF" % ffile] 
     220            for fname in fnames: 
     221                if verbose: 
     222                    print "reading %s..." % fname 
     223                try: 
     224                    _f = FilterFile(os.path.join(fdir, fname)) 
     225                    pz_vel = _f.parse() 
     226                    break 
     227                except IOError: 
     228                    continue 
     229            filters[ffile] = pz_vel 
    162230 
    163231        if start == "...": 
     
    167235 
    168236        if stop == "...": 
    169             stop = None 
     237            # dummy value 
     238            stop = FUTURE_DATE 
    170239        else: 
    171240            start = toUTCDateTime(stop) 
     
    201270 
    202271        if l[2] == "...": 
    203             stop = None 
     272            stop = FUTURE_DATE 
    204273        else: 
    205274            stop = toUTCDateTime(l[2]) 
     
    212281    return gains 
    213282 
    214 def merge(stainf, filt, sens): 
     283def merge(stainf, filt, sens, verbose): 
    215284    """ 
    216285    Merge information together. 
     
    220289    for s in stainf: 
    221290        # filters seems not to contain wildcards, so they reflect full streams, 
    222         # e.g. GRA1-BH-E 
     291        # e.g. GRA1-BH-E - if not: your're messed up ^o^ 
    223292        filters = {} 
    224         for f in filt: 
    225             if not f[0].startswith(s['station']): 
     293        for strm in filt: 
     294            if not strm[0].startswith(s['station']): 
    226295                continue 
     296            d =  list(strm[1:]) 
     297            d.append(filt[strm]) 
    227298            try: 
    228                 filters[f[0]].append(f[1:]) 
     299                filters[strm[0]].append(d) 
    229300            except KeyError: 
    230                 filters[f[0]] = [f[1:]] 
     301                filters[strm[0]] = [d] 
    231302 
    232303        sensitivities = {} 
     
    243314 
    244315        # loop through components 
    245         for f in filters: 
     316        for strm in filters: 
     317            fts = Timespans() 
     318            for _i in filters[strm]: 
     319                fts.add(*_i) 
     320 
    246321            # use fnmatch to respect wildcards 
    247322            for g in sensitivities: 
    248                 if not fnmatch.fnmatch(f, g): 
     323                if not fnmatch.fnmatch(strm, g): 
    249324                    continue 
    250325 
    251                 gains = sensitivities[g] 
     326                gts = Timespans() 
     327                for _i in sensitivities[g]: 
     328                    gts.add(*_i) 
    252329                break 
    253330 
    254331            # loop through time spans 
    255             for t in filters[f]: 
    256                 if len(gains) == 1: 
    257                     # if there's only one, take it 
    258                     # (assuming that input data is right) 
    259                     gain = gains[0] 
    260                 else: 
    261                     try: 
    262                         gain = gains[gains.index(list(t))] 
    263                     except ValueError: 
    264                         print "no gain found for %s in %s" % (f, t) 
    265                         import pdb; pdb.set_trace() 
    266                         continue 
     332            if s['station'].startswith("CLZ"): 
     333                import pdb; pdb.set_trace() 
     334 
     335            timespans = fts + gts 
     336            for t in timespans: 
     337                try: 
     338                    gain = gts.get_value(t[0]) 
     339                except IndexError: 
     340                    print "no gain defined for %s in %s" % (strm, t) 
     341                    import pdb; pdb.set_trace() 
     342                    continue 
    267343 
    268344                # use maximum start date 
    269345                # happens only if startdate is set at sensitivities and not set 
    270346                # at filters 
    271                 date = max(list(t), gain[:2]) 
    272                 _, stream, comp = f.split('-') 
     347                _, stream, comp = strm.split('-') 
     348 
     349                ffilter = fts.get_value(t[0]) 
    273350 
    274351                s["stream"] = stream 
    275352                s["component"] = comp 
    276                 s["gain"] = gain[-1] 
    277                 s["zeros"] = "zeros" 
    278                 s["poles"] = "poles" 
    279                 s["start"], s["end"] = date 
     353                s["gain"] = gain 
     354                s["zeros"] = str(ffilter[0].input_coeff) 
     355                s["poles"] = str(ffilter[0].output_coeff) 
     356                s["start"], s["end"] = t 
    280357 
    281358                try: 
    282359                    meta = ChannelMeta(**s) 
    283                 except: 
     360                except Exception as E: 
     361                    print E 
    284362                    import pdb; pdb.set_trace() 
    285363                    continue 
    286364 
    287                 Stations.add(meta, update=True) 
     365                try: 
     366                    if verbose: 
     367                        msg = "Save information for %s from %s..." 
     368                        print msg % (strm, t[0]) 
     369 
     370                    Stations.add(meta) 
     371                except Exception as E: 
     372                    print E 
     373                    import pdb; pdb.set_trace() 
     374                    print meta 
    288375 
    289376def main(): 
     
    310397    sens = parse_gain(args[2], options.verbose) 
    311398 
    312     merge(statinf, filters, sens) 
     399    merge(statinf, filters, sens, options.verbose) 
    313400 
    314401if __name__ == "__main__": 
Note: See TracChangeset for help on using the changeset viewer.