This notebook shows some quantities for the design of a new OpSim run in multiple Declination bands. The idea is to run only the WFD for the last nine years of LSST. We want to obtain

3 equal area declination bands covering the area surveyed by WFD
Obtain the FieldIDs corresponding to these different bands
Get an estimate of the time spent in WFD vs all proposals in the baseline cadence



In [1]:

    
import numpy as np



In [2]:

    
from opsimsummary import OpSimOutput



In [3]:

    
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt



In [4]:

    
opsout = OpSimOutput.fromOpSimDB('/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db',
                                subset='_all')









    



 reading from database sqlite:////Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db



In [5]:

    
opsout.summary.expMJD.size









    Out[5]:





2448282



In [6]:

    
opsout.summary.reset_index().drop_duplicates(subset='obsHistID').expMJD.size









    Out[6]:





2447931

Declination Boundaries for the WFD proposal



In [7]:

    
maxdec = opsout.summary.query('propID == 54').fieldDec.max()
mindec = opsout.summary.query('propID == 54').fieldDec.min()
maxtheta = -  mindec + np.pi /2.0
mintheta = -  maxdec + np.pi / 2.0

Total area covered, should be ~ 0.5 sky



In [8]:

    
area = 2 * np.pi * (np.cos(mintheta) - np.cos(maxtheta))



In [9]:

    
# Check that the sky fraction is ~ 0.5
area / 4 / np.pi









    Out[9]:





0.46500948237707568

Equations for the boundary values:

We want to divide the area covered by the by the WFD survey by two colatitudes into 3 bands of equal area. Therefore the equations are :

2 np.pi (np.cos(mintheta) - cosb1) == area / 3.

2 np.pi (cosb2 - np.cos(maxtheta)) == area / 3.



In [10]:

    
b1 = np.arccos((2 * np.pi * np.cos(mintheta) - area /3)/ 2 /np.pi)
b1dec = - b1 + np.pi/2



In [11]:

    
np.degrees(b1dec)









    Out[11]:





-15.254591573425083



In [12]:

    
b2 = np.arccos(( area /3 + 2 * np.pi * np.cos(maxtheta) )/ 2 /np.pi)
b2dec = -b2 + np.pi/2.



In [13]:

    
np.degrees(b2dec)









    Out[13]:





-34.967720698496933



In [14]:

    
def numFieldsinPatch(mindecval, maxdecval, patchName=None, propID=54):
    #print mindecval, maxdecval#, patchname
    if patchName is None:
        patchName = ''
    query_str ='fieldDec > {0} and fieldDec < {1}'.format(mindecval, maxdecval)
    print (query_str)
    dfpatch = opsout.summary.query(query_str)
    numFields = len(dfpatch.reset_index().drop_duplicates(subset='obsHistID').fieldID.unique())
    numFields_WFD = len(dfpatch.query('propID == {}'.format(propID)).fieldID.unique())
    print('Total Number of fields in {2} patch is {0} and {1} are in WFD'.format(numFields,
                                                                                 numFields_WFD,
                                                                                 patchName))
    return dfpatch, numFields, numFields_WFD



In [15]:

    
dfNorthPatch, numFieldsNorth, numFieldsNorthWFD = numFieldsinPatch('@b1dec', '@maxdec', patchName='North')









    



fieldDec > @b1dec and fieldDec < @maxdec
Total Number of fields in North patch is 832 and 795 are in WFD



In [16]:

    
dfMiddlePatch, numFieldsMiddle, numFieldsMiddleWFD = numFieldsinPatch('@b2dec', 
                                                                   '@b1dec', patchName='Middle')









    



fieldDec > @b2dec and fieldDec < @b1dec
Total Number of fields in Middle patch is 820 and 767 are in WFD



In [17]:

    
dfSouthPatch, numFieldsSouth, numFieldsSouthWFD = numFieldsinPatch('@mindec', 
                                                                   '@b2dec', patchName='South')









    



fieldDec > @mindec and fieldDec < @b2dec
Total Number of fields in South patch is 806 and 722 are in WFD



In [18]:

    
NorthFieldIDs = dfNorthPatch.query('propID == 54').fieldID
MiddleFielIDs = dfMiddlePatch.query('propID == 54').fieldID
SouthFielIDs = dfSouthPatch.query('propID == 54').fieldID



In [19]:

    
NorthFieldIDs.to_csv('NorthFieldIDs.csv', header=True, index=False)
MiddleFielIDs.to_csv('MiddleFieldIDs.csv', header=True, index=False)
SouthFielIDs.to_csv('SouthFieldIDs.csv', header=True, index=False)



In [ ]:

Plots



In [20]:

    
fig, ax = plt.subplots()
m = Basemap(#llcrnrlat=-32., llcrnrlon=48.,
            #urcrnrlat=-22., urcrnrlon=58,
            projection='moll', lon_0=0., lat_0=0., ax=ax)
xN, yN = m(dfNorthPatch.fieldRA.apply(np.degrees).values, 
           dfNorthPatch.fieldDec.apply(np.degrees).values)
xM, yM = m(dfMiddlePatch.fieldRA.apply(np.degrees).values, 
           dfMiddlePatch.fieldDec.apply(np.degrees).values)
xS, yS = m(dfSouthPatch.fieldRA.apply(np.degrees).values, 
           dfSouthPatch.fieldDec.apply(np.degrees).values)



m.scatter(xN, yN, marker='.', color='k')
m.scatter(xM, yM, marker='.', color='r')
m.scatter(xS, yS, marker='.', color='b')









    Out[20]:





<matplotlib.collections.PathCollection at 0x10e454310>



In [21]:

    
fig.savefig('DecBands.png')



In [ ]:



In [ ]:



In [ ]:

Checks

Check that fieldIDs don't overlap in the distinct bands



In [22]:

    
NF = set(NorthFieldIDs.values.tolist())
MF = set(MiddleFielIDs.values.tolist())
SF = set(SouthFielIDs.values.tolist())



In [23]:

    
print(NF.intersection(MF))
print(NF.intersection(SF))
print(SF.intersection(MF))









    



set([])
set([])
set([])

Quick approximation : Is the total number of fields the same in each of the patches?



In [24]:

    
numFieldsNorthWFD/np.float(numFieldsNorth), numFieldsMiddleWFD/np.float(numFieldsMiddle), numFieldsSouthWFD/np.float(numFieldsSouth),









    Out[24]:





(0.9555288461538461, 0.9353658536585366, 0.8957816377171216)



In [25]:

    
numFieldsMiddleWFD / np.float(numFieldsNorthWFD) , numFieldsSouthWFD/np.float(numFieldsNorthWFD)









    Out[25]:





(0.9647798742138365, 0.9081761006289308)



In [26]:

    
numFieldsMiddle / np.float(numFieldsNorth) , numFieldsSouth/np.float(numFieldsNorth)









    Out[26]:





(0.9855769230769231, 0.96875)

The number of total unique fields seems to decrease as we go South in bands despite the choice of wider bands in the South to keep the areas equal by about 3 percent (26 fields) in going from the Northern band to the Souther band
The number of WFD deep fields in the dec bands decrease faster as we go South, falling by about 10 percent (approximately 70 fields) as we go from the Northern most band to the Southern most band, implying that the visits in the South are to different proposals
As we check below, this is due to many of the visits being to other proposals. One that accounts for about half of the extra fields is the Galactic plane (propID 52)



In [27]:

    
x, nn1, nn2 = numFieldsinPatch('@b1dec', '@maxdec', patchName='North', propID=52)









    



fieldDec > @b1dec and fieldDec < @maxdec
Total Number of fields in North patch is 832 and 37 are in WFD



In [28]:

    
x, nm1, nm2 = numFieldsinPatch('@b2dec', '@b1dec', patchName='Middle', propID=52)









    



fieldDec > @b2dec and fieldDec < @b1dec
Total Number of fields in Middle patch is 820 and 53 are in WFD



In [29]:

    
x, ns1, ns2 = numFieldsinPatch('@mindec','@b2dec', patchName='South', propID=52)









    



fieldDec > @mindec and fieldDec < @b2dec
Total Number of fields in South patch is 806 and 84 are in WFD

We will ignore the smaller number of fields in the Southern bands for now

Time for OpSim

We will consider visits in the last nine years of the survey to account for the fact that the first year is likely different.



In [30]:

    
dfall = opsout.summary.query('night > 365')



In [31]:

    
wfd = dfall.query('propID == 54')
numWFD = len(wfd)



In [32]:

    
numOther = len(dfall) - len(wfd)



In [33]:

    
numOther/ np.float(numWFD)









    Out[33]:





0.11913749198460225



In [34]:

    
np.float(numWFD) / np.float(len(dfall))









    Out[34]:





0.8935452588820593

So, I suppose a quick approximaton will be to do only the WFD survey during the last three years, and downsample to 89 percent of the vists. Note that galactic plane, south Celestial cap are done in the first year itself. So these 9 years are shared between the WFD, DDF, and North Ecliptic Spur.

Random Fields



In [39]:

    
np.radians(1.75)









    Out[39]:





0.030543261909900768



In [35]:

    
opsout.summary.columns









    Out[35]:





Index([u'sessionID', u'propID', u'fieldID', u'fieldRA', u'fieldDec', u'filter',
       u'expDate', u'expMJD', u'night', u'visitTime', u'visitExpTime',
       u'finRank', u'FWHMeff', u'FWHMgeom', u'transparency', u'airmass',
       u'vSkyBright', u'filtSkyBrightness', u'rotSkyPos', u'rotTelPos', u'lst',
       u'altitude', u'azimuth', u'dist2Moon', u'solarElong', u'moonRA',
       u'moonDec', u'moonAlt', u'moonAZ', u'moonPhase', u'sunAlt', u'sunAz',
       u'phaseAngle', u'rScatter', u'mieScatter', u'moonIllum', u'moonBright',
       u'darkBright', u'rawSeeing', u'wind', u'humidity', u'slewDist',
       u'slewTime', u'fiveSigmaDepth', u'ditheredRA', u'ditheredDec'],
      dtype='object')



In [42]:

    
FieldIDs = opsout.summary.query('propID == 54').drop_duplicates(subset='fieldID')[['fieldID', 'fieldRA', 'fieldDec']]
FieldIDs[['fieldRA', 'fieldDec']] = FieldIDs[['fieldRA', 'fieldDec']].apply(np.degrees) 
FieldIDs['width'] = 0.03
rng = np.random.RandomState(0)
FieldIDs = FieldIDs.sample(frac=1, replace=False, random_state=rng)
groups = np.array_split(FieldIDs, 3)



In [43]:

    
groups[0].fieldID.size









    Out[43]:





765

Check



In [ ]:

    
- FieldIDs.fieldID.size



In [44]:

    
fig, ax = plt.subplots()
m = Basemap(#llcrnrlat=-32., llcrnrlon=48.,
            #urcrnrlat=-22., urcrnrlon=58,
            projection='moll', lon_0=0., lat_0=0., ax=ax)
xN, yN = m(groups[0].fieldRA.values, 
           groups[0].fieldDec.values)
xM, yM = m(groups[1].fieldRA.values, 
           groups[1].fieldDec.values)
xS, yS = m(groups[2].fieldRA.values, 
           groups[2].fieldDec.values)

m.scatter(xN, yN, marker='.', color='k')
m.scatter(xM, yM, marker='.', color='b')
m.scatter(xS, yS, marker='.', color='r')









    Out[44]:





<matplotlib.collections.PathCollection at 0x10fba51d0>



In [45]:

    
fig.savefig('random_fields.png')



In [46]:

    
groups[0].to_csv('Group_0.csv')
groups[1].to_csv('Group_1.csv')
groups[2].to_csv('Group_2.csv')

Scratch



In [ ]:

    
dfall.propID.unique()



In [ ]:

    
opsout.propIDDict



In [ ]:

    
opsout.summary.query('propID == 53').night.max()



In [ ]:

    
fig, ax = plt.subplots()
opsout.summary.query('propID == 55').night.hist(histtype='step', 
                                                bins=np.arange(0., 3650., 10.),
                                                lw=2., ax=ax)
opsout.summary.query('propID == 52').night.hist(histtype='step', 
                                                bins=np.arange(0., 3650., 10.),
                                                lw=2.,ax =ax)
opsout.summary.query('propID == 53').night.hist(histtype='step', 
                                                bins=np.arange(0., 3650., 10.),
                                                lw=2.,ax=ax)
opsout.summary.query('propID == 54').night.hist(histtype='step', 
                                                bins=np.arange(0., 3650., 10.),
                                                lw=2., ax=ax)
opsout.summary.query('propID == 56').night.hist(histtype='step', 
                                                bins=np.arange(0., 3650., 10.),
                                                lw=2., ax=ax)



In [ ]:

    
fig, ax = plt.subplots()
wfd.slewTime.hist(histtype='step', color='k', lw=2., alpha=1., ax=ax, bins=np.arange(0., 160., 5.), normed=1)
dfall.query('propID == 52').slewTime.hist(histtype='step', color='r', 
                                          ls='dashed', lw=2., ax=ax,
                                          bins=np.arange(0., 160., 5.), normed=1)



In [ ]:

    
np.pi * 8.4 * 1000 /2.



In [ ]:

    
opsout.summary.plot(x='slewDist', y='slewTime', kind='scatter')



In [ ]:

    
import seaborn as sns
sns.set_style('whitegrid')



In [ ]:

    
opsout.summary.query('slewTime == 0.')



In [40]:

    
np.degrees(-1.082474)









    Out[40]:





-62.021191632644268



In [ ]: