From 2a47f8d82d03c06d26e0d8a9b1121f485a47a401 Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 15:23:10 -0700 Subject: [PATCH 1/8] start working on csp --- sndatasets/loaders.py | 78 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 74 insertions(+), 4 deletions(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index 769d7ac..dcdec76 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -9,8 +9,9 @@ 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, join as tabjoin from .dlutils import download_file, download_sn_positions from .utils import (hms_to_deg, sdms_to_deg, pivot_table, @@ -20,7 +21,8 @@ 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(): @@ -140,8 +142,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 +249,72 @@ def load_krisciunas(): meta=snmeta) return sne + +def load_csp(): + """CSP DR2 Photometry from Stritzinger et al 2011 + http://adsabs.harvard.edu/abs/2011AJ....142..156S """ + + prefix = CDS_PREFIX + 'J/AJ/142/156/' + readme = download_file(prefix + 'ReadMe', 'csp') + table1 = download_file(prefix + 'table1.dat', 'csp') # general + table4 = download_file(prefix + 'table4.dat', 'csp') # optical + table5 = download_file(prefix + 'table5.dat', 'csp') # ir + + meta = ascii.read(table1, format='cds', readme=readme) + ra = hms_to_deg(meta['RAh'], meta['RAm'], meta['RAs']) + dec = sdms_to_deg(meta['DE-'], meta['DEd'], meta['DEm'], meta['DEs']) + + opt_data = ascii.read(table4, format='cds', readme=readme) + ir_data = ascii.read(table5, format='cds', readme=readme) + data = tabjoin(opt_data, ir_data, join_type='outer') + + data = data.filled(0.) # copying this from kyle + + data = pivot_table(data, 'band', ['{}mag', 'e_{}mag'], + ['u', 'g', 'r', 'i', 'B', 'V', + 'Y', 'J', 'H', 'K']) + + data = data[data['mag'] != 0] # eliminate missing values + + bandtel = zip(data['band'], data['Tel']) + ir = ['Y', 'J', 'H'] + + data['filter'] = ['csp' + b if b not in ir else 'csp' + b + t[0] \ + for (b, t) in bandtel] + data['filter'] = [f.lower() for f in data['filter']] + del data['Tel'], data['band'], data['---'], data['f_JD'] + + 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 + + data['filter'] = [_which_V(jd_to_mjd(t)) if 'v' in b else b \ + for (b, t) in zip(data['filter'], data['JD'])] + + sne = OrderedDict() + magsys = sncosmo.get_magsystem('csp') + for i in range(len(meta)): + name = meta['SN'][i] + sndata = data[data['SN'] == name] + snmeta = OrderedDict([('name', name), + ('dataset', 'csp'), + ('z_helio', meta['z'][i]), + ('ra', ra[i]), + ('dec', dec[i])]) + zpsys = len(sndata) * ['csp'] + zp = [magsys.offsets[magsys.bands.index(sncosmo.get_bandpass(b))] \ + for b in data['filter']] + flux, fluxerr = mag_to_flux(sndata['mag'], sndata['e_mag'], zp) + sne[name] = Table([jd_to_mjd(sndata['JD']), sndata['filter'], + flux, fluxerr, zp, zpsys], + names=('time', 'band', 'flux', 'fluxerr', 'zp', + 'zpsys'), + meta=snmeta) + return sne + From 3fb4cacfda260a1f150642e65c0c4d1d9116bd5e Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 16:04:40 -0700 Subject: [PATCH 2/8] fix a few bugs --- sndatasets/loaders.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index dcdec76..b1292ee 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -47,6 +47,7 @@ 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 @@ -54,6 +55,8 @@ def load_kowalski08(): np.char.add('_', data['band'])) del data['Tel'] + data['e_mag'] /= 1e3 # convert from mmag to mag + # Split up table into one table per SN and add metadata. sne = OrderedDict() for i in range(len(meta)): @@ -292,10 +295,13 @@ def _which_V(mjd): ans = '9844' else: ans = '3009' - return 'cspv' + ans - - data['filter'] = [_which_V(jd_to_mjd(t)) if 'v' in b else b \ - for (b, t) in zip(data['filter'], data['JD'])] + return 'cspv' + ans + + data['JD'] += 2453000 + data['band'] = [_which_V(jd_to_mjd(t)) if 'v' in b else b \ + for (b, t) in zip(data['filter'], data['JD'])] + + del data['filter'] sne = OrderedDict() magsys = sncosmo.get_magsystem('csp') @@ -308,10 +314,9 @@ def _which_V(mjd): ('ra', ra[i]), ('dec', dec[i])]) zpsys = len(sndata) * ['csp'] - zp = [magsys.offsets[magsys.bands.index(sncosmo.get_bandpass(b))] \ - for b in data['filter']] + zp = [2.5 * np.log10(magsys.zpbandflux(b)) for b in sndata['band']] flux, fluxerr = mag_to_flux(sndata['mag'], sndata['e_mag'], zp) - sne[name] = Table([jd_to_mjd(sndata['JD']), sndata['filter'], + sne[name] = Table([jd_to_mjd(sndata['JD']), sndata['band'], flux, fluxerr, zp, zpsys], names=('time', 'band', 'flux', 'fluxerr', 'zp', 'zpsys'), From f35e8ce11c653799a564f7cfd1bbec2d853e7ecd Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 16:14:46 -0700 Subject: [PATCH 3/8] working version of DR2 --- sndatasets/loaders.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index b1292ee..2db387f 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -54,9 +54,7 @@ def load_kowalski08(): data['band'] = np.char.add(np.char.replace(data['Tel'], ' ', '_'), np.char.add('_', data['band'])) del data['Tel'] - - data['e_mag'] /= 1e3 # convert from mmag to mag - + # Split up table into one table per SN and add metadata. sne = OrderedDict() for i in range(len(meta)): @@ -287,6 +285,12 @@ def load_csp(): data['filter'] = [f.lower() for f in data['filter']] del data['Tel'], data['band'], data['---'], data['f_JD'] + # convert from mmag to mag + data['magerr'] = data['e_mag'].astype(float) / 1e3 + del data['e_mag'] + data['e_mag'] = data['magerr'] + del data['magerr'] + def _which_V(mjd): # return the CSP V band that was in use on mjd. if mjd < 53748: From 8891ef48895f94bf5be1b9367520c9d76708eb38 Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 19:44:42 -0700 Subject: [PATCH 4/8] fix reference bug to float64 in utils --- sndatasets/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sndatasets/utils.py b/sndatasets/utils.py index a883123..9e529e2 100644 --- a/sndatasets/utils.py +++ b/sndatasets/utils.py @@ -3,6 +3,7 @@ import math import numpy as np +from numpy import float64 from astropy.table import Table From 8c2f91e18f9641a802bb2c451d7c6bed772f3615 Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 20:01:24 -0700 Subject: [PATCH 5/8] working versions of CSP dr1 and dr2 --- sndatasets/loaders.py | 185 ++++++++++++++++++++++++++---------------- 1 file changed, 114 insertions(+), 71 deletions(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index 2db387f..18e497a 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -8,14 +8,19 @@ import os from os.path import join +import glob + import numpy as np import sncosmo from astropy.io import ascii -from astropy.table import Table, join as tabjoin +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/" @@ -251,79 +256,117 @@ def load_krisciunas(): return sne -def load_csp(): - """CSP DR2 Photometry from Stritzinger et al 2011 - http://adsabs.harvard.edu/abs/2011AJ....142..156S """ - - prefix = CDS_PREFIX + 'J/AJ/142/156/' - readme = download_file(prefix + 'ReadMe', 'csp') - table1 = download_file(prefix + 'table1.dat', 'csp') # general - table4 = download_file(prefix + 'table4.dat', 'csp') # optical - table5 = download_file(prefix + 'table5.dat', 'csp') # ir - - meta = ascii.read(table1, format='cds', readme=readme) - ra = hms_to_deg(meta['RAh'], meta['RAm'], meta['RAs']) - dec = sdms_to_deg(meta['DE-'], meta['DEd'], meta['DEm'], meta['DEs']) - - opt_data = ascii.read(table4, format='cds', readme=readme) - ir_data = ascii.read(table5, format='cds', readme=readme) - data = tabjoin(opt_data, ir_data, join_type='outer') - - data = data.filled(0.) # copying this from kyle - - data = pivot_table(data, 'band', ['{}mag', 'e_{}mag'], - ['u', 'g', 'r', 'i', 'B', 'V', - 'Y', 'J', 'H', 'K']) - data = data[data['mag'] != 0] # eliminate missing values - - bandtel = zip(data['band'], data['Tel']) - ir = ['Y', 'J', 'H'] - - data['filter'] = ['csp' + b if b not in ir else 'csp' + b + t[0] \ - for (b, t) in bandtel] - data['filter'] = [f.lower() for f in data['filter']] - del data['Tel'], data['band'], data['---'], data['f_JD'] - - # convert from mmag to mag - data['magerr'] = data['e_mag'].astype(float) / 1e3 - del data['e_mag'] - data['e_mag'] = data['magerr'] - del data['magerr'] - - 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 +def load_csp(): + """CSP DR1 + DR2 Photometry from Contreras et al 2010 and Stritzinger + et al 2011 - data['JD'] += 2453000 - data['band'] = [_which_V(jd_to_mjd(t)) if 'v' in b else b \ - for (b, t) in zip(data['filter'], data['JD'])] + http://adsabs.harvard.edu/abs/2010AJ....139..519C + http://adsabs.harvard.edu/abs/2011AJ....142..156S + """ - del data['filter'] + 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') % ('d','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 + print(1) + 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() - magsys = sncosmo.get_magsystem('csp') - for i in range(len(meta)): - name = meta['SN'][i] - sndata = data[data['SN'] == name] - snmeta = OrderedDict([('name', name), - ('dataset', 'csp'), - ('z_helio', meta['z'][i]), - ('ra', ra[i]), - ('dec', dec[i])]) - zpsys = len(sndata) * ['csp'] - zp = [2.5 * np.log10(magsys.zpbandflux(b)) for b in sndata['band']] - flux, fluxerr = mag_to_flux(sndata['mag'], sndata['e_mag'], zp) - sne[name] = Table([jd_to_mjd(sndata['JD']), sndata['band'], - flux, fluxerr, zp, zpsys], - names=('time', 'band', 'flux', 'fluxerr', 'zp', - 'zpsys'), - meta=snmeta) + 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 From 73baa3e0f71e777d935d5143d745a216ec9a2463 Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 20:03:34 -0700 Subject: [PATCH 6/8] remove unused glob --- sndatasets/loaders.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index 18e497a..9a7729a 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -8,8 +8,6 @@ import os from os.path import join -import glob - import numpy as np import sncosmo from astropy.io import ascii From bccdcc970ca9d3d350454eca7ddd5fcb7df0c89f Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Sun, 31 Jul 2016 20:16:23 -0700 Subject: [PATCH 7/8] remove print statement --- sndatasets/loaders.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index 9a7729a..5a1b0c9 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -344,7 +344,6 @@ def _which_V(mjd): if not cond.any(): filt += 's' # the SN is not in the vizier catalog - print(1) elif tel[cond][0]['Tel'].lower().startswith('d'): filt += 'd' else: From afe84f2056b120d2bc8d0decd40e6fd5f17349c9 Mon Sep 17 00:00:00 2001 From: Danny Goldstein Date: Wed, 24 Aug 2016 11:20:34 -0700 Subject: [PATCH 8/8] load ra and dec properly --- sndatasets/loaders.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sndatasets/loaders.py b/sndatasets/loaders.py index 5a1b0c9..360b4db 100644 --- a/sndatasets/loaders.py +++ b/sndatasets/loaders.py @@ -310,7 +310,7 @@ def _which_V(mjd): if j == 2: sline = line[1:].split() meta['z_cmb'] = float(sline[2]) - ra = (sline[5].replace(':', '%s') + '%s') % ('d','m','s') + ra = (sline[5].replace(':', '%s') + '%s') % ('h','m','s') dec = (sline[8].replace(':', '%s') + '%s') % ('d','m','s') # convert from dms to degrees