Changeset 796 for SHX/trunk


Ignore:
Timestamp:
10/24/12 18:04:13 (8 years ago)
Author:
marcus
Message:
  • more work on inventory import
Location:
SHX/trunk/SeismicHandler
Files:
3 edited

Legend:

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

    r788 r796  
     1#! /usr/bin/env python 
    12# -*- coding: utf-8 -*- 
    23 
    3 #    This file is part of Seismic Handler eXtended (SHX). 
    4 # 
    5 #    SHX is free software: you can redistribute it and/or modify 
    6 #    it under the terms of the GNU Lesser General Public License as published 
    7 #    by the Free Software Foundation, either version 3 of the License, or 
    8 #    (at your option) any later version. 
    9 # 
    10 #    SHX is distributed in the hope that it will be useful, 
    11 #    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    12 #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
    13 #    GNU Lesser General Public License for more details. 
    14 # 
    15 #    You should have received a copy of the GNU Lesser General Public License 
    16 #    along with SHX.  If not, see <http://www.gnu.org/licenses/>. 
     4#    This file is part of Seismic Handler eXtended (SHX). For terms of use and 
     5#    license information please see license.txt and visit 
     6#    http://www.seismic-handler.org/portal/wiki/Shx/LicenseTerms 
    177 
    188""" 
    199Filter file access. Mainly taken from obshandler work meeting in 2010 at 
    2010Fuerstenfeldbruck observatory. 
     11 
     12Thanks to Moritz Beyreuther, Lion Krischer and Robert Barsch. 
    2113""" 
    22 # Thanks to Moritz Beyreuther, Lion Krischer and Robert Barsch. 
    2314 
    2415import os 
  • SHX/trunk/SeismicHandler/modules/stations.py

    r788 r796  
    6767 
    6868    def __init__(self, network, station, location, stream, component, \ 
    69                        arraycode, latitude, longitude, elevation, depth, \ 
    70                        offsetx, offsety, gain, poles, zeros, start, end=None, \ 
    71                        description=""): 
     69                       arraycode, latitude, longitude, elevation, offsetx, \ 
     70                       offsety, gain, poles, zeros, start, end=None, \ 
     71                       description="", depth="0."): 
    7272 
    7373        self.channel = ".".join([network, station, location, \ 
     
    302302        raise KeyError("no meta data of '%s' found for '%s'" % (code, date)) 
    303303 
    304     def add(self, station): 
     304    def add(self, station, update=False): 
    305305        """ 
    306306        Add station to database. 
  • SHX/trunk/SeismicHandler/tools/import_inventory.py

    r788 r796  
    1414import re 
    1515import sys 
     16import fnmatch 
    1617import optparse 
    1718from SeismicHandler.core import Stations 
    1819from SeismicHandler.config import Settings 
    19 from SeismicHandler.modules.stations import resolveStations 
     20from SeismicHandler.modules.stations import resolveStations, ChannelMeta 
    2021from SeismicHandler.modules.filters import FilterFile 
     22from obspy.core import UTCDateTime 
     23from obspy.sh.core import toUTCDateTime 
    2124 
    2225def parse_statinf(name, verbose=False): 
     
    2427    Parse basic station information. 
    2528    """ 
    26     f = open(name) 
    27  
    2829    # mapping to db columns 
    2930    map = { 
     
    3334        "yrel": "offsety", 
    3435        "name": "description", 
     36        "array": "arraycode", 
    3537    } 
    3638 
    3739    data = {} 
    38     for l in f: 
     40    with open(name) as f: 
     41        content = f.readlines() 
     42 
     43    for l in content: 
    3944        if l.startswith("#") or l.startswith("!"): 
    4045            continue 
    4146 
    42         l = l.strip() 
     47        l = l.strip().upper() 
    4348        l = ": ".join(l.split(":")) 
    4449 
     
    7782 
    7883            try: 
    79                 data[station][map[key]] = value 
     84                data[station][map[key.lower()]] = value 
    8085            except KeyError: 
    81                 data[station][key] = value 
    82  
    83     f.close() 
     86                data[station][key.lower()] = value 
    8487 
    8588    final = [] 
     
    116119        fdir = os.getenv(fdir[1:]) 
    117120 
    118     f = open(name) 
    119  
    120121    filters = {} 
    121122    assigns = {} 
    122     for l in f: 
    123         l = l.strip() 
     123    with open(name) as f: 
     124        content = f.readlines() 
     125 
     126    for l in content: 
     127        l = l.strip().upper() 
    124128 
    125129        if not l: 
     
    129133            continue 
    130134 
    131         import pdb; pdb.set_trace() 
    132         canal, ffile = l.split() 
    133  
    134         sta, str, comp = canal.split('-') 
     135        l = l.split() 
     136        start = "..." 
     137        stop = "..." 
     138        if len(l) == 2: 
     139            canal, ffile = l 
     140        elif len(l) == 4: 
     141            canal, ffile, start, stop = l 
     142        elif len(l) == 5: 
     143            # we only use FLF filters, so FORCEFFT has no effect 
     144            canal, ffile, start, stop, _ = l 
     145        else: 
     146            print >> sys.stderr, ( 
     147                "Invalid filter lookup information:\n" 
     148                " ".join(l) 
     149            ) 
     150            quit(1) 
    135151 
    136152        if ffile not in filters: 
    137             pz = FilterFile(os.path.join(fdir, "TF_VEL_%s.FLF" % ffile)) 
    138             filters[ffile] = pz 
    139  
    140     f.close() 
    141  
    142  
     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] 
     162 
     163        if start == "...": 
     164            start = UTCDateTime(0) 
     165        else: 
     166            start = toUTCDateTime(start) 
     167 
     168        if stop == "...": 
     169            stop = None 
     170        else: 
     171            start = toUTCDateTime(stop) 
     172 
     173        # assign to stream 
     174        if (canal, start, stop) in assigns.keys(): 
     175            print "Duplicate information for stream %s!" % canal 
     176 
     177        assigns[(canal, start, stop)] = filters[ffile] 
     178 
     179    return assigns 
    143180 
    144181def parse_gain(name, verbose=False): 
     
    146183    Parse sensitivity information. 
    147184    """ 
    148     pass 
    149  
    150 def merge(stainf, filt, sens, fdir): 
     185    gains = {} 
     186    with open(name) as f: 
     187        content = f.readlines() 
     188 
     189    for l in content: 
     190        l = l.strip().upper() 
     191 
     192        if not l or l.startswith("!"): 
     193            continue 
     194 
     195        l = l.split() 
     196 
     197        if l[1] == "...": 
     198            start = UTCDateTime(0) 
     199        else: 
     200            start = toUTCDateTime(l[1]) 
     201 
     202        if l[2] == "...": 
     203            stop = None 
     204        else: 
     205            stop = toUTCDateTime(l[2]) 
     206 
     207        gains[(l[0], start, stop)] = l[3] 
     208 
     209    if verbose: 
     210        print "%u gain values read" % len(gains) 
     211 
     212    return gains 
     213 
     214def merge(stainf, filt, sens): 
    151215    """ 
    152216    Merge information together. 
    153217    """ 
    154     pass 
     218     
     219    # loop through stations 
     220    for s in stainf: 
     221        # filters seems not to contain wildcards, so they reflect full streams, 
     222        # e.g. GRA1-BH-E 
     223        filters = {} 
     224        for f in filt: 
     225            if not f[0].startswith(s['station']): 
     226                continue 
     227            try: 
     228                filters[f[0]].append(f[1:]) 
     229            except KeyError: 
     230                filters[f[0]] = [f[1:]] 
     231 
     232        sensitivities = {} 
     233        for g in sens: 
     234            if not g[0].startswith(s['station']): 
     235                continue 
     236 
     237            d = list(g[1:]) 
     238            d.append(sens[g]) 
     239            try: 
     240                sensitivities[g[0]].append(d) 
     241            except KeyError: 
     242                sensitivities[g[0]] = [d] 
     243 
     244        # loop through components 
     245        for f in filters: 
     246            # use fnmatch to respect wildcards 
     247            for g in sensitivities: 
     248                if not fnmatch.fnmatch(f, g): 
     249                    continue 
     250 
     251                gains = sensitivities[g] 
     252                break 
     253 
     254            # 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 
     267 
     268                # use maximum start date 
     269                # happens only if startdate is set at sensitivities and not set 
     270                # at filters 
     271                date = max(list(t), gain[:2]) 
     272                _, stream, comp = f.split('-') 
     273 
     274                s["stream"] = stream 
     275                s["component"] = comp 
     276                s["gain"] = gain[-1] 
     277                s["zeros"] = "zeros" 
     278                s["poles"] = "poles" 
     279                s["start"], s["end"] = date 
     280 
     281                try: 
     282                    meta = ChannelMeta(**s) 
     283                except: 
     284                    import pdb; pdb.set_trace() 
     285                    continue 
     286 
     287                Stations.add(meta, update=True) 
    155288 
    156289def main(): 
Note: See TracChangeset for help on using the changeset viewer.