In [1]:
import time
import numpy as np
import pandas as pd
In [2]:
import lsst.sims.skybrightness as sb
from lsst.sims.utils import Site
import ephem
from obscond import SkyCalculations as sm
In [3]:
from opsimsummary import SynOpSim
In [4]:
from obscond.observingPotential import ObservationPotential
In [5]:
from lsst.sims.utils import angularSeparation
In [6]:
from scipy.interpolate import interp1d
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
sns.set_context('notebook')
In [8]:
import obscond
In [9]:
print(obscond.__version__)
In [10]:
# uDDF center Design DC2 Run1 Note
rauDDF = np.radians(53.125)
decuDDF = np.radians(-28.100)
In [11]:
op = ObservationPotential(rauDDF, decuDDF)
In [12]:
# 60 sec intervals of mjd values between 7 and 10 years
t0 = 59580 + 7 * 365
sec = 1.0 / 24.0 /60. /60.
hr = 60. * 60. * sec
t2 = np.arange(t0, t0 + 365 *3, 60 * sec)
In [13]:
opdf = op.potential_obscond(t2, fieldRA=rauDDF, fieldDec=decuDDF)
In [14]:
len(opdf)/ np.pi /1.0e7
Out[14]:
In [15]:
opdf.head()
Out[15]:
In [16]:
# new
constraints = 'alt > 30 and alt < 85 and sunAlt < -12 and moonDist > 30'
availabletimes = op.available_times(opdf, constraints=constraints)
In [17]:
#new
len(availabletimes)
Out[17]:
In [18]:
# new
fig, ax = plt.subplots()
#ax0, ax1 = ax
#dff = opdf.query('alt < 85. and alt > 20.')
ax.plot(opdf.mjd, opdf.sunAlt, label='sun Alt')
ax.plot(availabletimes.mjd, availabletimes.sunAlt, label='field constraints')
ax.axhline(-12.0, ls='dashed', lw=2)
ax.set_xlabel('MJD (day)')
ax.set_ylabel('sun Altitude (deg)')
fig.savefig('availabletimes_block2.png')
In [19]:
synopsim = SynOpSim.fromOpSimDB('/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db',
subset='ddf', opsimversion='lsstv3')
In [20]:
all_baseline_visits = synopsim.pointings.query('fieldID == 1427')
In [21]:
baseline_block2 = synopsim.pointings.query('fieldID == 1427 and night > 365*7')
In [22]:
unights = baseline_block2.query('filter == "u"').night.unique()
In [23]:
baseline_nights_2 = baseline_block2.query('filter != "u"').night.unique()
In [24]:
print('Number of nights using grizy filters is {}'.format(baseline_nights_2.size))
In [25]:
ratio_sampling = 0.5 # See DC2 design note (0.75)
required_nights = np.ceil(baseline_nights_2.size / ratio_sampling).astype(np.int)
In [26]:
print('Given the requirement of the cadence being {0}, the number of nights is {1}'.format(ratio_sampling, required_nights))
In [27]:
times = availabletimes.query('night in @baseline_nights_2').groupby('night').agg(dict(mjd=ObservationPotential.timerange)).rename(columns=dict(mjd='time'))
In [28]:
fig, ax = plt.subplots()
_ = times.time.hist(bins=np.arange(2., 10, 0.25), histtype='step',
lw=2, ax=ax)
ax.axvline(5.0, color='k', ls='dashed')
ax.set_ylabel('nights')
ax.set_xlabel('hours of night available per night')
ax.set_title('Hours/night for minion DDF nights in yrs 4-7')
fig.savefig('block_2_hours_pernight.png')
In [29]:
len(availabletimes)
Out[29]:
In [37]:
synopsim.pointings.index.values.max()
Out[37]:
In [29]:
shortnight = availabletimes.groupby('night').agg(dict(mjd=ObservationPotential.timerange))
shortnights = shortnight.query('mjd < 5').index.unique()
print(len(shortnights))
In [30]:
unights = baseline_block2.query('filter == "u"').night.unique()
In [31]:
target_nights = availabletimes.query('night not in @shortnights').query('night not in @unights').night.unique()
In [32]:
rng = np.random.RandomState(1)
target_nights.sort()
chosen_nights = rng.choice(target_nights, replace=False, size=required_nights)
chosen_nights.sort()
In [33]:
nstats = op.nightStats(availabletimes)
st_times = op.start_times(nstats, chosen_nights, rng)
In [34]:
dc2_visits_block2 = op.dc2_visits(st_times, year_block=2, pointings=all_baseline_visits)
In [35]:
dc2_visits_block2.groupby('filter').agg(dict(expMJD='count'))
Out[35]:
In [37]:
fig, ax = plt.subplots()
dc2_visits_block2.rawSeeing.hist(ax=ax,**dict(histtype='step', normed=1))
all_baseline_visits.rawSeeing.hist(ax=ax, **dict(histtype='step', normed=1))
Out[37]:
In [39]:
dc2_visits_block2['obsHistID'] = np.arange(len(dc2_visits_block2)) + 2500000 + 20000
In [63]:
dc2_visits_block2.to_csv('dc2_blocks2_v2.csv', index=False)
In [41]:
# 60 sec intervals of mjd values between 7 and 10 years
t0 = 59580 + 4 * 365
sec = 1.0 / 24.0 /60. /60.
hr = 60. * 60. * sec
t01 = np.arange(t0, t0 + 365 *3, 60 * sec)
In [42]:
opdf = op.potential_obscond(t01, fieldRA=rauDDF, fieldDec=decuDDF)
In [43]:
# new
constraints_1 = 'alt > 30 and alt < 85 and sunAlt < -12 and moonDist > 30'
availabletimes_1 = op.available_times(opdf, constraints=constraints)
In [44]:
(availabletimes_1.mjd.min() - 59580)/365
Out[44]:
In [45]:
shortnight_1 = availabletimes_1.groupby('night').agg(dict(mjd=ObservationPotential.timerange))
shortnights_1 = shortnight_1.query('mjd < 5').index.unique()
print(len(shortnights_1))
In [46]:
baseline_block1 = synopsim.pointings.query('fieldID == 1427 and night < 365*7 and night > 365*4')
In [47]:
unights_1 = baseline_block1.query('filter == "u"').night.unique()
In [48]:
target_nights_1 = availabletimes_1.query('night not in @shortnights_1').query('night not in @unights_1').night.unique()
In [49]:
baseline_nights_1 = baseline_block1.query('filter != "u"').night.unique()
In [50]:
ratio_sampling_1 = 0.75 # See DC2 design note (0.75)
required_nights_1 = np.ceil(baseline_nights_1.size / ratio_sampling).astype(np.int)
In [51]:
required_nights_1
Out[51]:
In [52]:
rng = np.random.RandomState(1)
target_nights_1.sort()
chosen_nights_1 = rng.choice(target_nights_1, replace=False, size=required_nights_1)
chosen_nights_1.sort()
In [53]:
nstats_1 = op.nightStats(availabletimes_1)
st_times_1 = op.start_times(nstats_1, chosen_nights_1, rng)
In [54]:
dc2_visits_block1 = op.dc2_visits(st_times_1, year_block=1, pointings=baseline_block1)
In [69]:
dc2_visits_block1.groupby('filter').agg(dict(expMJD='count'))
Out[69]:
In [58]:
dc2_visits_block1.columns
Out[58]:
In [60]:
dc2_visits_block1['obsHistID'] = np.arange(len(dc2_visits_block1)) + 2500000
In [64]:
dc2_visits_block1.to_csv('dc2_blocks1_v2.csv', index=False)
In [59]:
len(dc2_visits_block1)
Out[59]:
In [60]:
len(dc2_visits_block2)
Out[60]:
In [61]:
dc2_visits_block1.night = dc2_visits_block1.night.astype(int)
In [62]:
np.bincount(baseline_block1.sort_values(by='night').night.diff()[1:].astype(np.int))
Out[62]:
In [107]:
np.bincount(dc2_visits_block1.sort_values(by='night').night.diff()[1:].astype(np.int))
Out[107]:
In [71]:
!head dc2_block1.csv
In [ ]:
In [ ]:
In [182]:
dc2_visits_block1 = op.dc2_visits(st_times.iloc[0], 1)
In [ ]:
In [140]:
st_times.iloc[0]
Out[140]:
In [68]:
dc2_visits_block2.to_csv('dc2_visits_block2.csv')
In [138]:
dc2_visits_block2[0].sort_values(by='expMJD').expMJD
Out[138]:
In [41]:
dc2_visits_block2 = ObservationPotential.dc2_visits(st_times, year_block=2, delta=1)
In [42]:
dc2_visits_block2.head()
Out[42]:
In [48]:
dc2_visits_block2.sort_values(by='expMJD').expMJD.diff().min()*24*60*60
Out[48]:
In [50]:
dc2_visits_block2.groupby('filter').agg(dict(night='count'))
Out[50]:
In [62]:
minion_block2 = synopsim.pointings.query('night > 365 *7 and filter != "u" and fieldID == 1427')
In [63]:
minion_block2['filter'].unique().size
Out[63]:
In [67]:
dc2_visits_block2.groupby(['night', 'filter']).agg(dict(expMJD='count'))
Out[67]:
In [59]:
dc2_visits_block2.night.unique().size
Out[59]:
In [64]:
minion_block2.groupby('filter').agg(dict(night='count'))
Out[64]:
In [65]:
len(minion_block2)
Out[65]:
In [66]:
len(dc2_visits_block2)
Out[66]:
In [120]:
dc2_visits_block2.sort_values(by='expMJD').expMJD
In [99]:
vals = dc2_visits_block2.sort_values(by='expMJD').expMJD.diff()*24.0*60.* 60.
In [82]:
st_times
Out[82]:
In [ ]:
op.nightStats
In [31]:
availabletimes.night.unique().size
Out[31]:
In [32]:
availabletimes.query('night not in @shortnights').night.unique().size
Out[32]:
In [75]:
len(target_nights)
Out[75]:
In [77]:
expected_cadence = target_nights.size / np.float(required_nights)
print(expected_cadence)
In [37]:
targets.night.unique().size
Out[37]:
In [38]:
rng = np.random.RandomState(1)
target_nights.sort()
chosen_nights = rng.choice(target_nights, replace=False, size=required_nights)
chosen_nights.sort()
In [39]:
chosen_nights.size
Out[39]:
In [ ]:
chosen_nights
In [ ]:
st_times = op.start_times(nigh)
In [40]:
baseline_nights
Out[40]:
In [41]:
shortnights[:40]
Out[41]:
In [42]:
target_nights
Out[42]:
In [43]:
unights
Out[43]:
In [46]:
_ = plt.hist(np.diff(chosen_nights), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='designed')
_ = plt.hist(np.diff(np.sort(baseline_nights)), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='minion')
In [60]:
np.bincount(np.sort(np.diff(baseline_nights)))
Out[60]:
In [59]:
np.bincount(np.sort(np.diff(chosen_nights)))
Out[59]:
In [44]:
fig, ax = plt.subplots(1, 2)
_ = ax[0].hist(np.diff(chosen_nights), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='designed')
_ = ax[0].hist(np.diff(np.sort(baseline_nights)), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='minion')
_ = ax[1].hist(np.diff(chosen_nights), bins=np.arange(0, 80, 1),
histtype='step', lw=2, cumulative=1, label='designed', normed=0)
_ = ax[1].hist(np.diff(np.sort(baseline_nights)), bins=np.arange(0, 80, 1),
histtype='step', lw=2, cumulative=1, label='minion', normed=0)
ax[1].set_xlim(xmax=10)
plt.legend(loc='best')
fig.savefig('gap_block2.png')
In [ ]:
In [78]:
fig, ax = plt.subplots(1, 2)
_ = ax[0].hist(np.diff(chosen_nights), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='designed')
_ = ax[0].hist(np.diff(np.sort(baseline_nights)), bins=np.arange(0, 80, 1),
histtype='step', lw=2, normed=0, label='minion')
_ = ax[1].hist(np.diff(chosen_nights), bins=np.arange(0, 80, 1),
histtype='step', lw=2, cumulative=1, label='designed', normed=0)
_ = ax[1].hist(np.diff(np.sort(baseline_nights)), bins=np.arange(0, 80, 1),
histtype='step', lw=2, cumulative=1, label='minion', normed=0)
ax[1].set_xlim(xmax=10)
plt.legend(loc='best')
fig.savefig('gap_block2.png')
In [50]:
fig, ax = plt.subplots()
ax.plot(chosen_nights[:-1], np.diff(chosen_nights), '.')
for nt in unights:
ax.axvline(nt, ls='dashed', color='k', lw=1)
#for nt in shortnights:
# ax.axvline(nt, ls='dashed', color='r', lw=0.1)
In [62]:
availabletimes.describe()
Out[62]:
In [63]:
nstats = op.nightStats(availabletimes)
In [69]:
nstats.head()
Out[69]:
In [72]:
nstats.loc[chosen_nights].head()
Out[72]:
In [66]:
availabletimes.groupby('night').agg(dict(mjd='count')).describe()
Out[66]:
In [56]:
availabletimes.groupby('night').agg(dict(mjd='count')).apply(lambda x: x/60.).describe()
Out[56]:
In [57]:
nstats.loc[chosen_nights]
Out[57]:
In [73]:
st_times = op.start_times(nstats, chosen_nights, rng)
In [74]:
st_times
Out[74]:
In [63]:
v, vl = ObservationPotential.dc2_sequence(st_times.values[0], year_block=2)
In [64]:
dc2_visits = ObservationPotential.dc2_visits(st_times, year_block=1, delta=1)
In [65]:
dc2_visits.head()
Out[65]:
In [66]:
baseline_block2.expMJD.size
Out[66]:
In [67]:
(dc2_visits.expMJD.size + baseline_block2.query('filter == "u"').expMJD.size)
Out[67]:
In [68]:
dc2_visits['fieldRA'] = rauDDF
dc2_visits['fieldDec'] = decuDDF
In [69]:
from lsst.sims.photUtils import BandpassDict
import obscond
In [70]:
dc2_visits['FWHMeff'] = 1.0
In [71]:
TotbpDict, hwbpDict = BandpassDict.loadBandpassesFromFiles()
skycalc = obscond.SkyCalculations(photparams="LSST", hwBandpassDict=hwbpDict)
In [72]:
from lsst.sims.skybrightness import SkyModel
In [73]:
synopsim.pointings.iloc[0][['fieldRA', 'fieldDec', 'expMJD']]
Out[73]:
In [159]:
sm = SkyModel()
In [163]:
decuDDF
Out[163]:
In [164]:
rauDDF
Out[164]:
In [175]:
dc2_visits.iloc[1]
Out[175]:
In [ ]:
In [84]:
#sm.setRaDecMjd(lat=-0.4904, lon=0.927, mjd=61137.5)
In [74]:
dc2_visits.to_csv('dc2_visits_block2')
In [181]:
sm.returnMags(bandpasses=hwbpDict)
Out[181]:
In [208]:
dc2_visits.iloc[293]
Out[208]:
In [102]:
dc2_visits.to_csv('dc2_vists_block1.csv')
In [85]:
skycalc.calculatePointings(dc2_visits.head(293), calcDepths=False)
In [132]:
from obscond.historicalWeatherData import SeeingFile
In [143]:
seeing = synopsim.pointings[['expMJD', 'rawSeeing']].sort_values(by='expMJD')
In [142]:
from scipy.interpolate import interp1d
In [144]:
rawSeeing = interp1d(seeing.expMJD, seeing.rawSeeing)
In [149]:
seeing.iloc[5:10]['rawSeeing']
Out[149]:
In [145]:
rawSeeing(seeing.expMJD.iloc[5:10])
Out[145]:
In [140]:
dc2_visits.expMJD.values
Out[140]:
In [116]:
dc2_visits['night'] = np.floor(dc2_visits.expMJD - 59580)
In [117]:
dc2_visits
Out[117]:
In [282]:
df = availabletimes.groupby('night').agg(dict(mjd=(min, max, timerange))).reset_index()
In [285]:
df.values.shape
Out[285]:
In [286]:
df = pd.DataFrame(df.values, columns=('night', 'minmjd', 'maxmjd', 'timerange'))
In [287]:
df.timediff = (df.maxmjd - df.minmjd)* 24
In [289]:
all(df.timerange == df.timediff)
Out[289]:
In [ ]:
In [88]:
def mintime(series):
return np.remainder(min(series), 1.0) * 24.0
In [ ]:
x = dff.groupby('night').aggregate(dict(alt='count')).reset_index()
In [ ]:
In [ ]:
1480 +
In [102]:
np.degrees(op.sunAlt_singleTime(61040.0))
Out[102]:
In [103]:
op.obs
Out[103]:
In [104]:
op.doff
Out[104]:
In [107]:
s = ephem.Sun()
In [109]:
op.obs.date = 61040. - op.doff
In [ ]:
sun.compute(lsstObs)
return sun.alt
In [110]:
s.compute(op.obs)
In [111]:
s.alt
Out[111]:
In [140]:
op.sunAlt_singleTime(61040.)
Out[140]:
In [113]:
from lsst.sims.utils import Site
In [114]:
site = Site('LSST')
In [115]:
import ephem
In [116]:
sun = ephem.Sun()
In [117]:
lsstObs = ephem.Observer()
lsstObs.lat = np.radians(site.latitude)
lsstObs.lon = np.radians(site.longitude)
lsstObs.elevation = site.height
In [118]:
lsstObs.date = 61040 - doff
In [120]:
sun.compute(lsstObs)
In [122]:
sun.alt
Out[122]:
In [124]:
lsstObs.lat
Out[124]:
In [125]:
op.obs.lat
Out[125]:
In [126]:
lsstObs.long
Out[126]:
In [127]:
lsstObs.lat
Out[127]:
In [265]:
1480 + 59580
Out[265]:
In [266]:
opdf.query('night > 1479 and night < 1520')
Out[266]:
In [211]:
synopsim.pointings.airmass.hist()
Out[211]:
In [ ]: