<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  Declination Boundaries for the WFD proposal
Step2:  Total area covered, should be ~ 0.5 sky
Step3:  Equations for the boundary values
Step4:  Plots
Step5:  Checks
Check that fieldIDs don't overlap in the distinct bands
Step6:  Quick approximation 
Step7:  So
- 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)
Step8:  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.
Step9:  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
Step10:  Check
Step11:  Scratch
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from opsimsummary import OpSimOutput
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
opsout = OpSimOutput.fromOpSimDB('/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db',
                                subset='_all')
opsout.summary.expMJD.size
opsout.summary.reset_index().drop_duplicates(subset='obsHistID').expMJD.size
Explanation: 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
End of explanation
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
Explanation: Declination Boundaries for the WFD proposal
End of explanation
area = 2 * np.pi * (np.cos(mintheta) - np.cos(maxtheta))
# Check that the sky fraction is ~ 0.5
area / 4 / np.pi
Explanation: Total area covered, should be ~ 0.5 sky
End of explanation
b1 = np.arccos((2 * np.pi * np.cos(mintheta) - area /3)/ 2 /np.pi)
b1dec = - b1 + np.pi/2
np.degrees(b1dec)
b2 = np.arccos(( area /3 + 2 * np.pi * np.cos(maxtheta) )/ 2 /np.pi)
b2dec = -b2 + np.pi/2.
np.degrees(b2dec)
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
dfNorthPatch, numFieldsNorth, numFieldsNorthWFD = numFieldsinPatch('@b1dec', '@maxdec', patchName='North')
dfMiddlePatch, numFieldsMiddle, numFieldsMiddleWFD = numFieldsinPatch('@b2dec', 
                                                                   '@b1dec', patchName='Middle')
dfSouthPatch, numFieldsSouth, numFieldsSouthWFD = numFieldsinPatch('@mindec', 
                                                                   '@b2dec', patchName='South')
NorthFieldIDs = dfNorthPatch.query('propID == 54').fieldID
MiddleFielIDs = dfMiddlePatch.query('propID == 54').fieldID
SouthFielIDs = dfSouthPatch.query('propID == 54').fieldID
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)
Explanation: 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.
End of explanation
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')
fig.savefig('DecBands.png')
Explanation: Plots
End of explanation
NF = set(NorthFieldIDs.values.tolist())
MF = set(MiddleFielIDs.values.tolist())
SF = set(SouthFielIDs.values.tolist())
print(NF.intersection(MF))
print(NF.intersection(SF))
print(SF.intersection(MF))
Explanation: Checks
Check that fieldIDs don't overlap in the distinct bands
End of explanation
numFieldsNorthWFD/np.float(numFieldsNorth), numFieldsMiddleWFD/np.float(numFieldsMiddle), numFieldsSouthWFD/np.float(numFieldsSouth), 
numFieldsMiddleWFD / np.float(numFieldsNorthWFD) , numFieldsSouthWFD/np.float(numFieldsNorthWFD)
numFieldsMiddle / np.float(numFieldsNorth) , numFieldsSouth/np.float(numFieldsNorth)
Explanation: Quick approximation : Is the total number of fields the same in each of the patches?
End of explanation
x, nn1, nn2 = numFieldsinPatch('@b1dec', '@maxdec', patchName='North', propID=52)
x, nm1, nm2 = numFieldsinPatch('@b2dec', '@b1dec', patchName='Middle', propID=52)
x, ns1, ns2 = numFieldsinPatch('@mindec','@b2dec', patchName='South', propID=52)
Explanation: So
- 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)
End of explanation
dfall = opsout.summary.query('night > 365')
wfd = dfall.query('propID == 54')
numWFD = len(wfd)
numOther = len(dfall) - len(wfd)
numOther/ np.float(numWFD) 
np.float(numWFD) / np.float(len(dfall))
Explanation: 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.
End of explanation
np.radians(1.75)
opsout.summary.columns
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)
groups[0].fieldID.size
Explanation: 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
End of explanation
  - FieldIDs.fieldID.size
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')
fig.savefig('random_fields.png')
groups[0].to_csv('Group_0.csv')
groups[1].to_csv('Group_1.csv')
groups[2].to_csv('Group_2.csv')
Explanation: Check
End of explanation
dfall.propID.unique()
opsout.propIDDict
opsout.summary.query('propID == 53').night.max()
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)
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)
np.pi * 8.4 * 1000 /2.
opsout.summary.plot(x='slewDist', y='slewTime', kind='scatter')
import seaborn as sns
sns.set_style('whitegrid')
opsout.summary.query('slewTime == 0.')
np.degrees(-1.082474)
Explanation: Scratch
End of explanation
<END_TASK>