In [1]:
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
In [2]:
from lsst.sims.photUtils import Bandpass, BandpassDict
In [3]:
import lsst.sims.skybrightness as sb
from lsst.sims.photUtils import Sed, calcM5, PhotometricParameters
In [4]:
from lsst.utils import getPackageDir
In [5]:
import os
In [6]:
results = os.path.join('/Users/rbiswas/data/LSST/OpSimData/m5calcs/', 'combined_results.hdf')
df = pd.read_hdf(results)
In [7]:
df.head()
Out[7]:
In [8]:
def get_terms(i, numcols, ax=None):
row = i // numcols
col = i % numcols
if ax is None:
return row, col
else:
return ax[row, col]
In [44]:
mbins = np.arange(-3., 3., 0.1)
fig_m5, axm5 = plt.subplots(2, 3, sharex=True, sharey=True)
for i, b in enumerate('ugrizy'):
tmp = df.query('band == @b')
_x = get_terms(i, 3, axm5)
_ = _x.hist(tmp.fieldm5 - tmp.opsimm5, bins=mbins, histtype='step', lw=1, alpha=1,
color='k')
_ = _x.hist(tmp.fieldm5 - tmp.fiveSigmaDepth, bins=mbins, histtype='step', lw=0.7, alpha=0.5,
color='r')
In [45]:
mbins = np.arange(-3., 3., 0.1)
fig_m, axm = plt.subplots(2, 3, sharex=True, sharey=True)
for i, b in enumerate('ugrizy'):
tmp = df.query('band == @b')
_x = get_terms(i, 3, axm)
_ = _x.hist(tmp.skymags - tmp.filtSkyBrightness, bins=mbins, histtype='step', lw=1, alpha=1,
color='k')
axm[1, 0].set_xlabel('msky - OpSimmsky')
axm[1, 1].set_xlabel('msky - OpSimmsky')
axm[1, 2].set_xlabel('msky - OpSimmsky')
axm[0, 0].set_ylabel('PDF')
axm[1, 0].set_ylabel('PDF')
Out[45]:
In [ ]:
fig_m.savefig('msky.pdf')
def atmTransName(airmass): """ return filename for atmospheric transmission with aerosols for airmass closest to input
Parameters
----------
airmass : airmass
"""
l = np.arange(1.0, 2.51, 0.1)
idx = np.abs(l - airmass).argmin()
a = np.int(10*l[idx])
baseline = getPackageDir('THROUGHPUTS')
fname = os.path.join(baseline, 'atmos', 'atmos_{}_aerosol.dat'.format(a))
return fname
def mCalcs(airmass, bandName, ra, dec, expMJD, FWHMeff, hwbpdict, photparams=None, sm=None): """ sm : """ if photparams is None: photparams = PhotometricParameters() if sm is None: sm = SkyModel(observatory='LSST', mags=False, preciseAltAz=True)
# Obtain full sky transmission at airmass
# Note that this method is not interpolating but choosing the atmospheric transmission from
# Modtran simulations of the closest airmass in a sequence of np.arange(1., 2.51, 0.1)
fname = atmTransName(airmass)
print(fname)
atmTrans = np.loadtxt(fname)
wave, trans = hwbpdict[bandName].multiplyThroughputs(atmTrans[:, 0], atmTrans[:, 1])
bp = Bandpass(wavelen=wave, sb=trans)
# Set the observing condition
sm.setRaDecMjd(lon=[ra], lat=[dec], filterNames=[bandName],
mjd=expMJD, degrees=False, azAlt=False)
# Get the sky sed
wave, spec = sm.returnWaveSpec()
sed = Sed(wavelen=wave, flambda=spec[0])
m5 = calcM5(sed, bp, hwbpdict[bandName], photparams, FWHMeff)
# Get the sky magnitude only in the band concerned
m = sm.returnMags(bandpasses=hwbpdict)[bandName][0]
return m5, m
In [10]:
from brightness import mCalcs, atmTransName
In [11]:
df.head()
Out[11]:
I copied the values of expMJD, FWHMeff, fieldRA, fieldDec from a sqlite query by hand (will change in a later iteration). Comparison of m5 and msky values from the dataframe for the first two rows is quite good.
In [12]:
tot, hwbpdict = BandpassDict.loadBandpassesFromFiles()
In [13]:
from lsst.sims.skybrightness import SkyModel
In [14]:
mCalcs(1.464, 'y', 1.676483, -1.082473, 59580.033829, 1.263038, hwbpdict)
Out[14]:
In [15]:
mCalcs(1.454958, 'y', 1.69412, -1.033972, 59580.034275, 1.258561, hwbpdict)
Out[15]:
In [16]:
df['diffmsky'] = df['skymags'] - df['filtSkyBrightness']
In [17]:
df.query('diffmsky > 3.3')
Out[17]:
In [58]:
mCalcs(1.008652, 'g', 0.925184, -0.4789, 61044.077855, 1.086662, hwbpdict=hwbpdict)
Out[58]:
In [18]:
from opsimsummary import OpSimOutput
In [19]:
opsout = OpSimOutput.fromOpSimDB('/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db')
In [20]:
row = opsout.summary.ix[994149]
In [21]:
mCalcs(row.airmass, row['filter'], row.fieldRA, row.fieldDec, row.expMJD, row.FWHMeff, hwbpdict=hwbpdict)
Out[21]:
In [26]:
sb.__version__
Out[26]:
In [27]:
np.__version__
Out[27]:
In [ ]:
Sed.writeSED()
In [28]:
laptop = np.loadtxt('skySED_laptop.csv')
In [29]:
lsstuw = np.loadtxt('skySED.csv')
In [38]:
laptop.shape
Out[38]:
In [40]:
import healpy as hp
In [42]:
hp.__file__
Out[42]:
In [39]:
fig, ax = plt.subplots(1, 2)
ax[0].plot(laptop[:, 0], laptop[:, 1] - lsstuw[:,1], 'k-')
#ax[0].plot(lsstuw[:, 0], lsstuw[:, 1])
ax[1].plot(laptop[:, 0], laptop[:, 1]/lsstuw[:, 1])
Out[39]:
In [ ]: