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


/usr/local/manual/anaconda/lib/python2.7/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.
  "You should import from ipykernel or jupyter_client instead.", ShimWarning)

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)


u
g
r
i
z
y

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]:
fieldRA fieldDec
95688 6.097944 -1.10516
95703 6.097944 -1.10516
95704 6.097944 -1.10516
96854 6.097944 -1.10516
96863 6.097944 -1.10516

List of Positions to Try out


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]:
position RAdeg Decdeg datapts ra_radians dec_radians
0 lsstposnew1 280.000 -40.000 86 4.886922 -0.698132
1 lsstposnnew2 100.000 -20.000 53 1.745329 -0.349066
2 deep1 349.386 -63.321 2363 6.097936 -1.105160
3 lsstpos3 190.000 -83.000 23942 3.316126 -1.448623
4 lsstpos4 20.000 -83.000 25244 0.349066 -1.448623

Obtaining the Light Curves


In [15]:
posns.ix[0, 'ra_radians']


Out[15]:
4.8869219055841224

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)


metrics.py:151: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
  df.sort('SNR', ascending=False, inplace=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)


------------------------
('Entering iteration ', 0)
('Number of obs found ', 1112)
0 [ 846  868  938  948  968 1053]
---------------------------

In [23]:
writeResults(1)


------------------------
('Entering iteration ', 1)
('Number of obs found ', 1157)
1 [1629 1720 1752 1812 1845]
---------------------------

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 [ ]: