diff --git a/astroquery/alfa/__init__.py b/astroquery/alfa/__init__.py new file mode 100644 index 0000000000..23d87084a8 --- /dev/null +++ b/astroquery/alfa/__init__.py @@ -0,0 +1,8 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +ALFALFA Spectra Archive Query Tool +----------------------------------- + +:Author: Jordan Mirocha (mirochaj@gmail.com) +""" +from .core import * \ No newline at end of file diff --git a/astroquery/alfa/core.py b/astroquery/alfa/core.py new file mode 100644 index 0000000000..9ed1d70f1d --- /dev/null +++ b/astroquery/alfa/core.py @@ -0,0 +1,258 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" + +core.py + +Author: Jordan Mirocha +Affiliation: University of Colorado at Boulder +Created on: Fri May 3 09:45:13 2013 + +Description: + +""" + +import urllib +import numpy as np +import numpy.ma as ma +from astropy.io import fits +from astropy import coordinates as coord + +ascii_prefix = "http://arecibo.tc.cornell.edu/hiarchive/alfalfa/spectraASCII" +fits_prefix = "http://arecibo.tc.cornell.edu/hiarchive/alfalfa/spectraFITS" +propert_path = "http://egg.astro.cornell.edu/alfalfa/data/a40files/a40.datafile1.csv" + +placeholder = -999999 + +def get_catalog(): + """ + Download catalog of ALFALFA source properties. + + Notes + ----- + This catalog has ~15,000 entries, so after it's downloaded, it is made + global to save some time later. + + Returns + ------- + Dictionary of results, each element is a masked array. + + """ + + if 'ALFALFACAT' in globals(): + return ALFALFACAT + + f = urllib.urlopen(propert_path) + + # Read header + cols = f.readline().rstrip('\n').split(',') + + catalog = {} + for col in cols: + catalog[col] = [] + + # Parse result + for line in f: + l = line.rstrip('\n').split(',') + for i, col in enumerate(cols): + item = l[i].strip() + if item == '\"\"': + catalog[col].append(placeholder) + elif item.isdigit(): + catalog[col].append(int(item)) + elif item.replace('.', '').isdigit(): + catalog[col].append(float(item)) + else: + catalog[col].append(item) + + # Mask out blank elements + for col in cols: + mask = np.zeros(len(catalog[col])) + mask[np.array(catalog[col]) == placeholder] = 1 + catalog[col] = ma.array(catalog[col], mask=mask) + + # Make this globally available so we don't have to re-download it + # again in this session + global ALFALFACAT + ALFALFACAT = catalog + + return catalog + +def get_spectrum(agc=None, ra=None, dec=None, unit=None, counterpart=False, + ascii=False): + """ + Download spectrum from ALFALFA catalogue. + + Parameters + ---------- + agc : int + Identification number for object in ALFALFA catalog. + ra : float + Right ascension (degrees). + dec : float + Declination (degrees). + ascii : bool + Download spectrum from remote server in ASCII or FITS format? + counterpart : bool + Do supplied ra and dec refer to HI source or optical counterpart? + + Notes + ----- + If AGC number is not supplied, will download entire ALFALFA catalog and + do a cross-ID with supplied RA and DEC. + + Returns + ------- + If ascii == False, returns Spectrum instance, which contains FITS hdulist + and properties for x and y axes. If ascii == True, returns a tuple of 4 + arrays: velocity [km/s], frequency [MHz], flux density [mJy], and baseline + fit [mJy], respectively. + + See Also + -------- + get_catalog : method that downloads ALFALFA catalog + + """ + + if agc is not None: + agc = str(agc).zfill(6) + else: + if ra is None and dec is None: + raise ValueError('Must supply ra and dec if agc=None!') + + try: + cat = ALFALFACAT + except NameError: + cat = get_catalog() + + if not isinstance(ra, coord.angles.RA): + ra = coord.RA(ra, unit=unit) + if not isinstance(ra, coord.angles.Dec): + dec = coord.Dec(dec, unit=unit) + + # Use RA and DEC to find appropriate AGC + if counterpart: + ra_ref = cat['RAdeg_OC'] + dec_ref = cat['DECdeg_OC'] + else: + ra_ref = cat['RAdeg_HI'] + dec_ref = cat['Decdeg_HI'] + + dra = np.abs(ra_ref - ra.degrees) \ + * np.cos(dec.degrees * np.pi / 180.) + ddec = np.abs(dec_ref - dec.degrees) + dr = np.sqrt(dra**2 + ddec**2) + + agc = cat['AGCNr'][np.argmin(dr)] + + print 'Found HI source AGC #%i %g arcseconds from supplied position.' \ + % (agc, dr.min() * 3600.) + + if ascii: + link = "%s/A%s.txt" % (ascii_prefix, agc) + f = urllib.urlopen(link) + + vel = [] + freq = [] + flux = [] + baseline = [] + for i, line in enumerate(f): + if i <= 30: continue + + data = line.split() + + vel.append(float(data[0])) + freq.append(float(data[1])) + flux.append(float(data[2])) + baseline.append(float(data[3])) + + f.close() + + return np.array(vel), np.array(freq), np.array(flux), np.array(baseline) + + link = "%s/A%s.fits" % (fits_prefix, agc) + hdulist = fits.open(urllib.urlopen(link).url, ignore_missing_end=True) + return Spectrum(hdulist) + +def match_object(ra, dec, ra_ref, dec_ref): + """ + Assumes everything is in degrees. Supply ra and dec of single object being considered, as well as reference + arrays of RA and DEC for some sample. Returns index of match in reference arrays. + """ + + + if min(dr) < maxsep: + return np.argmin(dr) + else: + return placeholder + + +class Spectrum: + def __init__(self, hdulist): + self.hdulist = hdulist + + @property + def hdr(self): + if not hasattr(self, '_hdr'): + for item in self.hdulist: + if item.name == 'PRIMARY': + continue + break + + self._hdr = {} + for key in item.header.keys(): + self._hdr[key] = self.header.get(key) + + return self._hdr + + @property + def header(self): + return self.hdulist[0].header + + @property + def freq(self): + """ Return xarr in units of MHz. """ + if not hasattr(self, '_freq'): + for item in self.hdulist: + if item.name == 'PRIMARY': + continue + break + + for i, col in enumerate(item.columns): + if col.name == 'FREQ': + self._freq = item.data[0][i] + break + + return self._freq + + @property + def varr(self): + """ Return xarr in units of helio-centric velocity (km/s). """ + if not hasattr(self, '_varr'): + for item in self.hdulist: + if item.name == 'PRIMARY': + continue + break + + for i, col in enumerate(item.columns): + if col.name == 'VHELIO': + self._varr = item.data[0][i] + break + + return self._varr + + @property + def data(self): + if not hasattr(self, '_data'): + for item in self.hdulist: + if item.name == 'PRIMARY': + continue + break + + for i, col in enumerate(item.columns): + if col.name == 'FLUXDENS': + self._data = item.data[0][i] + break + + return self._data + + \ No newline at end of file diff --git a/astroquery/alfa/tests/__init__.py b/astroquery/alfa/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/astroquery/alfa/tests/test_alfa.py b/astroquery/alfa/tests/test_alfa.py new file mode 100644 index 0000000000..1066553250 --- /dev/null +++ b/astroquery/alfa/tests/test_alfa.py @@ -0,0 +1,16 @@ +from astroquery import alfa + +# Test Case: A Seyfert 1 galaxy +RA = '0h8m05.63s' +DEC = '14d50m23.3s' + +def test_alfa_catalog(): + cat = alfa.get_catalog() + +def test_alfa_spectrum(): + sp = alfa.get_spectrum(ra=RA, dec=DEC, counterpart=True) + +if __name__ == '__main__': + test_alfa_catalog() + test_alfa_spectrum() + \ No newline at end of file diff --git a/astroquery/sdss/__init__.py b/astroquery/sdss/__init__.py index 311403c6db..fe3bd70c93 100644 --- a/astroquery/sdss/__init__.py +++ b/astroquery/sdss/__init__.py @@ -5,4 +5,4 @@ :Author: Jordan Mirocha (mirochaj@gmail.com) """ -from core import * +from .core import * diff --git a/astroquery/sdss/core.py b/astroquery/sdss/core.py index dedfeba286..68295e233b 100644 --- a/astroquery/sdss/core.py +++ b/astroquery/sdss/core.py @@ -17,6 +17,7 @@ import astropy.wcs as wcs import os, re, math from astropy.io import fits +import re, sqlcl, math, urllib from astropy import coordinates as coord from . import sqlcl @@ -36,7 +37,7 @@ 'qso_bright': 32 } -# Some website prefixes we need +# Some website prefixes we need spectro1d_prefix = 'http://das.sdss.org/spectro/1d_26' images_prefix = 'http://das.sdss.org/www/cgi-bin/drC' template_prefix = 'http://www.sdss.org/dr5/algorithms/spectemplates/spDR2' @@ -162,10 +163,6 @@ def get_spectrum(crossID=None, plate=None, fiberID=None, mjd=None): well as the FITS header in dictionary form. """ - safe_to_rm = True - if os.path.exists('spectro'): - safe_to_rm = False - if crossID is not None: plate = crossID['plate'] fiberID = crossID['fiberID'] @@ -174,17 +171,11 @@ def get_spectrum(crossID=None, plate=None, fiberID=None, mjd=None): plate = str(plate).zfill(4) fiber = str(fiberID).zfill(3) mjd = str(mjd) - web = '%s/%s/1d/spSpec-%s-%s-%s.fit' % (spectro1d_prefix, plate, mjd, + link = '%s/%s/1d/spSpec-%s-%s-%s.fit' % (spectro1d_prefix, plate, mjd, plate, fiber) - - os.system('wget -x -nH -nv -q %s' % web) - - hdulist = fits.open('spectro/1d_26/%s/1d/spSpec-%s-%s-%s.fit' % (plate, - mjd, plate, fiber), ignore_missing_end=True) - - if safe_to_rm: - os.system('rm -rf spectro') - + + hdulist = fits.open(urllib.urlopen(link).url, ignore_missing_end=True) + return Spectrum(hdulist) def get_image(crossID=None, run=None, rerun=None, camcol=None, @@ -222,10 +213,6 @@ def get_image(crossID=None, run=None, rerun=None, camcol=None, header in dictionary form. """ - safe_to_rm = True - if os.path.exists('imaging') or os.path.exists('www'): - safe_to_rm = False - if crossID is not None: run = crossID['run'] rerun = crossID['rerun'] @@ -239,18 +226,9 @@ def get_image(crossID=None, run=None, rerun=None, camcol=None, # Download and read in image data link = '%s?RUN=%i&RERUN=%i&CAMCOL=%i&FIELD=%s&FILTER=%s' % (images_prefix, run, rerun, camcol, field, band) - path_to_img = 'www/cgi-bin/drC?RUN=%i&RERUN=%i&CAMCOL=%i&FIELD=%s&FILTER=%s' % (run, - rerun, camcol, field, band) - - os.system('wget -x -nH -nv -q \'%s\'' % link) - - hdulist = fits.open(path_to_img, ignore_missing_end=True) - - # Erase download directory tree - if safe_to_rm: - os.system('rm -rf www') - os.system('rm -rf imaging') - + + hdulist = fits.open(urllib.urlopen(link).url, ignore_missing_end=True) + return Image(hdulist) def get_spectral_template(kind='qso'): @@ -284,13 +262,8 @@ def get_spectral_template(kind='qso'): available for some spectral types. """ - safe_to_rm = True - if os.path.exists('dr5'): - safe_to_rm = False - if kind == 'all': indices = list(np.arange(33)) - else: indices = spec_templates[kind] if type(indices) is not list: @@ -299,15 +272,11 @@ def get_spectral_template(kind='qso'): spectra = [] for index in indices: name = str(index).zfill(3) - web = '%s-%s.fit' % (template_prefix, name) - os.system('wget -x -nH -nv -q %s' % web) - hdulist = fits.open('dr5/algorithms/spectemplates/spDR2-%s.fit' % name) + link = '%s-%s.fit' % (template_prefix, name) + hdulist = fits.open(urllib.urlopen(link).url, ignore_missing_end=True) spectra.append(Spectrum(hdulist)) del hdulist - if safe_to_rm: - os.system('rm -rf dr5') - return spectra class Spectrum: diff --git a/astroquery/sdss/tests/test_sdss.py b/astroquery/sdss/tests/test_sdss.py index 4ed08b52af..751ad484b0 100644 --- a/astroquery/sdss/tests/test_sdss.py +++ b/astroquery/sdss/tests/test_sdss.py @@ -12,7 +12,11 @@ def test_sdss_image(): xid = sdss.crossID(ra=RA, dec=DEC) img = sdss.get_image(crossID=xid[0]) +def test_sdss_template(): + template = sdss.get_spectral_template('qso') + if __name__ == '__main__': test_sdss_spectrum() test_sdss_image() + test_sdss_template() \ No newline at end of file