In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
from sqlalchemy import create_engine
import numpy as np
import pandas as pd
import seaborn as sb
import os
import healpy as hp
In [2]:
import lsst.sims.maf.db as db
In [3]:
from sqlalchemy import create_engine
In [337]:
engine = create_engine('sqlite:////Users/rbiswas/data/LSST/OpSimData/enigma_1189_sqlite.db')
In [338]:
Proposal = pd.read_sql_table('Proposal', engine)
In [339]:
Proposal
Out[339]:
In [324]:
opsdb = db.OpsimDatabase('sqlite:////Users/rbiswas/data/LSST/OpSimData/enigma_1189_sqlite.db')
In [325]:
propID, propTags = opsdb.fetchPropInfo()
In [328]:
propID
Out[328]:
In [4]:
Summary = pd.read_hdf('storage.h5', 'table')
In [5]:
Summary.info()
In [8]:
s = Summary[Summary['airmass'] <1.2]
In [ ]:
In [ ]:
plt.plot(Summary['expDate'], Summary['expMJD'])
In [331]:
fulltablesize = Summary.size
In [332]:
dd = Summary[Summary['propID']==366]
ddsize = dd.size
wfd = Summary[Summary['propID']==364]
wfdsize = wfd.size
In [335]:
from __future__ import division
ddsize/ fulltablesize , wfdsize / fulltablesize
Out[335]:
In [344]:
fseason = Summary[(Summary['night'] < 366) & (Summary['propID']==364)]
In [345]:
testdf = fseason.copy(deep=True)
In [346]:
add_simlibCols(testdf)
Out[346]:
In [48]:
fieldsimlibs = testdf.groupby(by='fieldID')
In [258]:
def dec(fieldID):
decvals = np.unique(fieldsimlibs.get_group(fieldID).fieldDec.values)
if len(decvals)==1:
return decvals[0]
else:
raise ValueError('The fieldDec of this group seems to not be unique\n')
In [319]:
class SummaryOpsim(object):
import os, subprocess
def __init__(self, summarydf, user=None, host=None, survey='LSST', telescope='LSST', pixSize=0.2):
self.df = summarydf.copy(deep=True)
self.fieldsimlibs = self.df.groupby(by='fieldID')
self.fieldIds = self.fieldsimlibs.groups.keys()
if user is None:
user = os.getlogin()
self.user = user
if host is None:
proc = subprocess.Popen('hostname', stdout=subprocess.PIPE)
host, err = proc.communicate()
self.host = host
self.survey = survey
self.host = host
self.telescope = telescope
def ra(self, fieldID):
ravals = np.unique(self.fieldsimlibs.get_group(fieldID).fieldRA.values)
if len(ravals)==1:
return ravals[0]
else:
raise ValueError('The fieldDec of this group seems to not be unique\n')
def dec(self, fieldID):
decvals = np.unique(self.fieldsimlibs.get_group(fieldID).fieldDec.values)
if len(decvals)==1:
return decvals[0]
else:
raise ValueError('The fieldDec of this group seems to not be unique\n')
def meta(self, fieldID):
meta = {}
meta['LIBID'] = fieldID
meta['dec'] = self.dec(fieldID)
return meta
def fieldheader(self, fieldID):
ra = self.ra(fieldID)
dec = self.dec(fieldID)
nobs = self.df.size
s = '# --------------------------------------------' +'\n'
s += 'LIBID: {0:10d}'.format(fieldID) +'\n'
s += 'RA: {0:+10.6f} DECL: {1:+10.6f} NOBS: {2:10d} MWEBV: {3:5.2f} PIXSIZE: {4:5.3f}'.format(ra, dec, nobs, 0.0, 0.2) + '\n'
s += 'LIBID: {0:10d}'.format(fieldID) + '\n'
s += '# CCD CCD PSF1 PSF2 PSF2/1' +'\n'
s += '# MJD IDEXPT FLT GAIN NOISE SKYSIG (pixels) RATIO ZPTAVG ZPTERR MAG' + '\n'
return s
def fieldfooter(self, fieldID):
s = 'END_LIBID: {0:10d}'.format(fieldID)
s += '\n'
return s
def formatSimLibField(self, fieldID, sep=' '):
opSimSummary = self.fieldsimlibs.get_group(fieldID)
y =''
for row in opSimSummary.iterrows():
data = row[1] # skip the index
# MJD EXPID FILTER
lst = ['S:', "{0:5.3f}".format(data.expMJD), "{0:10d}".format(data.obsHistID), data['filter'],
# CCD Gain CCD Noise
"{0:5.2f}".format(1.), "{0:5.2f}".format(0.25),
# SKYSIG
"{0:6.2f}".format(data.simLibSkySig),
# PSF1 PSF2 PSFRatio
"{0:4.2f}".format(data.simLibPsf), "{0:4.2f}".format(0.), "{0:4.3f}".format(0.),
# ZPTAVG MAG
"{0:6.2f}".format(data.simLibZPTAVG), "{0:6.3f}".format(0.005), "{0:+7.3f}".format(-99.) ]
s = sep.join(lst)
y += s + '\n'
return y
def writeSimLibField(self, fieldID):
s = self.fieldheader(fieldID)
s += self.formatSimLibField(fieldID, sep=' ')
s += self.footer(fieldID)
return s
def simLibheader(self): #, user=None, host=None, survey='LSST', telescope='LSST'):
s = 'SURVEY: {0:} FILTERS: ugrizY TELESCOPE: {1:}\n'.format(self.survey, self.telescope)
s += 'USER: {0:} HOST: {1:}'.format(self.user, self.host)
s +=
return s
def writeSimlib(self, filename, comments=''):
with open(filename, 'w') as fh:
simlib_header = self.simLibheader()
simlib_footer = self.simLibFooter()
fh.write(simlib_header)
fh.write(comments)
fh.write('BEGIN LIBGEN\n')
fh.write('\n')
for fieldID in self.fieldIds:
fh.write(self.fieldheader(fieldID))
fh.write(self.formatSimLibField(fieldID))
fh.write(self.fieldfooter(fieldID))
fh.write(simlib_footer)
In [347]:
SO = SummaryOpsim(testdf)
In [348]:
SO.writeSimlib('fullsimlib.out')
In [359]:
Summary.replace({'filter': {'y','Y'}})
Out[359]:
In [302]:
ss = SO.simLibheader()
In [135]:
xx = df[['expMJD', 'obsHistID','filter' ,'simLibSkySig', 'simLibPsf', 'simLibZPTAVG']].head()
In [136]:
xx
Out[136]:
In [107]:
xx.to_string(header=False, index=False, float_format="{0:.4f}".format())
In [160]:
xx['filter'] = xx['filter'].astype(str)
In [172]:
y = np.asarray(xx['filter'],dtype='S1')
In [173]:
Out[173]:
In [250]:
def formatSimLibField(opSimSummary, sep=' '):
y =''
for row in opSimSummary.iterrows():
data = row[1] # skip the index
# MJD EXPID FILTER
lst = ['S:', "{0:5.3f}".format(data.expMJD), "{0:10d}".format(data.obsHistID), data['filter'],
# CCD Gain CCD Noise
"{0:5.2f}".format(1.), "{0:5.2f}".format(0.25),
# SKYSIG
"{0:6.2f}".format(data.simLibSkySig),
# PSF1 PSF2 PSFRatio
"{0:4.2f}".format(data.simLibPsf), "{0:4.2f}".format(0.), "{0:4.3f}".format(0.),
# ZPTAVG MAG
"{0:6.2f}".format(data.simLibZPTAVG), "{0:6.3f}".format(0.005), "{0:+7.3f}".format(-99.) ]
s = sep.join(lst)
y += s + '\n'
return y
In [303]:
f = open('test.simlib', 'w')
f.write(ss)
f.close()
In [244]:
"{0:>6.3f}".format(1)
Out[244]:
In [235]:
ss
Out[235]:
In [210]:
li = [12,45,78,784,2,69,1254,4785,984]
print map('the number is {}'.format,li)
In [217]:
'{:0.2f}'.format(1.)
Out[217]:
In [182]:
lst = ['S:', '49353.176']
' '.join(lst)
Out[182]:
In [185]:
lst = ['S:', "{0:5.3f}".format(49353.176)]
' '.join(lst)
Out[185]:
In [68]:
Summary.ditheredRA.describe()
Out[68]:
In [53]:
df = fieldsimlibs.get_group(fieldIds[0])
In [78]:
-6.1 % (2*np.pi)
Out[78]:
In [79]:
plt.plot((Summary['ditheredRA'] - Summary['fieldRA']).values % (2. *np.pi), (Summary['ditheredDec'] - Summary['fieldDec']).values,'o')
Out[79]:
In [75]:
Summary[Summary['ditheredRA'] - Summary['fieldRA'] < -2.][['ditheredRA', 'fieldRA']]
Out[75]:
In [40]:
def add_simlibCols(opsimtable, pixSize=0.2):
opsim_seeing = opsimtable['finSeeing']
opsim_maglim = opsimtable['fiveSigmaDepth']
opsim_magsky = opsimtable['filtSkyBrightness']
# area = (opsim_seeing, pixsize=0.2):
opsimtable['simLibPsf'] = opsim_seeing /2.35 /pixSize
area = (1.51 * opsim_seeing)**2.
opsim_snr = 5.
arg = area * opsim_snr * opsim_snr
zpt_approx = 2.0 * opsim_maglim - opsim_magsky + 2.5 * np.log10(arg)
# ARG again in David Cinabro's code
val = -0.4 * (opsim_magsky - opsim_maglim)
tmp = 10.0 ** val
zpt_cor = 2.5 * np.log10(1.0 + 1.0 / (area * tmp))
simlib_zptavg = zpt_approx + zpt_cor
# ZERO PT CALCULATION
opsimtable['simLibZPTAVG'] = simlib_zptavg
#SKYSIG Calculation
npix_asec = 1./ pixSize**2.
opsimtable['simLibSkySig'] = ((1.0/ npix_asec)*10.0 **(-0.4 * (opsim_maglim - simlib_zptavg)))**0.5
return opsimtable
In [14]:
s = add_simlibCols(Summary)
In [19]:
groups = s.groupby(by='fieldID')
In [24]:
groups.get_group(1).night.max()
Out[24]:
In [30]:
In [ ]:
s.values
In [35]:
l = map(lambda x: len(np.unique(groups.get_group(x).fieldRA.values)), groups.groups.keys())
In [5]:
t = Summary.loc[Summary['obsHistID'].isin(range(5))]
In [36]:
type(l)
Out[36]:
In [39]:
np.unique(np.array(l))
Out[39]:
In [6]:
t
Out[6]:
In [7]:
type(t)
Out[7]:
In [10]:
y = add_simlibCols(t)
In [12]:
y[['simLibPsf', 'skysig', 'ZPTAVG']]
Out[12]:
In [ ]:
Summary.loc[Summary[]]
In [75]:
fig, ax = plt.subplots(2,3)
In [ ]:
In [76]:
Summary.fiveSigmaDepth.hist(by=Summary['filter'], ax=ax, histtype='step', color='b', lw=2.0)
Out[76]:
In [78]:
fig
Out[78]:
In [77]:
Summary.filtSkyBrightness.hist(by=Summary['filter'], ax=ax, color='r', lw=2.0, histtype='step')
Out[77]:
In [68]:
t = Summary.loc[Summary['night'].isin(range(365))].head(100)
In [74]:
np.unique(Summary.loc[Summary['night'].isin(range(365))]['filter'].values)
Out[74]:
In [65]:
t.plot('expMJD', 'filtSkyBrightness')
Out[65]:
In [80]: