ivs.catalogs.vizier

1 # -*- coding: utf-8 -*- 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 #-- standard libraries 86 import numpy as np 87 import urllib 88 import logging 89 import os 90 import itertools 91 import astropy.io.fits as pf 92 import tarfile 93 import tempfile 94 import shutil 95 import ConfigParser 96 from scipy.spatial import KDTree 97 98 #-- IvS repository 99 from ivs.inout import ascii 100 from ivs.inout import fits 101 from ivs.inout import http 102 from ivs.units import conversions 103 from ivs.aux import loggers 104 from ivs.aux import numpy_ext 105 from ivs.sed import filters 106 from ivs import config 107 108 logger = logging.getLogger("CAT.VIZIER") 109 logger.addHandler(loggers.NullHandler()) 110 111 basedir = os.path.dirname(os.path.abspath(__file__)) 112 113 #-- read in catalog information 114 cat_info = ConfigParser.ConfigParser() 115 cat_info.optionxform = str # make sure the options are case sensitive 116 cat_info.readfp(open(os.path.join(basedir,'vizier_cats_phot.cfg'))) 117 118 cat_info_fund = ConfigParser.ConfigParser() 119 cat_info_fund.optionxform = str # make sure the options are case sensitive 120 cat_info_fund.readfp(open(os.path.join(basedir,'vizier_cats_fund.cfg'))) 121 122 mirrors = {'cycle': itertools.cycle(['vizier.u-strasbg.fr', # France 123 'vizier.nao.ac.jp', # Japan 124 'vizier.hia.nrc.ca', # Canada 125 'vizier.ast.cam.ac.uk', # UK 126 'vizier.cfa.harvard.edu', # USA CFA/harvard 127 'www.ukirt.jach.hawaii.edu',# USA Ukirt 128 'vizier.inasan.ru', # Russia 129 'vizier.iucaa.ernet.in', # India 130 'data.bao.ac.cn'])} # China 131 mirrors['current'] = mirrors['cycle'].next() 132 133 #{ Basic interfaces 134

135 -def change_mirror():

136 """ 137 Cycle through the mirrors of ViZieR. 138 """ 139 mirrors['current'] = mirrors['cycle'].next() 140 logger.info("Changed cycle to {}".format(mirrors['current']))

141

142 -def search(name,filetype='tsv',filename=None,**kwargs):

143 """ 144 Search and retrieve information from a VizieR catalog. 145 146 Two ways to search for data within a catalog C{name}: 147 148 1. You're looking for info on B{one target}, then give the target's 149 C{ID} or coordinates (C{ra} and C{dec}), and a search C{radius}. 150 151 2. You're looking for information of B{a whole field}, then give the 152 field's coordinates (C{ra} and C{dec}), and C{radius}. 153 154 If you have a list of targets, you need to loop this function. 155 156 If you supply a filename, the results will be saved to that path, and you 157 will get the filename back as received from urllib.URLopener (should be the 158 same as the input name, unless something went wrong). 159 160 If you don't supply a filename, you should leave C{filetype} to the default 161 C{tsv}, and the results will be saved to a temporary 162 file and deleted after the function is finished. The content of the file 163 will be read into a dictionary, as well as the units (two separate 164 dictionaries with the same keys, depending on the column names in the 165 catalog). The entries in the dictionary are of type C{ndarray}, and will 166 be converted to a float-array (no integers, we need to support nans) if 167 possible. If not, the array will consist of strings. The comments are also 168 returned as a list of strings. 169 170 WARNING: when retrieving a FITS file, ViZieR sometimes puts weird formats 171 into the header ('3F10.6E' in the 2MASS catalog), which cannot be read by 172 the C{astropy.io.fits} module. These columns are actually multi-dimensional vectors. 173 One option is to download to another format, or to restrict the columns with 174 C{out_all=None}. 175 176 Example usage: 177 178 1. Look for the Geneva V magnitude and coordinates of Vega in the GENEVA 179 catalog of Rufener. 180 181 >>> results,units,comms = search('II/169/main',ID='vega',radius=60.) 182 >>> print "Vega: Vmag = %.3f %s, RA = %.2f %s, DEC = %.2f %s"%(results['Vmag'],units['Vmag'],results['_RAJ2000'],units['_RAJ2000'],results['_DEJ2000'],units['_DEJ2000']) 183 Vega: Vmag = 0.061 mag, RA = 279.24 deg, DEC = 38.77 deg 184 185 2. Search for all targets in the 2MASS catalog in a particular field. 186 Download the results to a FITS file, read the file, and plot the results 187 to the screen. 188 189 >>> #filename = search('II/246/out',ra=100.79,dec=0.70,radius=1000.,filetype='fits',filename='2mass_test',out_all=None) 190 191 Now read in the FITS-file and plot the contents 192 193 >>> #import pylab 194 >>> #import astropy.io.fits as pf 195 >>> #ff = pf.open('2mass_test.fits') 196 >>> #p = pylab.gcf().canvas.set_window_title('test of <search>') 197 >>> #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') 198 >>> #p = pylab.colorbar() 199 >>> #p = p.set_label('Jmag') 200 >>> #p,q = pylab.xlabel('RA [deg]'),pylab.ylabel('DEC [deg]') 201 >>> #ff.close() 202 >>> #os.remove('2mass_test.fits') 203 204 205 @param name: name of a ViZieR catalog (e.g. 'II/246/out') 206 @type name: str 207 @param filetype: type of the file to write the results to ('tsv' if no file desired) 208 @type filetype: string (one of 'tsv','fits','ascii','csv'... 209 @param filename: name of the file to write the results to (no extension) 210 @type filename: str 211 @return: filename / catalog data columns, units, comments 212 @rtype: str/ record array, dict, list of str 213 """ 214 #-- two ways of giving filename: either complete filename with extension, 215 # or filename without extension, but C{filetype} speficied. 216 if filename is not None and '.' in os.path.basename(filename): 217 filetype = os.path.splitext(filename)[1][1:] 218 elif filename is not None: 219 filename = '%s.%s'%(filename,filetype) 220 221 #-- gradually build URI 222 base_url = _get_URI(name=name,**kwargs) 223 224 #-- prepare to open URI 225 url = urllib.URLopener() 226 filen,msg = url.retrieve(base_url,filename=filename) 227 # maybe we are just interest in the file, not immediately in the content 228 if filename is not None: 229 logger.info('Querying ViZieR source %s and downloading to %s'%(name,filen)) 230 url.close() 231 return filen 232 233 # otherwise, we read everything into a dictionary 234 if filetype=='tsv': 235 try: 236 results,units,comms = tsv2recarray(filen) 237 #-- raise an exception when multiple catalogs were specified 238 except ValueError: 239 raise ValueError, "failed to read %s, perhaps multiple catalogs specified (e.g. III/168 instead of III/168/catalog)"%(name) 240 url.close() 241 logger.info('Querying ViZieR source %s (%d)'%(name,(results is not None and len(results) or 0))) 242 return results,units,comms

243 244

245 -def list_catalogs(ID,filename=None,filetype='tsv',**kwargs):

246 """ 247 Print and return all catalogs containing information on the star. 248 249 If you give C{filetype} and C{filename}, all information will be downloaded 250 to that file. 251 252 Extra kwargs: see L{_get_URI}. 253 254 @param ID: identification of the star 255 @type ID: str 256 @keyword filetype: type of the output file ('fits','tsv','csv'...) 257 @type filetype: str 258 @keyword filename: name of the output file 259 @type filename: str 260 @return: dictionary with keys the catalog ID from VizieR and entries the 261 names of the catalogs 262 @rtype: dictionary 263 """ 264 base_url = _get_URI(ID=ID,filetype='fits',**kwargs) 265 266 #-- download the file 267 url = urllib.URLopener() 268 filen,msg = url.retrieve(base_url,filename=filename) 269 270 #-- if it is a FITS file, we extract all catalog IDs. We download the 271 # individual catalogs to retrieve their title. 272 if filetype=='fits': 273 mycats = {} 274 ff = pf.open(filen) 275 for ext in range(1,len(ff)): 276 name = ff[ext].header['CDS-name'] 277 results,units,comms = search(name,ID=ID,**kwargs) 278 for line in comms: 279 if 'Title:' in line: 280 title = line.strip().split(':')[1] 281 break 282 mycats[name] = title 283 logger.info('%25s %s'%(name,title)) 284 285 photometry = [col for col in units.keys() if 'mag' in units[col]] 286 rv = [col for col in units.keys() if 'rv' in col.lower()] 287 vsini = [col for col in units.keys() if 'sin' in col.lower()] 288 sptype = [col for col in units.keys() if col.lower()=='sp' or col.lower()=='sptype'] 289 fund = [col for col in units.keys() if 'logg' in col.lower() or 'teff' in col.lower()] 290 291 ff.close() 292 url.close() 293 return mycats 294 else: 295 url.close() 296 return filen

297 298

299 -def xmatch(source1,source2,output_file=None,tol=1.,**kwargs):

300 """ 301 Crossmatch two vizier catalogs via a fast KDTree. 302 303 The limit for these catalogs is probably somewhere between ~100000 entries, 304 so make sure your catalogs do not contain to many targets. You can always 305 do a subselection via the keyword arguments (e.g. give ra, dec and radius). 306 307 An output tsv file will be written (by default named 'source1__source2', 308 which can be read in via C{tsv2recarray} for further analysis. 309 310 tolerance is in arcseconds. 311 312 Extra keyword arguments are passed to C{search}. Column names of second 313 source will be appended with postfix '_2', to avoid clashes of double-defined 314 column headers. 315 """ 316 #-- construct default filename. 317 if output_file is None: 318 output_file = "__".join([source1,source2]).replace('/','_').replace('+','')+'.tsv' 319 320 #-- download the two catalogs 321 cat1,units1,comms1 = search(source1,**kwargs) 322 cat2,units2,comms2 = search(source2,**kwargs) 323 324 logger.info('Start Vizier Xmatch') 325 coords1 = np.array([cat1['_RAJ2000'],cat1['_DEJ2000']]).T 326 coords2 = np.array([cat2['_RAJ2000'],cat2['_DEJ2000']]).T 327 328 logger.info('Building KDTree of shape %d,%d'%coords1.shape) 329 tree = KDTree(coords1) 330 331 logger.info('Querying KDTree with %d entries'%(len(coords2))) 332 distance,order = tree.query(coords2) 333 334 keep = distance<(tol/(60.)) 335 336 logger.info('Matched %d points (tol<%.3g arcsec)'%(sum(keep),tol)) 337 338 #-- this the subset matching both catalogues 339 cat1 = cat1[order[keep]] 340 cat2 = cat2[keep] 341 342 #-- now write it to a vizier-like file... 343 #-- first append '2' to each column name of the second source in the 344 # comments, to make sure there are no doubles. 345 for i,line in enumerate(comms2): 346 line = line.split('\t') 347 if line[0]=='Column': 348 line[1] = line[1]+'_2' 349 comms2[i] = '\t'.join(line) 350 #-- change unit dictionaries and dtypes 351 units2_ = {} 352 for key in units2: units2_[key+'_2'] = units2[key] 353 cat2.dtype.names = [name+'_2' for name in cat2.dtype.names] 354 355 356 ff = open(output_file,'w') 357 ff.write('\n#'.join(comms1)) 358 ff.write('\n#'.join(comms2)) 359 360 names1 = list(cat1.dtype.names) 361 names2 = list(cat2.dtype.names) 362 dtypes = [(name,cat1.dtype[names1.index(name)].str) for name in names1] 363 dtypes += [(name,cat2.dtype[names2.index(name)].str) for name in names2] 364 365 ff.write('\n') 366 for nr,i in enumerate(dtypes): 367 ff.write(str(i[0])) 368 if nr<(len(dtypes)-1): ff.write('\t') 369 ff.write('\n') 370 for nr,i in enumerate(dtypes): 371 if i[0] in units1: ff.write(units1[i[0]]) 372 elif i[0] in units2_: ff.write(units2_[i[0]]) 373 else: 374 raise ValueError,'this cannot be' 375 if nr<(len(dtypes)-1): ff.write('\t') 376 377 ff.write('\n') 378 ff.write('\t'.join(['---']*len(dtypes))) 379 ff.write('\n') 380 381 for row1,row2 in itertools.izip(cat1,cat2): 382 ff.write('\t'.join([str(x) for x in row1])+'\t') 383 ff.write('\t'.join([str(x) for x in row2])+'\n') 384 385 ff.close()

386 387 388 389 390 #} 391 392 #{ Interface to specific catalogs 393

394 -def get_IUE_spectra(ID=None,directory=None,unzip=True,cat_info=False,select='low',**kwargs):

395 """ 396 Download IUE spectra. 397 398 If you want to download all the spectral files, set C{directory='/home/user/'} 399 or whatever. All the tarfiles will be downloaded to this directory, they 400 will be untarred, science data extracted and all unnecessary files and 401 directories will be deleted. If you don't set a directory, it will default 402 to the CWD. 403 404 If you don't wish to unzip them, set unzip=False 405 406 DEPRECATED: If you don't give a directory, the function will return a list 407 of all extracted spectra (no data files are kept). 408 409 You can retrieve the contents of the vizier catalog via {cat_info=True}. The 410 files will not be downloaded in this case. 411 """ 412 if directory is None: 413 direc = os.getcwd() 414 directory = os.getcwd() 415 filename = None 416 else: 417 direc = directory 418 if not os.path.isdir(direc): 419 os.mkdir(direc) 420 421 output = [] 422 #-- construct the download link form the camera and image data 423 data,units,comments = search('VI/110/inescat/',ID=ID,**kwargs) 424 425 if cat_info: 426 return data,units,comments 427 428 if data is None: 429 return output 430 431 for spectrum in data: 432 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'])) 433 logger.info('IUE spectrum %s/%s: %s'%(spectrum['Camera'],spectrum['Image'],download_link)) 434 435 #-- prepare to download the spectra to a temparorary file 436 if directory is not None: 437 filename = download_link.split('iue_mark=')[1].split('&')[0] 438 filename = os.path.join(direc,filename) 439 #-- download file and retrieve the path to the downloaded file 440 mytarfile = http.download(download_link,filename=filename) 441 if filename is None: 442 mytarfile,url = mytarfile 443 #-- perhaps unzipping is not necessary 444 if directory is not None and not unzip: 445 output.append(mytarfile) 446 url.close() 447 continue 448 #-- else unpack tar files and copy spectra to the directory, 449 # remember the location of the files: 450 if not tarfile.is_tarfile(mytarfile): 451 logger.info("Not a tarfile: %s"%(mytarfile)) 452 continue 453 tarf = tarfile.open(mytarfile) 454 files = tarf.getmembers(),tarf.getnames() 455 deldirs = [] 456 outfile = None 457 for mem,name in zip(*files): 458 #-- first check if it is a spectrum, but skip low or high res if needed. 459 myname = name.lower() 460 skip = False 461 if 'lo' in select and not 'lo' in os.path.basename(myname): 462 skip = True 463 if 'hi' in select and not 'hi' in os.path.basename(myname): 464 skip = true 465 if ('lwp' in name or 'swp' in name or 'lwr' in name) and not skip: 466 #-- first unpack this file 467 tarf.extract(mem,path=direc) 468 #-- move it to the direc directory 469 outfile = os.path.join(direc,os.path.basename(name)) 470 dirname = os.path.dirname(name) 471 shutil.move(os.path.join(direc,name),outfile) 472 logger.info("Extracted %s to %s"%(name,outfile)) 473 #-- remove any left over empty directory 474 deldirs.append(os.path.join(direc,dirname)) 475 else: 476 logger.debug("Did not extract %s (probably not a spectrum)"%(name)) 477 #-- remove the tar file: 478 tarf.close() 479 os.unlink(mytarfile) 480 for dirname in deldirs: 481 if dirname and os.path.isdir(dirname) and not os.listdir(dirname): 482 os.rmdir(dirname) 483 logger.debug("Deleted left over (empty) directory %s"%(dirname)) 484 if filename is None: url.close() 485 486 #-- only read in the data if they need to be extracted 487 if directory is not None and outfile: 488 output.append(outfile) 489 continue 490 if outfile and os.path.isfile(outfile): 491 wavelength,flux,error,header = fits.read_iue(outfile,return_header=True) 492 os.unlink(outfile) 493 output.append([wavelength,flux,error,header]) 494 else: 495 logger.info('Unsuccesfull extraction of %s'%(outfile)) 496 497 return output

498

499 -def get_UVSST_spectrum(units='erg/s/cm2/AA',**kwargs):

500 """ 501 Get a spectrum from the UVSST spectrograph onboard TD1. 502 503 From vizier catalog III/39A. 504 505 Also have a look at II/86/suppl. 506 """ 507 kwargs.setdefault('out_max',10) 508 kwargs.setdefault('radius',60.) 509 data,units_o,comments = search('III/39A/catalog',**kwargs) 510 if data is None: return None,None 511 fields = sorted([field for field in data.dtype.names if (field[0]=='F' and len(field)==5)]) 512 spectrum = np.zeros((len(fields),3)) 513 units_ = [] 514 for i,field in enumerate(fields): 515 wave,flux = float(field[1:]),data[field][0] 516 if flux==0: continue 517 #-- there is an error in the units ***in vizier***! 518 flux = conversions.convert(units_o[field],units,100*flux,wave=(wave,'AA'),unpack=True) 519 e_flux = data['e_'+field][0]/100.*flux 520 spectrum[i][0] = wave 521 spectrum[i][1] = flux 522 spectrum[i][2] = e_flux 523 units_.append(units_o[field]) 524 return spectrum[spectrum[:,0]>0]

525 #} 526 527 528 #{ Convenience functions

529 -def get_photometry(ID=None,extra_fields=['_r','_RAJ2000','_DEJ2000'],take_mean=False,**kwargs):

530 """ 531 Download all available photometry from a star to a record array. 532 533 For extra kwargs, see L{_get_URI} and L{vizier2phot} 534 535 """ 536 kwargs['ID'] = ID 537 to_units = kwargs.pop('to_units','erg/s/cm2/AA') 538 sources = kwargs.get('sources',cat_info.sections()) 539 master_ = kwargs.get('master',None) 540 master = None 541 #-- retrieve all measurements 542 for source in sources: 543 results,units,comms = search(source,**kwargs) 544 if results is None: continue 545 master = vizier2phot(source,results,units,master,extra_fields=extra_fields,take_mean=take_mean) 546 #-- convert the measurement to a common unit. 547 if to_units and master is not None: 548 #-- prepare columns to extend to basic master 549 dtypes = [('cwave','f8'),('cmeas','f8'),('e_cmeas','f8'),('cunit','a50')] 550 cols = [[],[],[],[]] 551 #-- forget about 'nan' errors for the moment 552 no_errors = np.isnan(master['e_meas']) 553 master['e_meas'][no_errors] = 0. 554 #-- extend basic master 555 zp = filters.get_info(master['photband']) 556 for i in range(len(master)): 557 to_units_ = to_units+'' 558 try: 559 value,e_value = conversions.convert(master['unit'][i],to_units,master['meas'][i],master['e_meas'][i],photband=master['photband'][i]) 560 except ValueError: # calibrations not available, or its a color 561 # if it is a magnitude color, try converting it to a flux ratio 562 if 'mag' in master['unit'][i]: 563 try: 564 value,e_value = conversions.convert('mag_color','flux_ratio',master['meas'][i],master['e_meas'][i],photband=master['photband'][i]) 565 to_units_ = 'flux_ratio' 566 except ValueError: 567 value,e_value = np.nan,np.nan 568 # else, we are powerless... 569 else: 570 value,e_value = np.nan,np.nan 571 try: 572 eff_wave = filters.eff_wave(master['photband'][i]) 573 except IOError: 574 eff_wave = np.nan 575 cols[0].append(eff_wave) 576 cols[1].append(value) 577 cols[2].append(e_value) 578 cols[3].append(to_units_) 579 master = numpy_ext.recarr_addcols(master,cols,dtypes) 580 #-- reset errors 581 master['e_meas'][no_errors] = np.nan 582 master['e_cmeas'][no_errors] = np.nan 583 584 if master_ is not None and master is not None: 585 master = numpy_ext.recarr_addrows(master_,master.tolist()) 586 elif master is None: 587 master = master_ 588 589 #-- and return the results 590 return master

591 592

593 -def quality_check(master,ID=None,return_master=True,**kwargs):

594 """ 595 Perform quality checks on downloaded data. 596 597 This function translates flags in to words, and looks up additional 598 information in selected catalogs. 599 """ 600 messages = ['' for i in range(len(master))] 601 #-- for some, we can do a quality check given the information that is 602 # already available in the master record 603 #-- IRAS 604 logger.info('Checking flags') 605 for i,entry in enumerate(master): 606 if 'USNO' in str(entry['photband']): 607 messages[i] = '; '.join([messages[i],'unreliable zeropoint and transmission profile']) 608 continue 609 if 'COUSINS' in str(entry['photband']): 610 messages[i] = '; '.join([messages[i],'unreliable zeropoint']) 611 continue 612 flag = entry['flag'] 613 try: 614 flag = float(flag) 615 except: 616 continue 617 if np.isnan(flag): continue 618 if entry['source'].strip() in ['II/125/main','II/275/fsr']: 619 flag = float(int(flag)) 620 if flag==3: messages[i] = '; '.join([messages[i],'high quality']) 621 if flag==2: messages[i] = '; '.join([messages[i],'moderate quality']) 622 if flag==1: messages[i] = '; '.join([messages[i],'upper limit']) 623 if entry['source'].strip() in ['II/298/fis','II/297/irc']: 624 flag = float(int(flag)) 625 if flag==3: messages[i] = '; '.join([messages[i],'high quality']) 626 if flag==2: messages[i] = '; '.join([messages[i],'source is confirmed but the flux is not reliable']) 627 if flag==1: messages[i] = '; '.join([messages[i],'the source is not confirmed']) 628 if flag==0: messages[i] = '; '.join([messages[i],'not observed']) 629 if 'DIRBE' in entry['source'].strip(): 630 flag = '{0:03d}'.format(int(float(flag))) 631 if flag[0]=='1': messages[i] = '; '.join([messages[i],'IRAS/2MASS companion greater than DIRBE noise level']) 632 if flag[1]=='1': messages[i] = '; '.join([messages[i],'possibly extended emission or highly variable source']) # discrepancy between DIRBE and IRAS/2MASS flux density 633 if flag[2]=='1': messages[i] = '; '.join([messages[i],'possibly affected by nearby companion']) 634 if entry['source'].strip() in ['V/114/msx6_main']: 635 if flag==4: messages[i] = '; '.join([messages[i],'excellent']) 636 if flag==3: messages[i] = '; '.join([messages[i],'good']) 637 if flag==2: messages[i] = '; '.join([messages[i],'fair']) 638 if flag==1: messages[i] = '; '.join([messages[i],'on limit']) 639 if flag==0: messages[i] = '; '.join([messages[i],'not detected']) 640 641 #-- for other targets, we need to query additional information 642 sources_with_quality_check = ['B/denis/denis','II/311/wise','II/246/out'] 643 sources = set(list(master['source'])) & set(sources_with_quality_check) 644 645 denis_image_flags = {'01':'clouds during observation', 646 '02':'electronic Read-out problem', 647 '04':'internal temperature problem', 648 '08':'very bright star', 649 '10':'bright star', 650 '20':'stray light', 651 '40':'unknown problem'} 652 denis_source_flags = {'01':'source might be a dust on mirror', 653 '02':'source is a ghost detection of a bright star', 654 '04':'source is saturated', 655 '08':'source is multiple detect', 656 '10':'reserved'} 657 wise_conf_flag = {'d':'diffraction spike', 658 'p':'contaminated by latent image left by bright star', 659 'h':'halo of nearby bright source', 660 'o':'optical ghost'} 661 wise_var_flag = {'n':'too few measurements to decide if variable', 662 '0':'most likely not variable (0, 0=certainly not-5=probably not)', 663 '1':'most likely not variable (1, 0=certainly not-5=probably not)', 664 '2':'most likely not variable (2, 0=certainly not-5=probably not)', 665 '3':'most likely not variable (3, 0=certainly not-5=probably not)', 666 '4':'most likely not variable (4, 0=certainly not-5=probably not)', 667 '5':'most likely not variable (5, 0=certainly not-5=probably not)', 668 '6':'likely variable (6, 6=likely-7=more likely)', 669 '7':'likely variable (7, 6=likely-7=more likely)', 670 '8':'most likely variable (8, 8=most likely-9=almost certain)', 671 '9':'most likely variable (9, 8=most likely-9=almost certain)'} 672 twomass_qual_flag = {'X':'detection, but no valid brightness estimate (X)', 673 'U':'upper limit (U)', 674 'F':'error estimate not reliable (F)', 675 'E':'poor PSF fit (E)', 676 'A':'high quality (A)', 677 'B':'high quality (B)', 678 'C':'high quality (C)', 679 'D':'high quality (D)'} 680 681 682 indices = np.arange(len(master)) 683 logger.info('Checking source catalogs for additional information') 684 for source in sorted(sources): 685 results,units,comms = search(source,ID=ID,**kwargs) 686 if results is None: continue 687 #-- the DENIS database 688 if source=='B/denis/denis': 689 for iflag,photband in zip(['Iflg','Jflg','Kflg'],['DENIS.I','DENIS.J','DENIS.KS']): 690 index = indices[(master['source']==source) & (master['photband']==photband)] 691 if not len(index): 692 continue 693 flag = float(results[0][iflag]) 694 if np.isnan(flag): 695 messages[index[0]] = '; '.join([messages[index[0]],'high quality']) 696 continue 697 flag = '{0:04d}'.format(int(flag)) 698 image_flag = flag[:2] 699 source_flag = flag[2:] 700 #-- keep track for output 701 if image_flag in denis_image_flags: 702 messages[index[0]] = '; '.join([messages[index[0]],denis_image_flags[image_flag]]) 703 if source_flag in denis_source_flags: 704 messages[index[0]] = '; '.join([messages[index[0]],denis_source_flags[source_flag]]) 705 if source=='II/311/wise': 706 conf = results[0]['ccf'] 707 var = results[0]['var'] 708 ex = int(results[0]['ex']) 709 for i,photband in enumerate(['WISE.W1','WISE.W2','WISE.W3','WISE.W4']): 710 index = indices[(master['source']==source) & (master['photband']==photband)] 711 if len(index)!=1: 712 logger.warning("Skipping WISE flags, don't know what to do with {}".format(index)) 713 continue 714 if conf[i]!=' ' and conf[i] in wise_conf_flag: 715 messages[index[0]] = '; '.join([messages[index[0]],wise_conf_flag[conf[i].lower()]]) 716 if var[i]!=' ': 717 messages[index[0]] = '; '.join([messages[index[0]],wise_var_flag[var[i].lower()]]) 718 messages[index[0]] = '; '.join([messages[index[0]],(ex==0 and 'point source' or 'extended source')]) 719 if source=='II/246/out': 720 flag = results[0]['Qflg'].strip() 721 for i,photband in enumerate(['2MASS.J','2MASS.H','2MASS.KS']): 722 if flag[i] in twomass_qual_flag: 723 index = indices[(master['source']==source) & (master['photband']==photband)] 724 if not len(index): 725 continue 726 messages[index[0]] = '; '.join([messages[index[0]],twomass_qual_flag[flag[i]]]) 727 728 729 #-- strip first '; ': 730 for i in range(len(messages)): 731 if messages[i]: 732 messages[i] = messages[i][2:] 733 else: 734 messages[i] = '-' 735 messages[i] = messages[i].replace(' ','_') 736 737 #-- perhaps we want to return the master record array with an extra column 738 messages = np.rec.fromarrays([messages],names=['comments']) 739 if return_master and not 'comments' in master.dtype.names: 740 return numpy_ext.recarr_join(master,messages) 741 elif return_master: 742 master['comments'] = messages 743 return master 744 else: 745 return messages

746 747 748 749 750 751 752

753 -def tsv2recarray(filename):

754 """ 755 Read a Vizier tsv (tab-sep) file into a record array. 756 757 @param filename: name of the TSV file 758 @type filename: str 759 @return: catalog data columns, units, comments 760 @rtype: record array, dict, list of str 761 """ 762 data,comms = ascii.read2array(filename,dtype=np.str,splitchar='\t',return_comments=True) 763 results = None 764 units = {} 765 #-- retrieve the data and put it into a record array 766 if len(data)>0: 767 #-- now convert this thing into a nice dictionary 768 data = np.array(data) 769 #-- retrieve the format of the columns. They are given in the 770 # Fortran format. In rare cases, columns contain multiple values 771 # themselves (so called vectors). In those cases, we interpret 772 # the contents as a long string 773 formats = np.zeros_like(data[0]) 774 for line in comms: 775 line = line.split('\t') 776 if len(line)<3: continue 777 for i,key in enumerate(data[0]): 778 if key == line[1] and line[0]=='Column': # this is the line with information 779 formats[i] = line[2].replace('(','').replace(')','').lower() 780 if formats[i][0].isdigit(): formats[i] = 'a100' 781 elif 'f' in formats[i]: formats[i] = 'f8' # floating point 782 elif 'i' in formats[i]: formats[i] = 'f8' # integer, but make it float to contain nans 783 elif 'e' in formats[i]: formats[i] = 'f8' # exponential 784 #-- see remark about the nans a few lines down 785 if formats[i][0]=='a': 786 formats[i] = 'a'+str(int(formats[i][1:])+3) 787 #-- define dtypes for record array 788 dtypes = np.dtype([(i,j) for i,j in zip(data[0],formats)]) 789 #-- remove spaces or empty values 790 cols = [] 791 for i,key in enumerate(data[0]): 792 col = data[3:,i] 793 #-- fill empty values with nan, and make sure each string in the 794 # array is at least three spaces long (otherwise we cannot fit 795 # 'nan' in the row and the casting will not work properly in a 796 # few lines 797 cols.append([(row.isspace() or not row) and np.nan or 3*' '+row for row in col]) 798 #-- fix unit name 799 #if source in cat_info.sections() and cat_info.has_option(source,data[1,i]): 800 # units[key] = cat_info.get(source,data[1,i]) 801 #else: 802 units[key] = data[1,i] 803 #-- define columns for record array and construct record array 804 cols = [np.cast[dtypes[i]](cols[i]) for i in range(len(cols))] 805 results = np.rec.array(cols, dtype=dtypes) 806 return results,units,comms

807

808 -def vizier2phot(source,results,units,master=None,e_flag='e_',q_flag='q_',extra_fields=None,take_mean=False):

809 """ 810 Convert/combine VizieR record arrays to measurement record arrays. 811 812 Every line in the combined array represents a measurement in a certain band. 813 814 This is probably only useful if C{results} contains only information on 815 one target (or you have to give 'ID' as an extra field, maybe). 816 817 The standard columns are: 818 819 1. C{meas}: containing the photometric measurement 820 2. C{e_meas}: the error on the photometric measurement 821 3. C{flag}: an optional quality flag 822 4. C{unit}: the unit of the measurement 823 5. C{photband}: the photometric passband (FILTER.BAND) 824 6. C{source}: name of the source catalog 825 826 You can add extra information from the VizieR catalog via the list of keys 827 C{extra_fields}. 828 829 If you give a C{master}, the information will be added to a previous 830 record array. If not, a new master will be created. 831 832 Colors will be expanded, derived from the other columns and added to the 833 master. 834 835 The result is a record array with each row a measurement. 836 837 Example usage: 838 839 First look for all photometry of Vega in all VizieR catalogs: 840 841 >>> from ivs.sed import filters 842 >>> import pylab 843 >>> master = None 844 >>> for source in cat_info.sections(): 845 ... results,units,comms = search(source,ID='vega',radius=60.) 846 ... if results is not None: 847 ... master = vizier2phot(source,results,units,master,extra_fields=['_r','_RAJ2000','_DEJ2000']) 848 849 Keep only observations we have an measurement and error of, convert every 850 observation to 'Jy' and keep track of the results to plot. 851 852 >>> master = master[(-np.isnan(master['e_meas'])) & (-np.isnan(master['meas']))] 853 >>> eff_waves = filters.eff_wave(master['photband']) 854 >>> myvalue,e_myvalue = conversions.nconvert(master['unit'],'erg/s/cm2/AA',master['meas'],master['e_meas'],photband=master['photband']) 855 >>> for i in range(len(master)): 856 ... 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']) 857 GENEVA.V 6.100e-02+/- 2.500e-02 mag 3.620e-09 Jy 279.24 38.77 56.000 II/169/main 858 GENEVA.B -8.980e-01+/- 2.500e-02 mag 6.518e-09 Jy 279.24 38.77 56.000 II/169/main 859 GENEVA.U 6.070e-01+/- 2.500e-02 mag 3.223e-09 Jy 279.24 38.77 56.000 II/169/main 860 GENEVA.V1 7.640e-01+/- 2.500e-02 mag 3.782e-09 Jy 279.24 38.77 56.000 II/169/main 861 GENEVA.B1 2.000e-03+/- 2.500e-02 mag 6.584e-09 Jy 279.24 38.77 56.000 II/169/main 862 GENEVA.B2 6.120e-01+/- 2.500e-02 mag 6.208e-09 Jy 279.24 38.77 56.000 II/169/main 863 GENEVA.G 1.270e+00+/- 2.500e-02 mag 3.111e-09 Jy 279.24 38.77 56.000 II/169/main 864 2MASS.J -1.770e-01+/- 2.060e-01 mag 3.591e-10 Jy 279.23 38.78 0.034 II/246/out 865 2MASS.H -2.900e-02+/- 1.460e-01 mag 1.176e-10 Jy 279.23 38.78 0.034 II/246/out 866 2MASS.KS 1.290e-01+/- 1.860e-01 mag 3.799e-11 Jy 279.23 38.78 0.034 II/246/out 867 IRAS.F12 4.160e+01+/- 1.664e+00 Jy 1.024e-13 Jy 279.23 38.78 9.400 II/125/main 868 IRAS.F25 1.100e+01+/- 5.500e-01 Jy 6.195e-15 Jy 279.23 38.78 9.400 II/125/main 869 IRAS.F60 9.510e+00+/- 7.608e-01 Jy 8.420e-16 Jy 279.23 38.78 9.400 II/125/main 870 IRAS.F100 7.760e+00+/- 6.984e-01 Jy 2.349e-16 Jy 279.23 38.78 9.400 II/125/main 871 TD1.1565 5.689e-09+/- 1.700e-11 10mW/m2/nm 5.689e-09 Jy 279.23 38.78 18.510 II/59B/catalog 872 TD1.1965 4.928e-09+/- 1.300e-11 10mW/m2/nm 4.928e-09 Jy 279.23 38.78 18.510 II/59B/catalog 873 TD1.2365 3.700e-09+/- 1.000e-11 10mW/m2/nm 3.700e-09 Jy 279.23 38.78 18.510 II/59B/catalog 874 TD1.2740 3.123e-09+/- 9.000e-12 10mW/m2/nm 3.123e-09 Jy 279.23 38.78 18.510 II/59B/catalog 875 ANS.15N -4.910e-01+/- 1.000e-03 mag 5.707e-09 Jy 279.23 38.78 14.000 II/97/ans 876 ANS.15W -4.410e-01+/- 1.200e-02 mag 5.450e-09 Jy 279.23 38.78 14.000 II/97/ans 877 ANS.18 -4.620e-01+/- 3.000e-03 mag 5.556e-09 Jy 279.23 38.78 14.000 II/97/ans 878 ANS.25 4.600e-02+/- 4.000e-03 mag 3.480e-09 Jy 279.23 38.78 14.000 II/97/ans 879 ANS.33 1.910e-01+/- 3.000e-03 mag 3.045e-09 Jy 279.23 38.78 14.000 II/97/ans 880 HIPPARCOS.HP 8.680e-02+/- 2.100e-03 mag 3.840e-09 Jy 279.23 38.78 3.060 I/239/hip_main 881 MIPS.24 8.900e+03+/- 8.900e+01 mJy 4.628e-15 Jy 279.23 38.78 0.010 J/ApJ/653/675/table1 882 MIPS.70 1.142e+04+/- 2.283e+03 mJy 7.075e-16 Jy 279.23 38.78 0.010 J/ApJ/653/675/table1 883 JOHNSON.B 1.900e-02+/- 1.000e-02 mag 6.216e-09 Jy 279.23 38.78 3.060 I/280B/ascc 884 JOHNSON.V 7.400e-02+/- 2.000e-03 mag 3.428e-09 Jy 279.23 38.78 3.060 I/280B/ascc 885 JOHNSON.V 3.300e-02+/- 1.200e-02 mag 3.560e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans 886 JOHNSON.B-V -1.000e-03+/- 5.000e-03 mag nan Jy 279.23 38.78 0.010 II/168/ubvmeans 887 JOHNSON.U-B -6.000e-03+/- 6.000e-03 mag nan Jy 279.23 38.78 0.010 II/168/ubvmeans 888 JOHNSON.B 3.200e-02+/- 1.300e-02 mag 6.142e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans 889 JOHNSON.U 2.600e-02+/- 1.432e-02 mag 4.086e-09 Jy 279.23 38.78 0.010 II/168/ubvmeans 890 JOHNSON.K 1.300e-01+/- 1.900e-01 mag 3.764e-11 Jy 279.23 38.78 0.010 J/PASP/120/1128/catalog 891 AKARI.N60 6.582e+00+/- 2.090e-01 Jy 4.614e-16 Jy 279.23 38.78 3.400 II/298/fis 892 AKARI.WIDES 6.201e+00+/- 1.650e-01 Jy 2.566e-16 Jy 279.23 38.78 3.400 II/298/fis 893 AKARI.WIDEL 4.047e+00+/- 3.500e-01 Jy 5.658e-17 Jy 279.23 38.78 3.400 II/298/fis 894 AKARI.N160 3.221e+00+/- 2.550e-01 Jy 3.695e-17 Jy 279.23 38.78 3.400 II/298/fis 895 AKARI.S9W 5.670e+01+/- 4.010e-01 Jy 2.169e-13 Jy 279.24 38.78 2.550 II/297/irc 896 AKARI.L18W 1.254e+01+/- 7.770e-02 Jy 1.050e-14 Jy 279.24 38.78 2.550 II/297/irc 897 STROMGREN.B-Y 3.000e-03+/- 3.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog 898 STROMGREN.M1 1.570e-01+/- 3.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog 899 STROMGREN.C1 1.088e+00+/- 4.000e-03 mag nan Jy 279.23 38.78 22.000 II/215/catalog 900 STROMGREN.B 4.300e-02+/- 3.000e-03 mag 5.604e-09 Jy 279.23 38.78 22.000 II/215/catalog 901 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 902 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 903 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 904 DIRBE.F12 2.910e+01+/- 1.570e+01 Jy 6.014e-14 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 905 DIRBE.F25 1.630e+01+/- 3.120e+01 Jy 1.153e-14 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 906 DIRBE.F60 1.220e+01+/- 5.610e+01 Jy 1.195e-15 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 907 DIRBE.F100 -7.000e-01+/- 7.790e+01 Jy -2.238e-17 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 908 DIRBE.F140 2.557e+02+/- 5.223e+03 Jy 3.577e-15 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 909 DIRBE.F240 8.290e+01+/- 2.881e+03 Jy 4.152e-16 Jy 279.23 38.78 0.110 J/ApJS/154/673/DIRBE 910 WISE.W2 1.143e+00+/- 1.900e-02 mag 8.428e-13 Jy 279.24 38.78 3.276 II/311/wise 911 WISE.W3 -6.700e-02+/- 8.000e-03 mag 6.930e-14 Jy 279.24 38.78 3.276 II/311/wise 912 WISE.W4 -1.270e-01+/- 6.000e-03 mag 5.722e-15 Jy 279.24 38.78 3.276 II/311/wise 913 914 Make a quick plot: 915 916 >>> p = pylab.figure() 917 >>> p = pylab.loglog(eff_waves,myvalue,'ko') 918 >>> p = pylab.show() 919 920 @param source: name of the VizieR source 921 @type source: str 922 @param results: results from VizieR C{search} 923 @type results: record array 924 @param units: header of Vizier catalog with key name the column name and 925 key value the units of that column 926 @type units: dict 927 @param master: master record array to add information to 928 @type master: record array 929 @param e_flag: flag denoting the error on a column 930 @type e_flag: str 931 @param q_flag: flag denoting the quality of a measurement 932 @type q_flag: str 933 @param extra_fields: any extra columns you want to add information from 934 @type extra_fields: list of str 935 @return: array with each row a measurement 936 @rtype: record array 937 """ 938 if cat_info.has_option(source,'q_flag'): 939 q_flag = cat_info.get(source,'q_flag') 940 if cat_info.has_option(source,'e_flag'): 941 e_flag = cat_info.get(source,'e_flag') 942 943 #-- basic dtypes 944 dtypes = [('meas','f8'),('e_meas','f8'),('flag','a20'), 945 ('unit','a30'),('photband','a30'),('source','a50')] 946 947 #-- extra can be added: 948 names = list(results.dtype.names) 949 950 if extra_fields is not None: 951 for e_dtype in extra_fields: 952 if e_dtype in names: 953 dtypes.append((e_dtype,results.dtype[names.index(e_dtype)].str)) 954 else: 955 dtypes.append((e_dtype,'f8')) 956 957 #-- create empty master if not given 958 newmaster = False 959 if master is None or len(master)==0: 960 master = np.rec.array([tuple([('f' in dt[1]) and np.nan or 'none' for dt in dtypes])],dtype=dtypes) 961 newmaster = True 962 963 #-- add fluxes and magnitudes to the record array 964 cols_added = 0 965 for key in cat_info.options(source): 966 if key in ['e_flag','q_flag','mag','bibcode']: 967 continue 968 photband = cat_info.get(source,key) 969 key = key.replace('_[','[').replace(']_',']') 970 #-- contains measurement, error, quality, units, photometric band and 971 # source 972 if not take_mean or len(results)<=1: #take first 973 cols = [results[key][:1], 974 (e_flag+key in results.dtype.names and results[e_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 975 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 976 np.array(len(results[:1])*[units[key]],str), 977 np.array(len(results[:1])*[photband],str), 978 np.array(len(results[:1])*[source],str)] 979 else: 980 logger.warning('taking mean: {0} --> {1}+/-{2}'.format(results[key][0],results[key].mean(),results[key].std())) 981 cols = [[results[key].mean()],[results[key].std()], 982 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 983 np.array(len(results[:1])*[units[key]],str), 984 np.array(len(results[:1])*[photband],str), 985 np.array(len(results[:1])*[source],str)] 986 #-- correct errors given in percentage 987 if e_flag+key in results.dtype.names and units[e_flag+key]=='%': 988 cols[1] = cols[1]/100.*cols[0] 989 #-- add any extra fields if desired. 990 if extra_fields is not None: 991 for e_dtype in extra_fields: 992 cols += [e_dtype in results.dtype.names and results[e_dtype][:1] or np.ones(len(results[:1]))*np.nan] 993 #-- add to the master 994 rows = [] 995 for i in range(len(cols[0])): 996 rows.append(tuple([col[i] for col in cols])) 997 master = np.core.records.fromrecords(master.tolist()+rows,dtype=dtypes) 998 cols_added += 1 999 #-- fix colours: we have to run through it two times to be sure to have 1000 # all the colours 1001 N = len(master)-cols_added 1002 master_ = _breakup_colours(master[N:]) 1003 master_ = _breakup_colours(master_) 1004 #-- combine and return 1005 master = np.core.records.fromrecords(master.tolist()[:N]+master_.tolist(),dtype=dtypes) 1006 1007 #-- skip first line from building 1008 if newmaster: master = master[1:] 1009 return master

1010 1011 1012

1013 -def vizier2fund(source,results,units,master=None,e_flag='e_',q_flag='q_',extra_fields=None):

1014 """ 1015 Convert/combine VizieR record arrays to measurement record arrays. 1016 1017 This is probably only useful if C{results} contains only information on 1018 one target (or you have to give 'ID' as an extra field, maybe). 1019 1020 The standard columns are: 1021 1022 1. C{meas}: containing the measurement of a fundamental parameter 1023 2. C{e_meas}: the error on the measurement of a fundamental parameter 1024 3. C{flag}: an optional quality flag 1025 4. C{unit}: the unit of the measurement 1026 5. C{source}: name of the source catalog 1027 1028 If a target appears more than once in a catalog, only the first match will 1029 be added. 1030 1031 The result is a record array with each row a measurement. 1032 1033 Example usage: 1034 1035 >>> master = None 1036 >>> for source in cat_info_fund.sections(): 1037 ... results,units,comms = search(source,ID='AzV 79',radius=60.) 1038 ... if results is not None: 1039 ... master = vizier2fund(source,results,units,master,extra_fields=['_r','_RAJ2000','_DEJ2000']) 1040 >>> master = master[-np.isnan(master['meas'])] 1041 >>> for i in range(len(master)): 1042 ... 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']) 1043 Teff 7.304e+03+/- nan K 12.67 -72.83 0.002 B/pastel/pastel 1044 logg 2.000e+00+/- nan [cm/s2] 12.67 -72.83 0.002 B/pastel/pastel 1045 [Fe/H] -8.700e-01+/- nan [Sun] 12.67 -72.83 0.002 B/pastel/pastel 1046 1047 @param source: name of the VizieR source 1048 @type source: str 1049 @param results: results from VizieR C{search} 1050 @type results: record array 1051 @param units: header of Vizier catalog with key name the column name and 1052 key value the units of that column 1053 @type units: dict 1054 @param master: master record array to add information to 1055 @type master: record array 1056 @param e_flag: flag denoting the error on a column 1057 @type e_flag: str 1058 @param q_flag: flag denoting the quality of a measurement 1059 @type q_flag: str 1060 @param extra_fields: any extra columns you want to add information from 1061 @type extra_fields: list of str 1062 @return: array with each row a measurement 1063 @rtype: record array 1064 """ 1065 if cat_info_fund.has_option(source,'q_flag'): 1066 q_flag = cat_info_fund.get(source,'q_flag') 1067 if cat_info_fund.has_option(source,'e_flag'): 1068 e_flag = cat_info_fund.get(source,'e_flag') 1069 1070 #-- basic dtypes 1071 dtypes = [('meas','f8'),('e_meas','f8'),('q_meas','f8'),('unit','a30'), 1072 ('source','a50'),('name','a50')] 1073 1074 #-- extra can be added: 1075 names = list(results.dtype.names) 1076 if extra_fields is not None: 1077 for e_dtype in extra_fields: 1078 dtypes.append((e_dtype,results.dtype[names.index(e_dtype)].str)) 1079 1080 #-- create empty master if not given 1081 newmaster = False 1082 if master is None or len(master)==0: 1083 master = np.rec.array([tuple([('f' in dt[1]) and np.nan or 'none' for dt in dtypes])],dtype=dtypes) 1084 newmaster = True 1085 1086 #-- add fluxes and magnitudes to the record array 1087 for key in cat_info_fund.options(source): 1088 if key in ['e_flag','q_flag']: 1089 continue 1090 photband = cat_info_fund.get(source,key) 1091 #-- contains measurement, error, quality, units, photometric band and 1092 # source 1093 #if len(results[e_flag+key])>1: 1094 key = key.replace('_[','[').replace(']_',']') 1095 cols = [results[key][:1], 1096 (e_flag+key in results.dtype.names and results[e_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 1097 (q_flag+key in results.dtype.names and results[q_flag+key][:1] or np.ones(len(results[:1]))*np.nan), 1098 np.array(len(results[:1])*[units[key]],str), 1099 np.array(len(results[:1])*[source],str), 1100 np.array(len(results[:1])*[key],str)] 1101 #-- add any extra fields if desired. 1102 if extra_fields is not None: 1103 for e_dtype in extra_fields: 1104 cols.append(results[:1][e_dtype]) 1105 #-- add to the master 1106 rows = [] 1107 for i in range(len(cols[0])): 1108 rows.append(tuple([col[i] for col in cols])) 1109 #print master 1110 master = np.core.records.fromrecords(master.tolist()+rows,dtype=dtypes) 1111 1112 #-- skip first line from building 1113 if newmaster: master = master[1:] 1114 return master

1115 1116

1117 -def catalog2bibcode(catalog):

1118 """ 1119 Retrieve the ADS bibcode of a ViZieR catalog. 1120 1121 @param catalog: name of the catalog (e.g. II/306/sdss8) 1122 @type catalog: str 1123 @return: bibtex code 1124 @rtype: str 1125 """ 1126 catalog = catalog.split("/") 1127 N = max(2,len(catalog)-1) 1128 catalog = "/".join(catalog[:N]) 1129 base_url = "http://cdsarc.u-strasbg.fr/viz-bin/Cat?{0}".format(catalog) 1130 url = urllib.URLopener() 1131 filen,msg = url.retrieve(base_url) 1132 1133 code = None 1134 ff = open(filen,'r') 1135 for line in ff.readlines(): 1136 if 'Keywords' in line: break 1137 if '=<A HREF' in line: 1138 code = line.split('>')[1].split('<')[0] 1139 ff.close() 1140 url.close() 1141 return code

1142

1143 -def bibcode2bibtex(bibcode):

1144 """ 1145 Retrieve the bibtex entry of an ADS bibcode. 1146 1147 @param bibcode: bibcode (e.g. C{2011yCat.2306....0A}) 1148 @type bibcode: str 1149 @return: bibtex entry 1150 @rtype: str 1151 """ 1152 base_url = "http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode={0}&data_type=BIBTEX&db_key=AST&nocookieset=1".format(bibcode) 1153 url = urllib.URLopener() 1154 filen,msg = url.retrieve(base_url) 1155 1156 bibtex = [] 1157 ff = open(filen,'r') 1158 for line in ff.readlines(): 1159 if (not bibtex and '@' in line) or bibtex: 1160 bibtex.append(line) 1161 ff.close() 1162 url.close() 1163 return "".join(bibtex)

1164

1165 -def catalog2bibtex(catalog):

1166 """ 1167 Retrieve the bibtex entry of a catalog. 1168 1169 @param catalog: name of the catalog (e.g. II/306/sdss8) 1170 @type catalog: str 1171 @return: bibtex entry 1172 @rtype: str 1173 """ 1174 bibcode = catalog2bibcode(catalog) 1175 bibtex = bibcode2bibtex(bibcode) 1176 return bibtex

1177 1178 1179 #} 1180 1181 #{ Internal helper functions 1182

1183 -def _get_URI(name=None,ID=None,ra=None,dec=None,radius=20., 1184 oc='deg',oc_eq='J2000', 1185 out_all=True,out_max=1000000, 1186 filetype='tsv',sort='_r',constraints=None,**kwargs):

1187 """ 1188 Build Vizier URI from available options. 1189 1190 kwargs are to catch unused arguments. 1191 1192 @param name: name of a ViZieR catalog (e.g. 'II/246/out') 1193 @type name: str 1194 @param filetype: type of the retrieved file 1195 @type filetype: str (one of 'tsv','csv','ascii'... see ViZieR site) 1196 @param oc: coordinates 1197 @type oc: str (one of 'deg'...) 1198 @param out_all: retrieve all or basic information 1199 @type out_all: boolean (True = all, None = basic) 1200 @param out_max: maximum number of rows 1201 @type out_max: integer 1202 @param ID: target name 1203 @type ID: str 1204 @param ra: target's right ascension 1205 @type ra: float 1206 @param dec: target's declination 1207 @type dec: float 1208 @param radius: search radius (arcseconds) 1209 @type radius: float 1210 @return: url 1211 @rtype: str 1212 """ 1213 base_url = 'http://{0}/viz-bin/asu-{1}/VizieR?&-oc={2},eq={3}'.format(mirrors['current'],filetype,oc,oc_eq) 1214 if sort: 1215 base_url += '&-sort=%s'%(sort) 1216 if constraints is not None: 1217 for constr in constraints: 1218 base_url += '&%s'%(constr) 1219 1220 if ID is not None: 1221 #-- if the ID is given in the form 'J??????+??????', derive the 1222 # coordinates of the target from the name. 1223 if ID[0]=='J': 1224 ra = int(ID[1:3]),int(ID[3:5]),float(ID[5:10]) 1225 dec = int(ID[10:13]),int(ID[13:15]),float(ID[15:]) 1226 ra = '%02d+%02d+%.2f'%ra 1227 dec = '+%+02d+%02d+%.2f'%dec 1228 ID = None 1229 1230 1231 if name: base_url += '&-source=%s'%(name) 1232 if out_all: base_url += '&-out.all' 1233 if out_max: base_url += '&-out.max=%s'%(out_max) 1234 if radius: base_url += '&-c.rs=%s'%(radius) 1235 if ID is not None and ra is None: base_url += '&-c=%s'%(urllib.quote(ID)) 1236 if ra is not None: base_url += '&-c.ra=%s&-c.dec=%s'%(ra,dec) 1237 logger.debug(base_url) 1238 #print base_url 1239 return base_url

1240

1241 -def _breakup_colours(master):

1242 """ 1243 From colors and one magnitude measurement, derive the other magnitudes. 1244 1245 @param master: master record array from vizier2phot. 1246 @type master: record array 1247 @return: master with added magnitudes 1248 @rtype: record array 1249 """ 1250 names = list(master.dtype.names) 1251 photbands = list(master['photband']) 1252 for i,photband in enumerate(photbands): 1253 system,color = photband.split('.') 1254 1255 ######################################################################## 1256 #-- NORMAL COLORS 1257 ######################################################################## 1258 if '-' in color: # we have a colour 1259 #-- in which column are the measurements (and error) located? 1260 index_meas, index_emeas = names.index('meas'),names.index('e_meas') 1261 index_band = names.index('photband') 1262 row = list(master[i]) 1263 meas,e_meas = row[index_meas],row[index_emeas] 1264 1265 band1,band2 = ['%s.%s'%(system,band) for band in color.split('-')] 1266 band1_present = band1 in photbands 1267 band2_present = band2 in photbands 1268 1269 if band1_present and not band2_present: 1270 #-- which row do we need to compute the component of the colour? 1271 index1 = photbands.index(band1) 1272 row1 = list(master[index1]) 1273 row1[index_meas] = row1[index_meas] - row[index_meas] 1274 errs = np.array([row[index_emeas],row1[index_emeas]],float) 1275 #-- allow for nan errors if all errors on the photbands are nan 1276 if sum(np.isnan(errs))<len(errs): 1277 errs[np.isnan(errs)] = 0. 1278 row1[index_emeas]= np.sqrt(np.sum(errs**2)) 1279 row1[index_band] = band2 1280 logger.debug("Added band %s = %s - %s (a)"%(band2,band1,photband)) 1281 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype) 1282 elif band2_present and not band1_present: 1283 #-- which row do we need to compute the component of the colour? 1284 index1 = photbands.index(band2) 1285 row1 = list(master[index1]) 1286 row1[index_meas] = row[index_meas] + row1[index_meas] 1287 errs = np.array([row[index_emeas],row1[index_emeas]],float) 1288 #-- allow for nan errors if all errors on the photbands are nan 1289 if sum(np.isnan(errs))<len(errs): 1290 errs[np.isnan(errs)] = 0. 1291 row1[index_emeas]= np.sqrt(np.sum(errs**2)) 1292 row1[index_band] = band1 1293 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype) 1294 logger.debug("Added band %s = %s + %s (b)"%(band1,band2,photband)) 1295 1296 ######################################################################## 1297 #-- STROMGREN COLORS 1298 ######################################################################## 1299 #-- stromgren index M1 1300 elif color.upper()=='M1': 1301 # m1 = v - 2*b + y 1302 #-- in which column are the measurements (and error) located? 1303 index_meas, index_emeas = names.index('meas'),names.index('e_meas') 1304 index_band = names.index('photband') 1305 #-- this is the m1 row 1306 row = list(master[i]) 1307 meas,e_meas = row[index_meas],row[index_emeas] 1308 #-- retrieve the measurements we need to calculate band-magnitudes 1309 my_photbands = list(master['photband']) 1310 if 'STROMGREN.Y' in my_photbands and 'STROMGREN.B' in my_photbands and (not 'STROMGREN.V' in my_photbands): 1311 index_b = my_photbands.index('STROMGREN.B') 1312 index_y = my_photbands.index('STROMGREN.Y') 1313 rowb,rowy = list(master[index_b]),list(master[index_y]) 1314 b,e_b = rowb[index_meas],rowb[index_emeas] 1315 y,e_y = rowy[index_meas],rowy[index_emeas] 1316 #-- add extra row 1317 row1 = list(master[index_band]) 1318 row1[index_band] = 'STROMGREN.V' 1319 row1[index_meas] = meas + 2*b - y 1320 row1[index_emeas] = np.sqrt(e_meas**2 + 2*e_b**2 + e_y**2) 1321 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype) 1322 logger.debug("Added band STROMGREN.Y (m1)") 1323 #-- stromgren index C1 1324 elif color.upper()=='C1': 1325 # c1 = u - 2*v + b 1326 #-- in which column are the measurements (and error) located? 1327 index_meas, index_emeas = names.index('meas'),names.index('e_meas') 1328 index_band = names.index('photband') 1329 #-- this is the m1 row 1330 row = list(master[i]) 1331 meas,e_meas = row[index_meas],row[index_emeas] 1332 #-- retrieve the measurements we need to calculate band-magnitudes 1333 my_photbands = list(master['photband']) 1334 if 'STROMGREN.V' in my_photbands and 'STROMGREN.B' in my_photbands and (not 'STROMGREN.U' in my_photbands): 1335 index_v = my_photbands.index('STROMGREN.V') 1336 index_b = my_photbands.index('STROMGREN.B') 1337 rowb,rowv = list(master[index_b]),list(master[index_v]) 1338 b,e_b = rowb[index_meas],rowb[index_emeas] 1339 v,e_v = rowv[index_meas],rowv[index_emeas] 1340 #-- add extra row 1341 row1 = list(master[index_band]) 1342 row1[index_band] = 'STROMGREN.U' 1343 row1[index_meas] = meas + 2*v - b 1344 row1[index_emeas] = np.sqrt(e_meas**2 + 2*e_v**2 + e_b**2) 1345 master = np.core.records.fromrecords(master.tolist()+[tuple(row1)],dtype=master.dtype) 1346 logger.debug("Added band STROMGREN.U (c1)") 1347 1348 1349 1350 1351 return master

1352 1353

1354 -def test():

1355 """ 1356 Execute all docstrings. 1357 1358 >>> import pylab 1359 >>> p = pylab.show() 1360 """ 1361 import doctest 1362 doctest.testmod()

1363 1364 #} 1365 1366 if __name__=="__main__": 1367 test() 1368

Source Code for Module ivs.catalogs.vizier