In [1]:
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
import os
import numpy as np
import copy
In [2]:
from OpSimSummary import summarize_opsim as oss
import sncosmo
In [3]:
#import seaborn as sns
#sns.set()
In [4]:
from metrics import PerSNMetric
from efficiencyTable import EfficiencyTable
In [5]:
from lsst.sims.photUtils import BandpassDict
In [6]:
from lsst.sims.utils.coordinateTransformations import haversine
In [7]:
# Catsim bandpasses
lsst_bp = BandpassDict.loadTotalBandpassesFromFiles()
In [8]:
# sncosmo Bandpasses required for fitting
throughputsdir = os.getenv('THROUGHPUTS_DIR')
from astropy.units import Unit
bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
for band in bandPassList:
# setup sncosmo bandpasses
bandfname = banddir + "/total_" + band + '.dat'
# register the LSST bands to the SNCosmo registry
# Not needed for LSST, but useful to compare independent codes
# Usually the next two lines can be merged,
# but there is an astropy bug currently which affects only OSX.
numpyband = np.loadtxt(bandfname)
print band
sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
trans=numpyband[:, 1],
wave_unit=Unit('nm'),
name=band)
sncosmo.registry.register(sncosmoband, force=True)
In [9]:
def distance(ra1, dec1, ra2, dec2):
"""
Obtain the angular distance in radians between two points with given
ra, dec values in radians
Parameters
----------
ra1 : array like, mandatory
dec1 :
ra2 :
dec2 :
"""
cosDelta = np.cos(dec1)* np.cos(dec2) *np.cos(ra1 -ra2) + np.sin(dec1)*np.sin(dec2)
return np.arccos(cosDelta)
In [10]:
opsimHDF = os.path.join(os.getenv('HOME'), 'data', 'LSST', 'OpSimData', 'storage.h5')
summarydf = pd.read_hdf(opsimHDF, 'table')
# df = df.query('propID == [364, 366]')
In [11]:
summarydf.query('fieldID == 290')[['fieldRA', 'fieldDec']].head()
Out[11]:
In [12]:
posns = pd.read_csv('example_data/posns.txt', delim_whitespace=True)
In [13]:
posns['ra_radians']= np.radians(posns.RAdeg)
posns['dec_radians']= np.radians(posns.Decdeg)
In [14]:
posns.head()
Out[14]:
In [15]:
posns.ix[0, 'ra_radians']
Out[15]:
In [16]:
radius = np.radians(1.75)
In [17]:
def _obsForSN(summarydf, ind):
summarydf['observed'] = 0
summarydf['distance'] = distance(summarydf['ditheredRA'],
summarydf['ditheredDec'],
#ra, dec)
posns.ix[ind, 'ra_radians'],
posns.ix[ind, 'dec_radians'])
summarydf.ix[summarydf.distance < radius, 'observed'] = 1
return summarydf.groupby('observed').get_group(1)
In [18]:
qm_0 = PerSNMetric(t0=49484, summarydf=_obsForSN(summarydf, 0), lsst_bp=lsst_bp, raCol='ditheredRA', decCol='ditheredDec')
fig = qm_0.lcplot(nightlyCoadd=True)
fig.savefig('results/SN_0.pdf')
qm_0.writeLightCurve('results/SN_0.dat', nightlyCoadd=True)
In [20]:
def writeResults(ind, mjd=49484):
print '------------------------'
print('Entering iteration ', ind)
df = _obsForSN(summarydf, ind)
print('Number of obs found ', len(df))
tmp = PerSNMetric(t0=mjd, summarydf=df,
lsst_bp=lsst_bp, raCol='ditheredRA', decCol='ditheredDec')
fig = tmp.lcplot(nightlyCoadd=True)
fig.savefig('results/SN_' + str(ind) + '_lc.pdf')
tmp.writeLightCurve('results/SN_' +str(ind) + '.dat', nightlyCoadd=True)
x = tmp.showFields()
x.savefig('results/SN_' +str(ind) + '_map.pdf')
xx = tmp.cadence_plot(summarydf=df, racol='ditheredRA', deccol='ditheredDec', mjd_center=mjd)
xx[0].savefig('results/SN_' + str(ind) + '_obs.pdf')
print ind, np.unique(tmp.fieldIds)
print '---------------------------'
In [22]:
writeResults(0)
In [23]:
writeResults(1)
In [ ]:
writeResults(1)
In [ ]:
writeResults(2)
In [ ]:
writeResults(3)
In [ ]:
writeResults(4)
In [ ]:
for i in range(len(posns)):
writeResults(i)
In [ ]:
len(posns)
In [ ]:
writeResults(8, mjd=49603)
In [ ]:
writeResults(12, mjd=49603)
In [ ]:
ss = _obsForSN(summarydf, 8)
In [ ]:
tmp = PerSNMetric(t0=49480, summarydf=ss,
lsst_bp=lsst_bp, raCol='ditheredRA', decCol='ditheredDec')
In [ ]:
xx = tmp.cadence_plot(summarydf=ss, racol='ditheredRA', deccol='ditheredDec')
In [ ]:
so = oss.SummaryOpsim(summarydf)
In [ ]:
so.mjdvalfornight(250)
In [ ]:
for ind in range(8, len(posns)):
writeResults(ind, mjd=49603)
In [ ]:
so.fi
In [ ]:
writeResults(9, 49603)
In [ ]:
so.cadence_plot(summarydf=ss)
In [ ]:
posns.to_csv('results/PosnList.csv', index=True)
In [ ]:
posns
In [ ]: