Source Code for Module ivs.catalogs.gcpd

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  Interface to Geneva's Genaral Catalogue of Photometric Data 
  4  """ 
  5  import urllib 
  6  import logging 
  7  import os 
  8  import ConfigParser 
  9   
 10  import numpy as np 
 11   
 12  from ivs.aux import loggers 
 13  from ivs.aux import numpy_ext 
 14  from ivs.catalogs import sesame 
 15  from ivs.catalogs import vizier 
 16  from ivs.sed import filters 
 17  from ivs.units import conversions 
 18   
 19  #-- the photometric systems in GCPD are represented by a number 
 20  systems = {'JOHNSON':1, 
 21             'STROMGREN':4, 
 22             'ARGUE':10, 
 23             'GENEVA':13, 
 24             'KRON':19, 
 25             'VILNIUS':21, 
 26             'WOOD':22, 
 27             'WBVR':78, 
 28             'STRAIZYS':80, 
 29             'DDO':12} 
 30   
 31  logger = logging.getLogger("CAT.VIZIER") 
 32  logger.addHandler(loggers.NullHandler()) 
 33   
 34  basedir = os.path.dirname(os.path.abspath(__file__)) 
 35   
 36  #-- read in catalog information 
 37  cat_info = ConfigParser.ConfigParser() 
 38  cat_info.optionxform = str # make sure the options are case sensitive 
 39  cat_info.readfp(open(os.path.join(basedir,'gcpd_cats_phot.cfg'))) 
 40   
 41   
 42  #{ Basic interfaces 
 43   
 44 -def search(name,**kwargs): 
 45      """ 
 46      Search and retrieve information from the GCPD catalog. 
 47       
 48      @param name: name of photometric system 
 49      @type name: string 
 50      """ 
 51      base_url = _get_URI(name,**kwargs) 
 52       
 53      #-- the data is listed in two lines: one with the header, one with 
 54      #   the values 
 55      webpage = urllib.urlopen(base_url) 
 56      entries,values = None,None 
 57      log_message = '0' 
 58      start = -1 
 59      for line in webpage: 
 60          line = line.replace('\n','') 
 61          #-- change this marker to be much more general: now it only works for 
 62          #   geneva photometry 
 63          if r'<HR>' in line: start = 0 
 64          if start==1 and r'<B>' in line and not 'No values' in line: 
 65              entries = line.split('<B>')[1].split('</B>')[0].split('\t') 
 66              #-- sometimes there's two columns with the name N: make them unique 
 67              entries = [(entry=='N' and 'N%d'%(i) or entry) for i,entry in enumerate(entries)] 
 68              log_message = '1' 
 69              continue 
 70          if entries is not None: 
 71              values = line.split('\t') 
 72              break 
 73          start += 1 
 74      logger.info('Querying GCPD for %s (%s)'%(name,log_message)) 
 75      webpage.close() 
 76       
 77      #-- assume that all entries are floats, and are values are magnitudes 
 78      if entries is not None: 
 79          if len(values) < len(entries): 
 80              values += ['nan']*(len(entries)-len(values)) 
 81          #-- in some columns, there are '*' or 'STD' signs 
 82          values = [value.replace('*','').replace('STD','').replace('/','') for value in values] 
 83          #-- some columns have no values 
 84          values = [(value and value or 'nan') for value in values] 
 85          #-- some columns are upper limits (or whatever ':' means) 
 86          values = [(':' in value and 'nan' or value) for value in values] 
 87          dtypes = np.dtype([(i,'f8') for i in entries]) 
 88          units = {} 
 89          for entry in entries: 
 90              units[entry] = 'mag' 
 91          cols = [np.cast[dtypes[i]](np.array([values[i]])) for i in range(len(values))] 
 92          results = np.rec.array(cols, dtype=dtypes) 
 93      else: 
 94          results,units,comms = None,None,None 
 95       
 96      return results, units, None 
 97       
 98       
 99 -def gcpd2phot(source,results,units,master=None,e_flag='e_',q_flag='q_',extra_fields=None): 
100      """ 
101      Convert/combine GCPD record arrays to measurement record arrays. 
102       
103      Every line in the combined array represents a measurement in a certain band. 
104       
105      The standard columns are: 
106       
107          1. C{meas}: containing the photometric measurement 
108          2. C{e_meas}: the error on the photometric measurement 
109          3. C{flag}: an optional quality flag 
110          4. C{unit}: the unit of the measurement 
111          5. C{photband}: the photometric passband (FILTER.BAND) 
112          6. C{source}: name of the source catalog 
113       
114      If you give a C{master}, the information will be added to a previous 
115      record array. If not, a new master will be created. 
116       
117      Colors will be expanded, derived from the other columns and added to the 
118      master. 
119       
120      The result is a record array with each row a measurement. 
121       
122      Extra fields are not available for the GCPD, they will be filled in with 
123      nans. 
124       
125      @param source: name of the VizieR source 
126      @type source: str 
127      @param results: results from VizieR C{search} 
128      @type results: record array 
129      @param units: header of Vizier catalog with key name the column name and 
130      key value the units of that column 
131      @type units: dict 
132      @param master: master record array to add information to 
133      @type master: record array 
134      @param e_flag: flag denoting the error on a column 
135      @type e_flag: str 
136      @param q_flag: flag denoting the quality of a measurement 
137      @type q_flag: str 
138      @param extra_fields: any extra columns you want to add information from 
139      @type extra_fields: list of str 
140      @return: array with each row a measurement 
141      @rtype: record array 
142      """ 
143      if cat_info.has_option(source,'e_flag'): 
144          e_flag = cat_info.get(source,'e_flag') 
145       
146      #-- basic dtypes 
147      dtypes = [('meas','f8'),('e_meas','f8'),('flag','a20'), 
148                    ('unit','a30'),('photband','a30'),('source','a50')] 
149       
150      #-- extra can be added, but only if a master is already given!! The reason 
151      #   is that thre GCPD actually does not contain any extra information, so 
152      #   we will never be able to add it and will not know what dtype the extra 
153      #   columns should be    
154      #-- extra can be added: 
155      names = list(results.dtype.names) 
156      if extra_fields is not None: 
157          for e_dtype in extra_fields: 
158              dtypes.append((e_dtype,'f8')) 
159        
160      #-- create empty master if not given 
161      newmaster = False 
162      if master is None or len(master)==0: 
163          master = np.rec.array([tuple([('f' in dt[1]) and np.nan or 'none' for dt in dtypes])],dtype=dtypes) 
164          newmaster = True 
165       
166      #-- add fluxes and magnitudes to the record array 
167      cols_added = 0 
168      for key in cat_info.options(source): 
169          if key[:2] =='e_' or key=='bibcode': 
170              continue 
171          photband = cat_info.get(source,key) 
172          #-- contains measurement, error, quality, units, photometric band and 
173          #   source 
174          cols = [results[key][:1], 
175                  (e_flag+key in cat_info.options(source) and results[cat_info.get(source,e_flag+key)][:1] or np.ones(len(results[:1]))*np.nan), 
176                  (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 
177                  np.array(len(results[:1])*[units[key]],str), 
178                  np.array(len(results[:1])*[photband],str), 
179                  np.array(len(results[:1])*['GCPD'],str)] 
180          #-- add any extra fields if desired. 
181          if extra_fields is not None: 
182              for e_dtype in extra_fields: 
183                  cols += [e_dtype in results.dtype.names and results[e_dtype][:1] or np.ones(len(results[:1]))*np.nan] 
184          #-- add to the master 
185          rows = [] 
186          for i in range(len(cols[0])): 
187              rows.append(tuple([col[i] for col in cols])) 
188          master = np.core.records.fromrecords(master.tolist()+rows,dtype=dtypes) 
189          cols_added += 1 
190       
191      #-- fix colours: we have to run through it two times to be sure to have 
192      #   all the colours 
193      N = len(master)-cols_added 
194      master_ = vizier._breakup_colours(master[N:]) 
195      for i in range(5): 
196          master_ = vizier._breakup_colours(master_) 
197      #-- combine and return 
198      master = np.core.records.fromrecords(master.tolist()[:N]+master_.tolist(),dtype=dtypes) 
199       
200      #-- skip first line from building  
201      if newmaster: master = master[1:] 
202      return master 
203   
204 -def get_photometry(ID=None,extra_fields=[],**kwargs): 
205      """ 
206      Download all available photometry from a star to a record array. 
207       
208      Extra fields will not be useful probably. 
209       
210      For extra kwargs, see L{_get_URI} and L{gcpd2phot} 
211      """ 
212      to_units = kwargs.pop('to_units','erg/s/cm2/AA') 
213      master_ = kwargs.get('master',None) 
214      master = None 
215      #-- retrieve all measurements 
216      for source in cat_info.sections(): 
217          results,units,comms = search(source,ID=ID,**kwargs) 
218          if results is not None: 
219              master = gcpd2phot(source,results,units,master,extra_fields=extra_fields) 
220       
221      #-- convert the measurement to a common unit. 
222      if to_units and master is not None: 
223          #-- prepare columns to extend to basic master 
224          dtypes = [('cwave','f8'),('cmeas','f8'),('e_cmeas','f8'),('cunit','a50')] 
225          cols = [[],[],[],[]] 
226          #-- forget about 'nan' errors for the moment 
227          no_errors = np.isnan(master['e_meas']) 
228          master['e_meas'][no_errors] = 0. 
229          #-- extend basic master 
230          try: 
231              zp = filters.get_info(master['photband']) 
232          except: 
233              print master['photband'] 
234              raise 
235          for i in range(len(master)): 
236              to_units_ = to_units+'' 
237              try: 
238                  value,e_value = conversions.convert(master['unit'][i],to_units,master['meas'][i],master['e_meas'][i],photband=master['photband'][i]) 
239              except ValueError: # calibrations not available 
240                  # if it is a magnitude color, try converting it to a flux ratio 
241                  if 'mag' in master['unit'][i]: 
242                      try: 
243                          value,e_value = conversions.convert('mag_color','flux_ratio',master['meas'][i],master['e_meas'][i],photband=master['photband'][i]) 
244                          to_units_ = 'flux_ratio' 
245                      except ValueError: 
246                          value,e_value = np.nan,np.nan 
247                  # else, we are powerless... 
248                  else: 
249                      value,e_value = np.nan,np.nan 
250              try: 
251                  eff_wave = filters.eff_wave(master['photband'][i]) 
252              except IOError: 
253                  eff_wave = np.nan 
254              cols[0].append(eff_wave) 
255              cols[1].append(value) 
256              cols[2].append(e_value) 
257              cols[3].append(to_units_) 
258          master = numpy_ext.recarr_addcols(master,cols,dtypes) 
259          #-- reset errors 
260          master['e_meas'][no_errors] = np.nan 
261          master['e_cmeas'][no_errors] = np.nan 
262       
263      #-- if a master is given as a keyword, and data is found in this module, 
264      #   append the two 
265      if master_ is not None and master is not None: 
266          master = numpy_ext.recarr_addrows(master_,master.tolist()) 
267      elif master is None: 
268          master = master_ 
269      #-- and return the results 
270      return master 
271   
272  #} 
273  #{ Internal helper functions 
274   
275 -def _get_URI(name='GENEVA',ID=None,**kwargs): 
276      """ 
277      Build GCPD URI from available options. 
278       
279      kwargs are to catch unused arguments. 
280       
281      @param name: photometric system name (E.g. JOHNSON, STROMGREN, GENEVA...) 
282      @type name: string (automatically uppercased) 
283      @keyword ID: star name (should be resolvable by SIMBAD) 
284      @type ID: string 
285      @return:  
286      """ 
287      #-- GCPD is poor at recognising aliases: therefore we try different 
288      #   identifiers retrieved from Sesame that GCPD understands 
289      recognized_alias = ['HD','BD',"CD"] 
290       
291      try: 
292          aliases = sesame.search(ID)['alias'] 
293          for alias in aliases: 
294              if alias[:2] in recognized_alias: 
295                  ID = alias[:2]+' '+alias[2:] 
296                  break 
297          else: 
298              logger.error('Star %s has no aliases recognised by GCPD: query will not return results'%(ID)) 
299      except KeyError: 
300          logger.error('Unknown star %s: GCPD query will not return results'%(ID)) 
301       
302       
303      base_url = 'http://obswww.unige.ch/gcpd/cgi-bin/photoSys.cgi?phot=%02d&type=original&refer=with&mode=starno&ident=%s'%(systems[name],urllib.quote(ID)) 
304      logger.debug(base_url) 
305      return base_url 
306       
307  #} 
308   
309  if __name__=="__main__": 
310      results,units,comms = search('GENEVA',ID='HD180642') 
311      master = gcpd2phot('GENEVA',results,units) 
312      print master 
313      print "" 
314      master = get_photometry(ID='vega') 
315      print master 
316