1
2 """
3 Interface to the VizieR website.
4
5 Download or retrieve VizieR catalogs.
6
7 The basic interface C{search} lets you download B{entire catalogs or parts} of
8 them. The structure array containts then one row per target, the columns
9 denoting the different columns of the catalogs. You can also specify two
10 catalogs in C{xmatch}; the second will then be cross-matched against the first.
11
12 The convenience functions (C{get_photometry},...) aim to collect information
13 from different catalogs on B{one target} in one array. Each row represents
14 one measurement from one catalog (typically you'll have many rows both from one
15 catalog but also from different catalogs). The columns then denote the contents
16 of each row (e.g. the magnitude, photometric passband etc).
17
18 Section 1. Download catalogs
19 ============================
20
21 Section 1.1. To a file
22 ----------------------
23
24 Download the entire Van Leeuwen Hipparcos New Reduction catalog to a file. The
25 filename is returned as a check for success.
26
27 >>> filename = search('I/311/hip2',filename='vanleeuwen.tsv')
28
29 Download a 60 arcsec circular area from the catalog around the coordinates
30 ra=237.1, dec=-10.10
31
32 >>> filename = search('I/311/hip2',ra=237.1,dec=-10.10,radius=60.,filename='vanleeuwen.tsv')
33
34 Search for the presence of a target in the catalog. The downloaded file will
35 contain no rows if the target is not in the catalog. If more than one target are
36 in the search radius around the target, there will be more than one row. They
37 are ordered via distance to the target, so it's probably the first one you need.
38
39 >>> filename = search('I/311/hip2',ID='vega',filename='vanleeuwen.tsv')
40 >>> filename = search('I/311/hip2',ID='vega',filename='vanleeuwen.tsv',radius=60.)
41
42 Section 1.2 To a RecordArray
43 ----------------------------
44
45 Instead of downloading to a file and then reading in the file for further
46 analysis, you can download the contents of the file directly to a record array,
47 retrieving the units and comments from the catalog in one go. The equivalent of
48 the third example above becomes
49
50 >>> rec_arr,unit_dict,comment_str = search('I/311/hip2',ID='vega')
51
52 With these record arrays, its very easy to select targets based on criteria.
53 For example, if you want to extract 2MASS targets in a certain area with a
54 negative H-K index, you can do
55
56 >>> data,units,comms = search('II/246/out',ra=237.1,dec=-10.10,radius=600.)
57 >>> selection = (data['Hmag'] - data['Kmag']) < 0
58 >>> data = data[selection]
59
60 You can also read in a data file you've previously downloaded via
61
62 >>> data,units,comms = tsv2recarray('vanleeuwen.tsv')
63
64 Section 1.3 List relevant catalogs
65 ----------------------------------
66
67 To know in which catalogs your target is present, list them all via
68
69 >>> my_cats = list_catalogs('vega')
70
71 Now you could iterate over them and download them to a file or whatever.
72
73 Section 2. Convenience functions
74 ================================
75
76 You can define 'standard' photometric catalogs in the C{vizier_cats.cfg} file.
77 This file is basically a translator for VizieR column headers to photometric
78 passbands (and colors). For examples, see the file itself.
79
80 You can add catalogs on the fly via
81
82 >>> cat_info.add_section('my_new_catalog')
83 >>> cat_info.set('my_new_catalog','Bmag','JOHNSON.B')
84 """
85
86 import numpy as np
87 import urllib
88 import logging
89 import os
90 import itertools
91 from astropy.io import fits as pyfits
92 import tarfile
93 import tempfile
94 import shutil
95 import ConfigParser
96 from scipy.spatial import KDTree
97
98
99 from cc.ivs.io import ascii
100 from cc.ivs.io import fits
101 from cc.ivs.io import http
102 from cc.ivs.units import conversions
103 from cc.ivs.aux import loggers
104 from cc.ivs.aux import numpy_ext
105 from cc.ivs.sed import filters
106
107 logger = logging.getLogger("CAT.VIZIER")
108 logger.addHandler(loggers.NullHandler())
109
110 basedir = os.path.dirname(os.path.abspath(__file__))
111
112
113 cat_info = ConfigParser.ConfigParser()
114 cat_info.optionxform = str
115 cat_info.readfp(open(os.path.join(basedir,'vizier_cats_phot.cfg')))
116
117 cat_info_fund = ConfigParser.ConfigParser()
118 cat_info_fund.optionxform = str
119 cat_info_fund.readfp(open(os.path.join(basedir,'vizier_cats_fund.cfg')))
120
121 mirrors = {'cycle': itertools.cycle(['vizier.u-strasbg.fr',
122 'vizier.nao.ac.jp',
123 'vizier.hia.nrc.ca',
124 'vizier.ast.cam.ac.uk',
125 'vizier.cfa.harvard.edu',
126 'www.ukirt.jach.hawaii.edu',
127 'vizier.inasan.ru',
128 'vizier.iucaa.ernet.in',
129 'data.bao.ac.cn'])}
130 mirrors['current'] = mirrors['cycle'].next()
131
132
133
135 """
136 Cycle through the mirrors of ViZieR.
137 """
138 mirrors['current'] = mirrors['cycle'].next()
139 logger.info("Changed cycle to {}".format(mirrors['current']))
140
141 -def search(name,filetype='tsv',filename=None,**kwargs):
142 """
143 Search and retrieve information from a VizieR catalog.
144
145 Two ways to search for data within a catalog C{name}:
146
147 1. You're looking for info on B{one target}, then give the target's
148 C{ID} or coordinates (C{ra} and C{dec}), and a search C{radius}.
149
150 2. You're looking for information of B{a whole field}, then give the
151 field's coordinates (C{ra} and C{dec}), and C{radius}.
152
153 If you have a list of targets, you need to loop this function.
154
155 If you supply a filename, the results will be saved to that path, and you
156 will get the filename back as received from urllib.URLopener (should be the
157 same as the input name, unless something went wrong).
158
159 If you don't supply a filename, you should leave C{filetype} to the default
160 C{tsv}, and the results will be saved to a temporary
161 file and deleted after the function is finished. The content of the file
162 will be read into a dictionary, as well as the units (two separate
163 dictionaries with the same keys, depending on the column names in the
164 catalog). The entries in the dictionary are of type C{ndarray}, and will
165 be converted to a float-array (no integers, we need to support nans) if
166 possible. If not, the array will consist of strings. The comments are also
167 returned as a list of strings.
168
169 WARNING: when retrieving a FITS file, ViZieR sometimes puts weird formats
170 into the header ('3F10.6E' in the 2MASS catalog), which cannot be read by
171 the C{pyfits} module. These columns are actually multi-dimensional vectors.
172 One option is to download to another format, or to restrict the columns with
173 C{out_all=None}.
174
175 Example usage:
176
177 1. Look for the Geneva V magnitude and coordinates of Vega in the GENEVA
178 catalog of Rufener.
179
180 >>> results,units,comms = search('II/169/main',ID='vega',radius=60.)
181 >>> print "Vega: Vmag = %.3f %s, RA = %.2f %s, DEC = %.2f %s"%(results['Vmag'],units['Vmag'],results['_RAJ2000'],units['_RAJ2000'],results['_DEJ2000'],units['_DEJ2000'])
182 Vega: Vmag = 0.061 mag, RA = 279.24 deg, DEC = 38.77 deg
183
184 2. Search for all targets in the 2MASS catalog in a particular field.
185 Download the results to a FITS file, read the file, and plot the results
186 to the screen.
187
188 >>> #filename = search('II/246/out',ra=100.79,dec=0.70,radius=1000.,filetype='fits',filename='2mass_test',out_all=None)
189
190 Now read in the FITS-file and plot the contents
191
192 >>> #from astropy.io import fits as pyfits
193 >>> #import pylab
194 >>> #ff = pyfits.open('2mass_test.fits')
195 >>> #p = pylab.gcf().canvas.set_window_title('test of <search>')
196 >>> #p = pylab.scatter(ff[1].data.field('_RAJ2000'),ff[1].data.field('_DEJ2000'),c=ff[1].data.field('Jmag'),s=(20-ff[1].data.field('Jmag'))**2,cmap=pylab.cm.hot_r,edgecolors='none')
197 >>> #p = pylab.colorbar()
198 >>> #p = p.set_label('Jmag')
199 >>> #p,q = pylab.xlabel('RA [deg]'),pylab.ylabel('DEC [deg]')
200 >>> #ff.close()
201 >>> #os.remove('2mass_test.fits')
202
203
204 @param name: name of a ViZieR catalog (e.g. 'II/246/out')
205 @type name: str
206 @param filetype: type of the file to write the results to ('tsv' if no file desired)
207 @type filetype: string (one of 'tsv','fits','ascii','csv'...
208 @param filename: name of the file to write the results to (no extension)
209 @type filename: str
210 @return: filename / catalog data columns, units, comments
211 @rtype: str/ record array, dict, list of str
212 """
213
214
215 if filename is not None and '.' in os.path.basename(filename):
216 filetype = os.path.splitext(filename)[1][1:]
217 elif filename is not None:
218 filename = '%s.%s'%(filename,filetype)
219
220
221 base_url = _get_URI(name=name,**kwargs)
222
223
224 url = urllib.URLopener()
225 filen,msg = url.retrieve(base_url,filename=filename)
226
227 if filename is not None:
228 logger.info('Querying ViZieR source %s and downloading to %s'%(name,filen))
229 url.close()
230 return filen
231
232
233 if filetype=='tsv':
234 try:
235 results,units,comms = tsv2recarray(filen)
236
237 except ValueError:
238 raise ValueError, "failed to read %s, perhaps multiple catalogs specified (e.g. III/168 instead of III/168/catalog)"%(name)
239 url.close()
240 logger.info('Querying ViZieR source %s (%d)'%(name,(results is not None and len(results) or 0)))
241 return results,units,comms
242
243
245 """
246 Print and return all catalogs containing information on the star.
247
248 If you give C{filetype} and C{filename}, all information will be downloaded
249 to that file.
250
251 Extra kwargs: see L{_get_URI}.
252
253 @param ID: identification of the star
254 @type ID: str
255 @keyword filetype: type of the output file ('fits','tsv','csv'...)
256 @type filetype: str
257 @keyword filename: name of the output file
258 @type filename: str
259 @return: dictionary with keys the catalog ID from VizieR and entries the
260 names of the catalogs
261 @rtype: dictionary
262 """
263 base_url = _get_URI(ID=ID,filetype='fits',**kwargs)
264
265
266 url = urllib.URLopener()
267 filen,msg = url.retrieve(base_url,filename=filename)
268
269
270
271 if filetype=='fits':
272 mycats = {}
273 ff = pyfits.open(filen)
274 for ext in range(1,len(ff)):
275 name = ff[ext].header['CDS-name']
276 results,units,comms = search(name,ID=ID,**kwargs)
277 for line in comms:
278 if 'Title:' in line:
279 title = line.strip().split(':')[1]
280 break
281 mycats[name] = title
282 logger.info('%25s %s'%(name,title))
283
284 photometry = [col for col in units.keys() if 'mag' in units[col]]
285 rv = [col for col in units.keys() if 'rv' in col.lower()]
286 vsini = [col for col in units.keys() if 'sin' in col.lower()]
287 sptype = [col for col in units.keys() if col.lower()=='sp' or col.lower()=='sptype']
288 fund = [col for col in units.keys() if 'logg' in col.lower() or 'teff' in col.lower()]
289
290 ff.close()
291 url.close()
292 return mycats
293 else:
294 url.close()
295 return filen
296
297
298 -def xmatch(source1,source2,output_file=None,tol=1.,**kwargs):
299 """
300 Crossmatch two vizier catalogs via a fast KDTree.
301
302 The limit for these catalogs is probably somewhere between ~100000 entries,
303 so make sure your catalogs do not contain to many targets. You can always
304 do a subselection via the keyword arguments (e.g. give ra, dec and radius).
305
306 An output tsv file will be written (by default named 'source1__source2',
307 which can be read in via C{tsv2recarray} for further analysis.
308
309 tolerance is in arcseconds.
310
311 Extra keyword arguments are passed to C{search}. Column names of second
312 source will be appended with postfix '_2', to avoid clashes of double-defined
313 column headers.
314 """
315
316 if output_file is None:
317 output_file = "__".join([source1,source2]).replace('/','_').replace('+','')+'.tsv'
318
319
320 cat1,units1,comms1 = search(source1,**kwargs)
321 cat2,units2,comms2 = search(source2,**kwargs)
322
323 logger.info('Start Vizier Xmatch')
324 coords1 = np.array([cat1['_RAJ2000'],cat1['_DEJ2000']]).T
325 coords2 = np.array([cat2['_RAJ2000'],cat2['_DEJ2000']]).T
326
327 logger.info('Building KDTree of shape %d,%d'%coords1.shape)
328 tree = KDTree(coords1)
329
330 logger.info('Querying KDTree with %d entries'%(len(coords2)))
331 distance,order = tree.query(coords2)
332
333 keep = distance<(tol/(60.))
334
335 logger.info('Matched %d points (tol<%.3g arcsec)'%(sum(keep),tol))
336
337
338 cat1 = cat1[order[keep]]
339 cat2 = cat2[keep]
340
341
342
343
344 for i,line in enumerate(comms2):
345 line = line.split('\t')
346 if line[0]=='Column':
347 line[1] = line[1]+'_2'
348 comms2[i] = '\t'.join(line)
349
350 units2_ = {}
351 for key in units2: units2_[key+'_2'] = units2[key]
352 cat2.dtype.names = [name+'_2' for name in cat2.dtype.names]
353
354
355 ff = open(output_file,'w')
356 ff.write('\n#'.join(comms1))
357 ff.write('\n#'.join(comms2))
358
359 names1 = list(cat1.dtype.names)
360 names2 = list(cat2.dtype.names)
361 dtypes = [(name,cat1.dtype[names1.index(name)].str) for name in names1]
362 dtypes += [(name,cat2.dtype[names2.index(name)].str) for name in names2]
363
364 ff.write('\n')
365 for nr,i in enumerate(dtypes):
366 ff.write(str(i[0]))
367 if nr<(len(dtypes)-1): ff.write('\t')
368 ff.write('\n')
369 for nr,i in enumerate(dtypes):
370 if i[0] in units1: ff.write(units1[i[0]])
371 elif i[0] in units2_: ff.write(units2_[i[0]])
372 else:
373 raise ValueError,'this cannot be'
374 if nr<(len(dtypes)-1): ff.write('\t')
375
376 ff.write('\n')
377 ff.write('\t'.join(['---']*len(dtypes)))
378 ff.write('\n')
379
380 for row1,row2 in itertools.izip(cat1,cat2):
381 ff.write('\t'.join([str(x) for x in row1])+'\t')
382 ff.write('\t'.join([str(x) for x in row2])+'\n')
383
384 ff.close()
385
386
387
388
389
390
391
392
393 -def get_IUE_spectra(ID=None,directory=None,unzip=True,cat_info=False,select='low',**kwargs):
394 """
395 Download IUE spectra.
396
397 If you want to download all the spectral files, set C{directory='/home/user/'}
398 or whatever. All the tarfiles will be downloaded to this directory, they
399 will be untarred, science data extracted and all unnecessary files and
400 directories will be deleted. If you don't set a directory, it will default
401 to the CWD.
402
403 If you don't wish to unzip them, set unzip=False
404
405 DEPRECATED: If you don't give a directory, the function will return a list
406 of all extracted spectra (no data files are kept).
407
408 You can retrieve the contents of the vizier catalog via {cat_info=True}. The
409 files will not be downloaded in this case.
410 """
411 if directory is None:
412 direc = os.getcwd()
413 directory = os.getcwd()
414 filename = None
415 else:
416 direc = directory
417 if not os.path.isdir(direc):
418 os.mkdir(direc)
419
420 output = []
421
422 data,units,comments = search('VI/110/inescat/',ID=ID,**kwargs)
423
424 if cat_info:
425 return data,units,comments
426
427 if data is None:
428 return output
429
430 for spectrum in data:
431 download_link = "http://archive.stsci.edu/cgi-bin/iue_retrieve?iue_mark=%s%05d&mission=iue&action=Download_MX"%(spectrum['Camera'].strip(),int(spectrum['Image']))
432 logger.info('IUE spectrum %s/%s: %s'%(spectrum['Camera'],spectrum['Image'],download_link))
433
434
435 if directory is not None:
436 filename = download_link.split('iue_mark=')[1].split('&')[0]
437 filename = os.path.join(direc,filename)
438
439 mytarfile = http.download(download_link,filename=filename)
440 if filename is None:
441 mytarfile,url = mytarfile
442
443 if directory is not None and not unzip:
444 output.append(mytarfile)
445 url.close()
446 continue
447
448
449 if not tarfile.is_tarfile(mytarfile):
450 logger.info("Not a tarfile: %s"%(mytarfile))
451 continue
452 tarf = tarfile.open(mytarfile)
453 files = tarf.getmembers(),tarf.getnames()
454 deldirs = []
455 outfile = None
456 for mem,name in zip(*files):
457
458 myname = name.lower()
459 skip = False
460 if 'lo' in select and not 'lo' in os.path.basename(myname):
461 skip = True
462 if 'hi' in select and not 'hi' in os.path.basename(myname):
463 skip = true
464 if ('lwp' in name or 'swp' in name or 'lwr' in name) and not skip:
465
466 tarf.extract(mem,path=direc)
467
468 outfile = os.path.join(direc,os.path.basename(name))
469 dirname = os.path.dirname(name)
470 shutil.move(os.path.join(direc,name),outfile)
471 logger.info("Extracted %s to %s"%(name,outfile))
472
473 deldirs.append(os.path.join(direc,dirname))
474 else:
475 logger.debug("Did not extract %s (probably not a spectrum)"%(name))
476
477 tarf.close()
478 os.unlink(mytarfile)
479 for dirname in deldirs:
480 if dirname and os.path.isdir(dirname) and not os.listdir(dirname):
481 os.rmdir(dirname)
482 logger.debug("Deleted left over (empty) directory %s"%(dirname))
483 if filename is None: url.close()
484
485
486 if directory is not None and outfile:
487 output.append(outfile)
488 continue
489 if outfile and os.path.isfile(outfile):
490 wavelength,flux,error,header = fits.read_iue(outfile,return_header=True)
491 os.unlink(outfile)
492 output.append([wavelength,flux,error,header])
493 else:
494 logger.info('Unsuccesfull extraction of %s'%(outfile))
495
496 return output
497
499 """
500 Get a spectrum from the UVSST spectrograph onboard TD1.
501
502 From vizier catalog III/39A.
503
504 Also have a look at II/86/suppl.
505 """
506 kwargs.setdefault('out_max',10)
507 kwargs.setdefault('radius',60.)
508 data,units_o,comments = search('III/39A/catalog',**kwargs)
509 if data is None: return None,None
510 fields = sorted([field for field in data.dtype.names if (field[0]=='F' and len(field)==5)])
511 spectrum = np.zeros((len(fields),3))
512 units_ = []
513 for i,field in enumerate(fields):
514 wave,flux = float(field[1:]),data[field][0]
515 if flux==0: continue
516
517 flux = conversions.convert(units_o[field],units,100*flux,wave=(wave,'AA'),unpack=True)
518 e_flux = data['e_'+field][0]/100.*flux
519 spectrum[i][0] = wave
520 spectrum[i][1] = flux
521 spectrum[i][2] = e_flux
522 units_.append(units_o[field])
523 return spectrum[spectrum[:,0]>0]
524
525
526
527
528 -def get_photometry(ID=None,extra_fields=['_r','_RAJ2000','_DEJ2000'],take_mean=False,**kwargs):
529 """
530 Download all available photometry from a star to a record array.
531
532 For extra kwargs, see L{_get_URI} and L{vizier2phot}
533
534 """
535 kwargs['ID'] = ID
536 to_units = kwargs.pop('to_units','erg/s/cm2/AA')
537 sources = kwargs.get('sources',cat_info.sections())
538 master_ = kwargs.get('master',None)
539 master = None
540
541 for source in sources:
542 results,units,comms = search(source,**kwargs)
543 if results is None: continue
544 master = vizier2phot(source,results,units,master,extra_fields=extra_fields,take_mean=take_mean)
545
546 if to_units and master is not None:
547
548 dtypes = [('cwave','f8'),('cmeas','f8'),('e_cmeas','f8'),('cunit','a50')]
549 cols = [[],[],[],[]]
550
551 no_errors = np.isnan(master['e_meas'])
552 master['e_meas'][no_errors] = 0.
553
554 zp = filters.get_info(master['photband'])
555 for i in range(len(master)):
556 to_units_ = to_units+''
557 try:
558 value,e_value = conversions.convert(master['unit'][i],to_units,master['meas'][i],master['e_meas'][i],photband=master['photband'][i])
559 except ValueError:
560
561 if 'mag' in master['unit'][i]:
562 try:
563 value,e_value = conversions.convert('mag_color','flux_ratio',master['meas'][i],master['e_meas'][i],photband=master['photband'][i])
564 to_units_ = 'flux_ratio'
565 except ValueError:
566 value,e_value = np.nan,np.nan
567
568 else:
569 value,e_value = np.nan,np.nan
570 try:
571 eff_wave = filters.eff_wave(master['photband'][i])
572 except IOError:
573 eff_wave = np.nan
574 cols[0].append(eff_wave)
575 cols[1].append(value)
576 cols[2].append(e_value)
577 cols[3].append(to_units_)
578 master = numpy_ext.recarr_addcols(master,cols,dtypes)
579
580 master['e_meas'][no_errors] = np.nan
581 master['e_cmeas'][no_errors] = np.nan
582
583 if master_ is not None and master is not None:
584 master = numpy_ext.recarr_addrows(master_,master.tolist())
585 elif master is None:
586 master = master_
587
588
589 return master
590
591
593 """
594 Perform quality checks on downloaded data.
595
596 This function translates flags in to words, and looks up additional
597 information in selected catalogs.
598 """
599 messages = ['' for i in range(len(master))]
600
601
602
603 logger.info('Checking flags')
604 for i,entry in enumerate(master):
605 if 'USNO' in str(entry['photband']):
606 messages[i] = '; '.join([messages[i],'unreliable zeropoint and transmission profile'])
607 continue
608 if 'COUSINS' in str(entry['photband']):
609 messages[i] = '; '.join([messages[i],'unreliable zeropoint'])
610 continue
611 flag = entry['flag']
612 try:
613 flag = float(flag)
614 except:
615 continue
616 if np.isnan(flag): continue
617 if entry['source'].strip() in ['II/125/main','II/275/fsr']:
618 flag = float(int(flag))
619 if flag==3: messages[i] = '; '.join([messages[i],'high quality'])
620 if flag==2: messages[i] = '; '.join([messages[i],'moderate quality'])
621 if flag==1: messages[i] = '; '.join([messages[i],'upper limit'])
622 if entry['source'].strip() in ['II/298/fis','II/297/irc']:
623 flag = float(int(flag))
624 if flag==3: messages[i] = '; '.join([messages[i],'high quality'])
625 if flag==2: messages[i] = '; '.join([messages[i],'source is confirmed but the flux is not reliable'])
626 if flag==1: messages[i] = '; '.join([messages[i],'the source is not confirmed'])
627 if flag==0: messages[i] = '; '.join([messages[i],'not observed'])
628 if 'DIRBE' in entry['source'].strip():
629 flag = '{0:03d}'.format(int(float(flag)))
630 if flag[0]=='1': messages[i] = '; '.join([messages[i],'IRAS/2MASS companion greater than DIRBE noise level'])
631 if flag[1]=='1': messages[i] = '; '.join([messages[i],'possibly extended emission or highly variable source'])
632 if flag[2]=='1': messages[i] = '; '.join([messages[i],'possibly affected by nearby companion'])
633 if entry['source'].strip() in ['V/114/msx6_main']:
634 if flag==4: messages[i] = '; '.join([messages[i],'excellent'])
635 if flag==3: messages[i] = '; '.join([messages[i],'good'])
636 if flag==2: messages[i] = '; '.join([messages[i],'fair'])
637 if flag==1: messages[i] = '; '.join([messages[i],'on limit'])
638 if flag==0: messages[i] = '; '.join([messages[i],'not detected'])
639
640
641 sources_with_quality_check = ['B/denis/denis','II/311/wise','II/246/out']
642 sources = set(list(master['source'])) & set(sources_with_quality_check)
643
644 denis_image_flags = {'01':'clouds during observation',
645 '02':'electronic Read-out problem',
646 '04':'internal temperature problem',
647 '08':'very bright star',
648 '10':'bright star',
649 '20':'stray light',
650 '40':'unknown problem'}
651 denis_source_flags = {'01':'source might be a dust on mirror',
652 '02':'source is a ghost detection of a bright star',
653 '04':'source is saturated',
654 '08':'source is multiple detect',
655 '10':'reserved'}
656 wise_conf_flag = {'d':'diffraction spike',
657 'p':'contaminated by latent image left by bright star',
658 'h':'halo of nearby bright source',
659 'o':'optical ghost'}
660 wise_var_flag = {'n':'too few measurements to decide if variable',
661 '0':'most likely not variable (0, 0=certainly not-5=probably not)',
662 '1':'most likely not variable (1, 0=certainly not-5=probably not)',
663 '2':'most likely not variable (2, 0=certainly not-5=probably not)',
664 '3':'most likely not variable (3, 0=certainly not-5=probably not)',
665 '4':'most likely not variable (4, 0=certainly not-5=probably not)',
666 '5':'most likely not variable (5, 0=certainly not-5=probably not)',
667 '6':'likely variable (6, 6=likely-7=more likely)',
668 '7':'likely variable (7, 6=likely-7=more likely)',
669 '8':'most likely variable (8, 8=most likely-9=almost certain)',
670 '9':'most likely variable (9, 8=most likely-9=almost certain)'}
671 twomass_qual_flag = {'X':'detection, but no valid brightness estimate (X)',
672 'U':'upper limit (U)',
673 'F':'error estimate not reliable (F)',
674 'E':'poor PSF fit (E)',
675 'A':'high quality (A)',
676 'B':'high quality (B)',
677 'C':'high quality (C)',
678 'D':'high quality (D)'}
679
680
681 indices = np.arange(len(master))
682 logger.info('Checking source catalogs for additional information')
683 for source in sorted(sources):
684 results,units,comms = search(source,ID=ID,**kwargs)
685 if results is None: continue
686
687 if source=='B/denis/denis':
688 for iflag,photband in zip(['Iflg','Jflg','Kflg'],['DENIS.I','DENIS.J','DENIS.KS']):
689 index = indices[(master['source']==source) & (master['photband']==photband)]
690 if not len(index):
691 continue
692 flag = float(results[0][iflag])
693 if np.isnan(flag):
694 messages[index] = '; '.join([messages[index],'high quality'])
695 continue
696 flag = '{0:04d}'.format(int(flag))
697 image_flag = flag[:2]
698 source_flag = flag[2:]
699
700 if image_flag in denis_image_flags:
701 messages[index] = '; '.join([messages[index],denis_image_flags[image_flag]])
702 if source_flag in denis_source_flags:
703 messages[index] = '; '.join([messages[index],denis_source_flags[source_flag]])
704 if source=='II/311/wise':
705 conf = results[0]['ccf']
706 var = results[0]['var']
707 ex = int(results[0]['ex'])
708 for i,photband in enumerate(['WISE.W1','WISE.W2','WISE.W3','WISE.W4']):
709 index = indices[(master['source']==source) & (master['photband']==photband)]
710 if len(index)!=1:
711 logger.warning("Skipping WISE flags, don't know what to do with {}".format(index))
712 continue
713 if conf[i]!=' ' and conf[i] in wise_conf_flag:
714 messages[index] = '; '.join([messages[index],wise_conf_flag[conf[i].lower()]])
715 if var[i]!=' ':
716 messages[index] = '; '.join([messages[index],wise_var_flag[var[i].lower()]])
717 messages[index] = '; '.join([messages[index],(ex==0 and 'point source' or 'extended source')])
718 if source=='II/246/out':
719 flag = results[0]['Qflg'].strip()
720 for i,photband in enumerate(['2MASS.J','2MASS.H','2MASS.KS']):
721 if flag[i] in twomass_qual_flag:
722 index = indices[(master['source']==source) & (master['photband']==photband)]
723 if not len(index):
724 continue
725 messages[index] = '; '.join([messages[index],twomass_qual_flag[flag[i]]])
726
727
728
729 for i in range(len(messages)):
730 if messages[i]:
731 messages[i] = messages[i][2:]
732 else:
733 messages[i] = '-'
734 messages[i] = messages[i].replace(' ','_')
735
736
737 messages = np.rec.fromarrays([messages],names=['comments'])
738 if return_master and not 'comments' in master.dtype.names:
739 return numpy_ext.recarr_join(master,messages)
740 elif return_master:
741 master['comments'] = messages
742 return master
743 else:
744 return messages
745
746
747
748
749
750
751
753 """
754 Read a Vizier tsv (tab-sep) file into a record array.
755
756 @param filename: name of the TSV file
757 @type filename: str
758 @return: catalog data columns, units, comments
759 @rtype: record array, dict, list of str
760 """
761 data,comms = ascii.read2array(filename,dtype=np.str,splitchar='\t',return_comments=True)
762 results = None
763 units = {}
764
765 if len(data)>0:
766
767 data = np.array(data)
768
769
770
771
772 formats = np.zeros_like(data[0])
773 for line in comms:
774 line = line.split('\t')
775 if len(line)<3: continue
776 for i,key in enumerate(data[0]):
777 if key == line[1] and line[0]=='Column':
778 formats[i] = line[2].replace('(','').replace(')','').lower()
779 if formats[i][0].isdigit(): formats[i] = 'a100'
780 elif 'f' in formats[i]: formats[i] = 'f8'
781 elif 'i' in formats[i]: formats[i] = 'f8'
782 elif 'e' in formats[i]: formats[i] = 'f8'
783
784 if formats[i][0]=='a':
785 formats[i] = 'a'+str(int(formats[i][1:])+3)
786
787 dtypes = np.dtype([(i,j) for i,j in zip(data[0],formats)])
788
789 cols = []
790 for i,key in enumerate(data[0]):
791 col = data[3:,i]
792
793
794
795
796 cols.append([(row.isspace() or not row) and np.nan or 3*' '+row for row in col])
797
798
799
800
801 units[key] = data[1,i]
802
803 cols = [np.cast[dtypes[i]](cols[i]) for i in range(len(cols))]
804 results = np.rec.array(cols, dtype=dtypes)
805 return results,units,comms
806
807 -def vizier2phot(source,results,units,master=None,e_flag='e_',q_flag='q_',extra_fields=None,take_mean=False):
808 """
809 Convert/combine VizieR record arrays to measurement record arrays.
810
811 Every line in the combined array represents a measurement in a certain band.
812
813 This is probably only useful if C{results} contains only information on
814 one target (or you have to give 'ID' as an extra field, maybe).
815
816 The standard columns are:
817
818 1. C{meas}: containing the photometric measurement
819 2. C{e_meas}: the error on the photometric measurement
820 3. C{flag}: an optional quality flag
821 4. C{unit}: the unit of the measurement
822 5. C{photband}: the photometric passband (FILTER.BAND)
823 6. C{source}: name of the source catalog
824
825 You can add extra information from the VizieR catalog via the list of keys
826 C{extra_fields}.
827
828 If you give a C{master}, the information will be added to a previous
829 record array. If not, a new master will be created.
830
831 Colors will be expanded, derived from the other columns and added to the
832 master.
833
834 The result is a record array with each row a measurement.
835
836 Example usage:
837
838 First look for all photometry of Vega in all VizieR catalogs:
839
840 >>> from cc.ivs.sed import filters
841 >>> import pylab
842 >>> master = None
843 >>> for source in cat_info.sections():
844 ... results,units,comms = search(source,ID='vega',radius=60.)
845 ... if results is not None:
846 ... master = vizier2phot(source,results,units,master,extra_fields=['_r','_RAJ2000','_DEJ2000'])
847
848 Keep only observations we have an measurement and error of, convert every
849 observation to 'Jy' and keep track of the results to plot.
850
851 >>> master = master[(-np.isnan(master['e_meas'])) & (-np.isnan(master['meas']))]
852 >>> eff_waves = filters.eff_wave(master['photband'])
853 >>> myvalue,e_myvalue = conversions.nconvert(master['unit'],'erg/s/cm2/AA',master['meas'],master['e_meas'],photband=master['photband'])
854 >>> for i in range(len(master)):
855 ... print '%15s %10.3e+/-%10.3e %11s %10.3e %3s %6.2f %6.2f %6.3f %23s'%(master[i]['photband'],master[i]['meas'],master[i]['e_meas'],master[i]['unit'],myvalue[i],'Jy',master[i]['_RAJ2000'],master[i]['_DEJ2000'],master[i]['_r'],master[i]['source'])
856 GENEVA.V 6.100e-02+/- 2.500e-02 mag 3.620e-09 Jy 279.24 38.77 56.000 II/169/main
857 GENEVA.B -8.980e-01+/- 2.500e-02 mag 6.518e-09 Jy 279.24 38.77 56.000 II/169/main
858 GENEVA.U 6.070e-01+/- 2.500e-02 mag 3.223e-09 Jy 279.24 38.77 56.000 II/169/main
859 GENEVA.V1 7.640e-01+/- 2.500e-02 mag 3.782e-09 Jy 279.24 38.77 56.000 II/169/main
860 GENEVA.B1 2.000e-03+/- 2.500e-02 mag 6.584e-09 Jy 279.24 38.77 56.000 II/169/main
861 GENEVA.B2 6.120e-01+/- 2.500e-02 mag 6.208e-09 Jy 279.24 38.77 56.000 II/169/main
862 GENEVA.G 1.270e+00+/- 2.500e-02 mag 3.111e-09 Jy 279.24 38.77 56.000 II/169/main
863 2MASS.J -1.770e-01+/- 2.060e-01 mag 3.591e-10 Jy 279.23 38.78 0.034 II/246/out
864 2MASS.H -2.900e-02+/- 1.460e-01 mag 1.176e-10 Jy 279.23 38.78 0.034 II/246/out
865 2MASS.KS 1.290e-01+/- 1.860e-01 mag 3.799e-11 Jy 279.23 38.78 0.034 II/246/out
866 IRAS.F12 4.160e+01+/- 1.664e+00 Jy 1.024e-13 Jy 279.23 38.78 9.400 II/125/main
867 IRAS.F25 1.100e+01+/- 5.500e-01 Jy 6.195e-15 Jy 279.23 38.78 9.400 II/125/main
868 IRAS.F60 9.510e+00+/- 7.608e-01 Jy 8.420e-16 Jy 279.23 38.78 9.400 II/125/main
869 IRAS.F100 7.760e+00+/- 6.984e-01 Jy 2.349e-16 Jy 279.23 38.78 9.400 II/125/main
870 TD1.1565 5.689e-09+/- 1.700e-11 10mW/m2/nm 5.689e-09 Jy 279.23 38.78 18.510 II/59B/catalog
871 TD1.1965 4.928e-09+/- 1.300e-11 10mW/m2/nm 4.928e-09 Jy 279.23 38.78 18.510 II/59B/catalog
872 TD1.2365 3.700e-09+/- 1.000e-11 10mW/m2/nm 3.700e-09 Jy 279.23 38.78 18.510 II/59B/catalog
873 TD1.2740 3.123e-09+/- 9.000e-12 10mW/m2/nm 3.123e-09 Jy 279.23 38.78 18.510 II/59B/catalog
874 ANS.15N -4.910e-01+/- 1.000e-03 mag 5.707e-09 Jy 279.23 38.78 14.000 II/97/ans
875 ANS.15W -4.410e-01+/- 1.200e-02 mag 5.450e-09 Jy 279.23 38.78 14.000 II/97/ans
876 ANS.18 -4.620e-01+/- 3.000e-03 mag 5.556e-09 Jy 279.23 38.78 14.000 II/97/ans
877 ANS.25 4.600e-02+/- 4.000e-03 mag 3.480e-09 Jy 279.23 38.78 14.000 II/97/ans
878 ANS.33 1.910e-01+/- 3.000e-03 mag 3.045e-09 Jy 279.23 38.78 14.000 II/97/ans
879 HIPPARCOS.HP 8.680e-02+/- 2.100e-03 mag 3.840e-09 Jy 279.23 38.78 3.060 I/239/hip_main
880 MIPS.24 8.900e+03+/- 8.900e+01 mJy 4.628e-15 Jy 279.23 38.78 0.010 J/ApJ/653/675/table1
881 MIPS.70 1.142e+04+/- 2.283e+03 mJy 7.075e-16 Jy 279.23 38.78 0.010 J/ApJ/653/675/table1
882 JOHNSON.B 1.900e-02+/- 1.000e-02 mag 6.216e-09 Jy 279.23 38.78 3.060 I/280B/ascc
883 JOHNSON.V 7.400e-02+/- 2.000e-03 mag 3.428e-09 Jy 279.23 38.78 3.060 I/280B/ascc
884 JOHNSON.V 3.300e-02+/- 1.200e-02 mag 3.560e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans
885 JOHNSON.B-V -1.000e-03+/- 5.000e-03 mag nan Jy 279.23 38.78 0.010 II/168/ubvmeans
886 JOHNSON.U-B -6.000e-03+/- 6.000e-03 mag nan Jy 279.23 38.78 0.010 II/168/ubvmeans
887 JOHNSON.B 3.200e-02+/- 1.300e-02 mag 6.142e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans
888 JOHNSON.U 2.600e-02+/- 1.432e-02 mag 4.086e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans
889 JOHNSON.K 1.300e-01+/- 1.900e-01 mag 3.764e-11 Jy 279.23 38.78 0.010 J/PASP/120/1128/catalog
890 AKARI.N60 6.582e+00+/- 2.090e-01 Jy 4.614e-16 Jy 279.23 38.78 3.400 II/298/fis
891 AKARI.WIDES 6.201e+00+/- 1.650e-01 Jy 2.566e-16 Jy 279.23 38.78 3.400 II/298/fis
892 AKARI.WIDEL 4.047e+00+/- 3.500e-01 Jy 5.658e-17 Jy 279.23 38.78 3.400 II/298/fis
893 AKARI.N160 3.221e+00+/- 2.550e-01 Jy 3.695e-17 Jy 279.23 38.78 3.400 II/298/fis
894 AKARI.S9W 5.670e+01+/- 4.010e-01 Jy 2.169e-13 Jy 279.24 38.78 2.550 II/297/irc
895 AKARI.L18W 1.254e+01+/- 7.770e-02 Jy 1.050e-14 Jy 279.24 38.78 2.550 II/297/irc
896 STROMGREN.B-Y 3.000e-03+/- 3.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog
897 STROMGREN.M1 1.570e-01+/- 3.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog
898 STROMGREN.C1 1.088e+00+/- 4.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog
899 STROMGREN.B 4.300e-02+/- 3.000e-03 mag 5.604e-09 Jy 279.23 38.78 22.000 II/215/catalog
900 DIRBE.F2_2 6.217e+02+/- 9.500e+00 Jy 3.791e-11 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
901 DIRBE.F3_5 2.704e+02+/- 1.400e+01 Jy 6.534e-12 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
902 DIRBE.F4_9 1.504e+02+/- 6.200e+00 Jy 1.895e-12 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
903 DIRBE.F12 2.910e+01+/- 1.570e+01 Jy 6.014e-14 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
904 DIRBE.F25 1.630e+01+/- 3.120e+01 Jy 1.153e-14 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
905 DIRBE.F60 1.220e+01+/- 5.610e+01 Jy 1.195e-15 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
906 DIRBE.F100 -7.000e-01+/- 7.790e+01 Jy -2.238e-17 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
907 DIRBE.F140 2.557e+02+/- 5.223e+03 Jy 3.577e-15 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
908 DIRBE.F240 8.290e+01+/- 2.881e+03 Jy 4.152e-16 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE
909 WISE.W2 1.143e+00+/- 1.900e-02 mag 8.428e-13 Jy 279.24 38.78 3.276 II/311/wise
910 WISE.W3 -6.700e-02+/- 8.000e-03 mag 6.930e-14 Jy 279.24 38.78 3.276 II/311/wise
911 WISE.W4 -1.270e-01+/- 6.000e-03 mag 5.722e-15 Jy 279.24 38.78 3.276 II/311/wise
912
913 Make a quick plot:
914
915 >>> p = pylab.figure()
916 >>> p = pylab.loglog(eff_waves,myvalue,'ko')
917 >>> p = pylab.show()
918
919 @param source: name of the VizieR source
920 @type source: str
921 @param results: results from VizieR C{search}
922 @type results: record array
923 @param units: header of Vizier catalog with key name the column name and
924 key value the units of that column
925 @type units: dict
926 @param master: master record array to add information to
927 @type master: record array
928 @param e_flag: flag denoting the error on a column
929 @type e_flag: str
930 @param q_flag: flag denoting the quality of a measurement
931 @type q_flag: str
932 @param extra_fields: any extra columns you want to add information from
933 @type extra_fields: list of str
934 @return: array with each row a measurement
935 @rtype: record array
936 """
937 if cat_info.has_option(source,'q_flag'):
938 q_flag = cat_info.get(source,'q_flag')
939 if cat_info.has_option(source,'e_flag'):
940 e_flag = cat_info.get(source,'e_flag')
941
942
943 dtypes = [('meas','f8'),('e_meas','f8'),('flag','a20'),
944 ('unit','a30'),('photband','a30'),('source','a50')]
945
946
947 names = list(results.dtype.names)
948
949 if extra_fields is not None:
950 for e_dtype in extra_fields:
951 if e_dtype in names:
952 dtypes.append((e_dtype,results.dtype[names.index(e_dtype)].str))
953 else:
954 dtypes.append((e_dtype,'f8'))
955
956
957 newmaster = False
958 if master is None or len(master)==0:
959 master = np.rec.array([tuple([('f' in dt[1]) and np.nan or 'none' for dt in dtypes])],dtype=dtypes)
960 newmaster = True
961
962
963 cols_added = 0
964 for key in cat_info.options(source):
965 if key in ['e_flag','q_flag','mag','bibcode']:
966 continue
967 photband = cat_info.get(source,key)
968 key = key.replace('_[','[').replace(']_',']')
969
970
971 if not take_mean or len(results)<=1:
972 cols = [results[key][:1],
973 (e_flag+key in results.dtype.names and results[e_flag+key][:1] or np.ones(len(results[:1]))*np.nan),
974 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan),
975 np.array(len(results[:1])*[units[key]],str),
976 np.array(len(results[:1])*[photband],str),
977 np.array(len(results[:1])*[source],str)]
978 else:
979 logger.warning('taking mean: {0} --> {1}+/-{2}'.format(results[key][0],results[key].mean(),results[key].std()))
980 cols = [[results[key].mean()],[results[key].std()],
981 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan),
982 np.array(len(results[:1])*[units[key]],str),
983 np.array(len(results[:1])*[photband],str),
984 np.array(len(results[:1])*[source],str)]
985
986 if e_flag+key in results.dtype.names and units[e_flag+key]=='%':
987 cols[1] = cols[1]/100.*cols[0]
988
989 if extra_fields is not None:
990 for e_dtype in extra_fields:
991 cols += [e_dtype in results.dtype.names and results[e_dtype][:1] or np.ones(len(results[:1]))*np.nan]
992
993 rows = []
994 for i in range(len(cols[0])):
995 rows.append(tuple([col[i] for col in cols]))
996 master = np.core.records.fromrecords(master.tolist()+rows,dtype=dtypes)
997 cols_added += 1
998
999
1000 N = len(master)-cols_added
1001 master_ = _breakup_colours(master[N:])
1002 master_ = _breakup_colours(master_)
1003
1004 master = np.core.records.fromrecords(master.tolist()[:N]+master_.tolist(),dtype=dtypes)
1005
1006
1007 if newmaster: master = master[1:]
1008 return master
1009
1010
1011
1012 -def vizier2fund(source,results,units,master=None,e_flag='e_',q_flag='q_',extra_fields=None):
1013 """
1014 Convert/combine VizieR record arrays to measurement record arrays.
1015
1016 This is probably only useful if C{results} contains only information on
1017 one target (or you have to give 'ID' as an extra field, maybe).
1018
1019 The standard columns are:
1020
1021 1. C{meas}: containing the measurement of a fundamental parameter
1022 2. C{e_meas}: the error on the measurement of a fundamental parameter
1023 3. C{flag}: an optional quality flag
1024 4. C{unit}: the unit of the measurement
1025 5. C{source}: name of the source catalog
1026
1027 If a target appears more than once in a catalog, only the first match will
1028 be added.
1029
1030 The result is a record array with each row a measurement.
1031
1032 Example usage:
1033
1034 >>> master = None
1035 >>> for source in cat_info_fund.sections():
1036 ... results,units,comms = search(source,ID='AzV 79',radius=60.)
1037 ... if results is not None:
1038 ... master = vizier2fund(source,results,units,master,extra_fields=['_r','_RAJ2000','_DEJ2000'])
1039 >>> master = master[-np.isnan(master['meas'])]
1040 >>> for i in range(len(master)):
1041 ... print '%10s %10.3e+/-%10.3e %11s %6.2f %6.2f %6.3f %23s'%(master[i]['name'],master[i]['meas'],master[i]['e_meas'],master[i]['unit'],master[i]['_RAJ2000'],master[i]['_DEJ2000'],master[i]['_r'],master[i]['source'])
1042 Teff 7.304e+03+/- nan K 12.67 -72.83 0.002 B/pastel/pastel
1043 logg 2.000e+00+/- nan [cm/s2] 12.67 -72.83 0.002 B/pastel/pastel
1044 [Fe/H] -8.700e-01+/- nan [Sun] 12.67 -72.83 0.002 B/pastel/pastel
1045
1046 @param source: name of the VizieR source
1047 @type source: str
1048 @param results: results from VizieR C{search}
1049 @type results: record array
1050 @param units: header of Vizier catalog with key name the column name and
1051 key value the units of that column
1052 @type units: dict
1053 @param master: master record array to add information to
1054 @type master: record array
1055 @param e_flag: flag denoting the error on a column
1056 @type e_flag: str
1057 @param q_flag: flag denoting the quality of a measurement
1058 @type q_flag: str
1059 @param extra_fields: any extra columns you want to add information from
1060 @type extra_fields: list of str
1061 @return: array with each row a measurement
1062 @rtype: record array
1063 """
1064 if cat_info_fund.has_option(source,'q_flag'):
1065 q_flag = cat_info_fund.get(source,'q_flag')
1066 if cat_info_fund.has_option(source,'e_flag'):
1067 e_flag = cat_info_fund.get(source,'e_flag')
1068
1069
1070 dtypes = [('meas','f8'),('e_meas','f8'),('q_meas','f8'),('unit','a30'),
1071 ('source','a50'),('name','a50')]
1072
1073
1074 names = list(results.dtype.names)
1075 if extra_fields is not None:
1076 for e_dtype in extra_fields:
1077 dtypes.append((e_dtype,results.dtype[names.index(e_dtype)].str))
1078
1079
1080 newmaster = False
1081 if master is None or len(master)==0:
1082 master = np.rec.array([tuple([('f' in dt[1]) and np.nan or 'none' for dt in dtypes])],dtype=dtypes)
1083 newmaster = True
1084
1085
1086 for key in cat_info_fund.options(source):
1087 if key in ['e_flag','q_flag']:
1088 continue
1089 photband = cat_info_fund.get(source,key)
1090
1091
1092
1093 key = key.replace('_[','[').replace(']_',']')
1094 cols = [results[key][:1],
1095 (e_flag+key in results.dtype.names and results[e_flag+key][:1] or np.ones(len(results[:1]))*np.nan),
1096 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan),
1097 np.array(len(results[:1])*[units[key]],str),
1098 np.array(len(results[:1])*[source],str),
1099 np.array(len(results[:1])*[key],str)]
1100
1101 if extra_fields is not None:
1102 for e_dtype in extra_fields:
1103 cols.append(results[:1][e_dtype])
1104
1105 rows = []
1106 for i in range(len(cols[0])):
1107 rows.append(tuple([col[i] for col in cols]))
1108
1109 master = np.core.records.fromrecords(master.tolist()+rows,dtype=dtypes)
1110
1111
1112 if newmaster: master = master[1:]
1113 return master
1114
1115
1117 """
1118 Retrieve the ADS bibcode of a ViZieR catalog.
1119
1120 @param catalog: name of the catalog (e.g. II/306/sdss8)
1121 @type catalog: str
1122 @return: bibtex code
1123 @rtype: str
1124 """
1125 catalog = catalog.split("/")
1126 N = max(2,len(catalog)-1)
1127 catalog = "/".join(catalog[:N])
1128 base_url = "http://cdsarc.u-strasbg.fr/viz-bin/Cat?{0}".format(catalog)
1129 url = urllib.URLopener()
1130 filen,msg = url.retrieve(base_url)
1131
1132 code = None
1133 ff = open(filen,'r')
1134 for line in ff.readlines():
1135 if 'Keywords' in line: break
1136 if '=<A HREF' in line:
1137 code = line.split('>')[1].split('<')[0]
1138 ff.close()
1139 url.close()
1140 return code
1141
1143 """
1144 Retrieve the bibtex entry of an ADS bibcode.
1145
1146 @param bibcode: bibcode (e.g. C{2011yCat.2306....0A})
1147 @type bibcode: str
1148 @return: bibtex entry
1149 @rtype: str
1150 """
1151 base_url = "http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode={0}&data_type=BIBTEX&db_key=AST&nocookieset=1".format(bibcode)
1152 url = urllib.URLopener()
1153 filen,msg = url.retrieve(base_url)
1154
1155 bibtex = []
1156 ff = open(filen,'r')
1157 for line in ff.readlines():
1158 if (not bibtex and '@' in line) or bibtex:
1159 bibtex.append(line)
1160 ff.close()
1161 url.close()
1162 return "".join(bibtex)
1163
1165 """
1166 Retrieve the bibtex entry of a catalog.
1167
1168 @param catalog: name of the catalog (e.g. II/306/sdss8)
1169 @type catalog: str
1170 @return: bibtex entry
1171 @rtype: str
1172 """
1173 bibcode = catalog2bibcode(catalog)
1174 bibtex = bibcode2bibtex(bibcode)
1175 return bibtex
1176
1177
1178
1179
1180
1181
1182 -def _get_URI(name=None,ID=None,ra=None,dec=None,radius=20.,
1183 oc='deg',oc_eq='J2000',
1184 out_all=True,out_max=1000000,
1185 filetype='tsv',sort='_r',constraints=None,**kwargs):
1186 """
1187 Build Vizier URI from available options.
1188
1189 kwargs are to catch unused arguments.
1190
1191 @param name: name of a ViZieR catalog (e.g. 'II/246/out')
1192 @type name: str
1193 @param filetype: type of the retrieved file
1194 @type filetype: str (one of 'tsv','csv','ascii'... see ViZieR site)
1195 @param oc: coordinates
1196 @type oc: str (one of 'deg'...)
1197 @param out_all: retrieve all or basic information
1198 @type out_all: boolean (True = all, None = basic)
1199 @param out_max: maximum number of rows
1200 @type out_max: integer
1201 @param ID: target name
1202 @type ID: str
1203 @param ra: target's right ascension
1204 @type ra: float
1205 @param dec: target's declination
1206 @type dec: float
1207 @param radius: search radius (arcseconds)
1208 @type radius: float
1209 @return: url
1210 @rtype: str
1211 """
1212 base_url = 'http://{0}/viz-bin/asu-{1}/VizieR?&-oc={2},eq={3}'.format(mirrors['current'],filetype,oc,oc_eq)
1213 if sort:
1214 base_url += '&-sort=%s'%(sort)
1215 if constraints is not None:
1216 for constr in constraints:
1217 base_url += '&%s'%(constr)
1218
1219 if ID is not None:
1220
1221
1222 if ID[0]=='J':
1223 ra = int(ID[1:3]),int(ID[3:5]),float(ID[5:10])
1224 dec = int(ID[10:13]),int(ID[13:15]),float(ID[15:])
1225 ra = '%02d+%02d+%.2f'%ra
1226 dec = '+%+02d+%02d+%.2f'%dec
1227 ID = None
1228
1229
1230 if name: base_url += '&-source=%s'%(name)
1231 if out_all: base_url += '&-out.all'
1232 if out_max: base_url += '&-out.max=%s'%(out_max)
1233 if radius: base_url += '&-c.rs=%s'%(radius)
1234 if ID is not None and ra is None: base_url += '&-c=%s'%(urllib.quote(ID))
1235 if ra is not None: base_url += '&-c.ra=%s&-c.dec=%s'%(ra,dec)
1236 logger.debug(base_url)
1237
1238 return base_url
1239
1241 """
1242 From colors and one magnitude measurement, derive the other magnitudes.
1243
1244 @param master: master record array from vizier2phot.
1245 @type master: record array
1246 @return: master with added magnitudes
1247 @rtype: record array
1248 """
1249 names = list(master.dtype.names)
1250 photbands = list(master['photband'])
1251 for i,photband in enumerate(photbands):
1252 system,color = photband.split('.')
1253
1254
1255
1256
1257 if '-' in color:
1258
1259 index_meas, index_emeas = names.index('meas'),names.index('e_meas')
1260 index_band = names.index('photband')
1261 row = list(master[i])
1262 meas,e_meas = row[index_meas],row[index_emeas]
1263
1264 band1,band2 = ['%s.%s'%(system,band) for band in color.split('-')]
1265 band1_present = band1 in photbands
1266 band2_present = band2 in photbands
1267
1268 if band1_present and not band2_present:
1269
1270 index1 = photbands.index(band1)
1271 row1 = list(master[index1])
1272 row1[index_meas] = row1[index_meas] - row[index_meas]
1273 errs = np.array([row[index_emeas],row1[index_emeas]],float)
1274
1275 if sum(np.isnan(errs))<len(errs):
1276 errs[np.isnan(errs)] = 0.
1277 row1[index_emeas]= np.sqrt(np.sum(errs**2))
1278 row1[index_band] = band2
1279 logger.debug("Added band %s = %s - %s (a)"%(band2,band1,photband))
1280 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype)
1281 elif band2_present and not band1_present:
1282
1283 index1 = photbands.index(band2)
1284 row1 = list(master[index1])
1285 row1[index_meas] = row[index_meas] + row1[index_meas]
1286 errs = np.array([row[index_emeas],row1[index_emeas]],float)
1287
1288 if sum(np.isnan(errs))<len(errs):
1289 errs[np.isnan(errs)] = 0.
1290 row1[index_emeas]= np.sqrt(np.sum(errs**2))
1291 row1[index_band] = band1
1292 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype)
1293 logger.debug("Added band %s = %s + %s (b)"%(band1,band2,photband))
1294
1295
1296
1297
1298
1299 elif color.upper()=='M1':
1300
1301
1302 index_meas, index_emeas = names.index('meas'),names.index('e_meas')
1303 index_band = names.index('photband')
1304
1305 row = list(master[i])
1306 meas,e_meas = row[index_meas],row[index_emeas]
1307
1308 my_photbands = list(master['photband'])
1309 if 'STROMGREN.Y' in my_photbands and 'STROMGREN.B' in my_photbands and (not 'STROMGREN.V' in my_photbands):
1310 index_b = my_photbands.index('STROMGREN.B')
1311 index_y = my_photbands.index('STROMGREN.Y')
1312 rowb,rowy = list(master[index_b]),list(master[index_y])
1313 b,e_b = rowb[index_meas],rowb[index_emeas]
1314 y,e_y = rowy[index_meas],rowy[index_emeas]
1315
1316 row1 = list(master[index_band])
1317 row1[index_band] = 'STROMGREN.V'
1318 row1[index_meas] = meas + 2*b - y
1319 row1[index_emeas] = np.sqrt(e_meas**2 + 2*e_b**2 + e_y**2)
1320 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype)
1321 logger.debug("Added band STROMGREN.Y (m1)")
1322
1323 elif color.upper()=='C1':
1324
1325
1326 index_meas, index_emeas = names.index('meas'),names.index('e_meas')
1327 index_band = names.index('photband')
1328
1329 row = list(master[i])
1330 meas,e_meas = row[index_meas],row[index_emeas]
1331
1332 my_photbands = list(master['photband'])
1333 if 'STROMGREN.V' in my_photbands and 'STROMGREN.B' in my_photbands and (not 'STROMGREN.U' in my_photbands):
1334 index_v = my_photbands.index('STROMGREN.V')
1335 index_b = my_photbands.index('STROMGREN.B')
1336 rowb,rowv = list(master[index_b]),list(master[index_v])
1337 b,e_b = rowb[index_meas],rowb[index_emeas]
1338 v,e_v = rowv[index_meas],rowv[index_emeas]
1339
1340 row1 = list(master[index_band])
1341 row1[index_band] = 'STROMGREN.U'
1342 row1[index_meas] = meas + 2*v - b
1343 row1[index_emeas] = np.sqrt(e_meas**2 + 2*e_v**2 + e_b**2)
1344 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype)
1345 logger.debug("Added band STROMGREN.U (c1)")
1346
1347
1348
1349
1350 return master
1351
1352
1354 """
1355 Execute all docstrings.
1356
1357 >>> import pylab
1358 >>> p = pylab.show()
1359 """
1360 import doctest
1361 doctest.testmod()
1362
1363
1364
1365 if __name__=="__main__":
1366 test()
1367