CSP DR1+DR2

sndatasets/loaders.py

@@ -9,18 +9,23 @@
 from os.path import join
 
 import numpy as np
+import sncosmo
 from astropy.io import ascii
-from astropy.table import Table
+from astropy.table import Table, vstack
+from astropy.coordinates import SkyCoord
+import tarfile
 
 from .dlutils import download_file, download_sn_positions
 from .utils import (hms_to_deg, sdms_to_deg, pivot_table,
-                    mag_to_flux, jd_to_mjd, sxhr_to_deg, sx_to_deg)
+                    mag_to_flux, jd_to_mjd, sxhr_to_deg, sx_to_deg,
+                    cmb_to_helio)
 
 
 CDS_PREFIX = "http://cdsarc.u-strasbg.fr/vizier/ftp/cats/"
 
 
-__all__ = ["load_kowalski08", "load_hamuy96", "load_krisciunas"]
+__all__ = ["load_kowalski08", "load_hamuy96", "load_krisciunas",
+           "load_csp"]
 
 
 def load_kowalski08():
@@ -45,13 +50,14 @@ def load_kowalski08():
 
     data = pivot_table(data, 'band', ['{}mag', 'e_{}mag'],
                        ['B', 'V', 'R', 'I'])
+
     data = data[data['mag'] != 0.]  # eliminate missing values
 
     # Join telescope and band into one column
     data['band'] = np.char.add(np.char.replace(data['Tel'], ' ', '_'),
                                np.char.add('_', data['band']))
     del data['Tel']
-
+    
     # Split up table into one table per SN and add metadata.
     sne = OrderedDict()
     for i in range(len(meta)):
@@ -140,8 +146,7 @@ def load_hamuy96():
                               ('ra', sxhr_to_deg(meta[name][0])),
                               ('dec', sx_to_deg(meta[name][1]))])
 
-        sndata = data[data['SN'] == name]
-
+        sndata = data[data['SN'] == name]            
         time = jd_to_mjd(sndata['HJD'])
         band = np.char.add('bessell', np.char.lower(sndata['band']))
         zp = 29. * np.ones(len(sndata), dtype=np.float64)
@@ -248,3 +253,117 @@ def load_krisciunas():
                           meta=snmeta)
 
     return sne
+
+
+def load_csp():
+    """CSP DR1 + DR2 Photometry from Contreras et al 2010 and Stritzinger
+    et al 2011
+
+    http://adsabs.harvard.edu/abs/2010AJ....139..519C
+    http://adsabs.harvard.edu/abs/2011AJ....142..156S
+    """
+
+    dataurl = 'http://csp.obs.carnegiescience.edu/data/CSP_Photometry_DR2.tar.gz'
+    tarball = download_file(dataurl, 'csp')
+    tar = tarfile.open(tarball)
+    path = '/'.join(tarball.split('/')[:-1])
+    tar.extractall(path=path)
+    names = filter(lambda fn: fn.endswith('dat') and not '._' in fn,
+                   tar.getnames())
+    names = map(lambda f: join(path, f), names)
+    tar.close()
+
+    # Keep track of the telescopes used for IR observations. 
+    ir1 = download_file(CDS_PREFIX + 'J/AJ/139/519/table5.dat', 'csp_ir1')
+    readme1 = download_file(CDS_PREFIX + 'J/AJ/139/519/ReadMe', 'csp_ir1')
+    ir2 = download_file(CDS_PREFIX + 'J/AJ/142/156/table5.dat', 'csp_ir2')
+    readme2 = download_file(CDS_PREFIX + 'J/AJ/142/156/ReadMe', 'csp_ir2')
+    tel1 = ascii.read(ir1, format='cds', readme=readme1)[['SN', 'JD', 'Tel']]
+    tel2 = ascii.read(ir2, format='cds', readme=readme2)[['SN', 'JD', 'Tel']]
+    tel2['JD'] += 2453000.  # JD column has an offset of 2453000
+    tel = vstack((tel1, tel2))
+    magsys = sncosmo.get_magsystem('csp')    
+
+    def _read_csp(f, **kwargs):
+
+        meta = OrderedDict()
+        data = []
+        colnames = ['mjd', 'filter', 'flux', 'fluxerr', 'zp', 'magsys']
+        readingdata = False
+
+        def _which_V(mjd):
+            # return the CSP V band that was in use on mjd.
+            if mjd < 53748:
+                ans = '3014'
+            elif mjd > 53761:
+                ans = '9844'
+            else:
+                ans = '3009'
+            return 'cspv' + ans
+
+        for j, line in enumerate(f):
+            if not readingdata:
+                if j == 4:
+                    filts = [n.lower() for n in line[1:].strip().split()[1:] if n != '+/-']
+                    readingdata = True
+
+                if j == 2:
+                    sline = line[1:].split()
+                    meta['z_cmb'] = float(sline[2])
+                    ra  = (sline[5].replace(':', '%s') + '%s') % ('h','m','s')
+                    dec = (sline[8].replace(':', '%s') + '%s') % ('d','m','s')
+
+                    # convert from dms to degrees
+                    coord = SkyCoord(ra, dec)
+
+                    meta['ra']  = coord.ra.value
+                    meta['dec'] = coord.dec.value
+                    meta['z_helio'] = cmb_to_helio(meta['z_cmb'],
+                                                   meta['ra'],
+                                                   meta['dec'])
+
+                if j == 0:
+                    meta['name'] = line.split()[-1]
+
+            else:
+                d = line.split()
+                mjd = float(d[0])
+                for i in range(1, len(d), 2):
+                    if d[i] != '99.900':
+                        if filts[i / 2] == 'v':
+                            # figure out which V
+                            filt = _which_V(mjd)
+                        else:
+                            filt = 'csp' + filts[i / 2]
+                        for b in ['y', 'j', 'h']:
+                            if b in filt:
+
+                                cond1 = tel['SN'] == meta['name'][2:]
+                                cond2 = np.isclose(jd_to_mjd(tel['JD']),mjd)
+                                cond = np.logical_and(cond1, cond2)
+
+                                if not cond.any():
+                                    filt += 's'  # the SN is not in the vizier catalog
+                                elif tel[cond][0]['Tel'].lower().startswith('d'):
+                                    filt += 'd'
+                                else:
+                                    filt += 's'
+                                break
+
+                        mag = float(d[i])
+                        magerr = float(d[i + 1])
+                        zpsys = 'csp'
+                        zp = 2.5 * np.log10(magsys.zpbandflux(filt))
+                        flux, fluxerr = mag_to_flux(mag, magerr, zp)
+                        data.append((mjd, filt, flux, fluxerr, zp, zpsys))
+
+        data = dict(zip(colnames, zip(*data)))
+        return Table(data, meta=meta)
+
+    sne = OrderedDict()
+    for name in names:
+        with open(name, 'r') as f:
+            snname = name.split('/')[-1][2:].split('opt')[0]
+            sne[snname] = _read_csp(f)
+    return sne
+

sndatasets/utils.py

@@ -3,6 +3,7 @@
 import math
 
 import numpy as np
+from numpy import float64
 from astropy.table import Table