<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 is a notebook to explore opSim outputs in different ways, mostly useful to supernova analysis. We will look at the opsim output called Enigma_1189.
Using the notebook requires installing
- lsst_sims
-
Step1:  Read in OpSim output for modern versions
Step2:  Read in the OpSim DataBase into a pandas dataFrame
Step3:  The opsim database is a large file (approx 4.0 GB), but still possible to read into memory on new computers. You usually only need the Summary Table, which is about 900 MB. If you are only interested in the Deep Drilling Fields, you can use the read_sql_query to only select information pertaining to Deep Drilling Observations. This has a memory footprint of about 40 MB.
Obviously, you can reduce this further by narrowing down the columns to those of interest only. For the entire Summary Table, this step takes a few minutes on my computer. 
If you are going to do the read from disk step very often, you can further reduce the time used by storing the output on disk as a hdf5 file and reading that into memory
We will look at three different Summaries of OpSim Runs. A summary of the 
1. Deep Drilling fields
Step4:  If we knew ahead of time the proposal ID, then we could have done this quicker using
Step5:  We can also sub-select this from the all-encompassing Summay Table. This can be done in two way
Step6:  First Season
We can visualize the cadence during the first season using the cadence plot for a particular field
Step7:  Suppose we have a supernova with a peak around a particular MJD of 49540, and we want to see what the observations happened around it
Step8:  Scratch
<ASSISTANT_TASK:>
Python Code:
import numpy as np
%matplotlib inline 
import matplotlib.pyplot as plt
import sncosmo
import os
import OpSimSummary
print OpSimSummary.__file__
# Required packages sqlachemy, pandas (both are part of anaconda distribution, or can be installed with a python installer)
# One step requires the LSST stack, can be skipped for a particular OPSIM database in question
import OpSimSummary.summarize_opsim as so
from sqlalchemy import create_engine
import pandas as pd
print so.__file__
# This step requires LSST SIMS package MAF. The main goal of this step is to set DD and WFD to integer keys that 
# label an observation as Deep Drilling or for Wide Fast Deep.
# If you want to skip this step, you can use the next cell by uncommenting it, and commenting out this cell, if all you
# care about is the database used in this example. But there is no guarantee that the numbers in the cell below will work
# on other versions of opsim database outputs
from lsst.sims.maf import db
from lsst.sims.maf.utils import opsimUtils
from LSSTmetrics import PerSNMetric
import LSSTmetrics
print LSSTmetrics.__file__
from lsst.sims.photUtils import BandpassDict
# DD = 366
# WFD = 364
Explanation: This is a notebook to explore opSim outputs in different ways, mostly useful to supernova analysis. We will look at the opsim output called Enigma_1189.
Using the notebook requires installing
- lsst_sims
-
End of explanation
# Change dbname to point at your own location of the opsim output
dbname = '/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db'
opsdb = db.OpsimDatabase(dbname)
propID, propTags = opsdb.fetchPropInfo()
DD = propTags['DD'][0]
WFD = propTags['WFD'][0]
print("The propID for the Deep Drilling Field {0:2d}".format(DD))
print("The propID for the Wide Fast Deep Field {0:2d}".format(WFD))
Explanation: Read in OpSim output for modern versions: (sqlite formats)
Description of OpSim outputs are available on the page https://confluence.lsstcorp.org/display/SIM/OpSim+Datasets+for+Cadence+Workshop+LSST2015http://tusken.astro.washington.edu:8080
Here we will use the opsim output http://ops2.tuc.noao.edu/runs/enigma_1189/data/enigma_1189_sqlite.db.gz
I have downloaded this database, unzipped and use the variable dbname to point to its location
End of explanation
engine = create_engine('sqlite:///' + dbname)
Explanation: Read in the OpSim DataBase into a pandas dataFrame
End of explanation
# Load to a dataframe
# Summary = pd.read_hdf('storage.h5', 'table')
# Summary = pd.read_sql_table('Summary', engine, index_col='obsHistID')
# EnigmaDeep  = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is 366', engine)
# EnigmaD  = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is 366', engine)
Explanation: The opsim database is a large file (approx 4.0 GB), but still possible to read into memory on new computers. You usually only need the Summary Table, which is about 900 MB. If you are only interested in the Deep Drilling Fields, you can use the read_sql_query to only select information pertaining to Deep Drilling Observations. This has a memory footprint of about 40 MB.
Obviously, you can reduce this further by narrowing down the columns to those of interest only. For the entire Summary Table, this step takes a few minutes on my computer. 
If you are going to do the read from disk step very often, you can further reduce the time used by storing the output on disk as a hdf5 file and reading that into memory
We will look at three different Summaries of OpSim Runs. A summary of the 
1. Deep Drilling fields: These are the observations corresponding to propID of the variable DD above, and are restricted to a handful of fields
2. WFD (Main) Survey: These are the observations corresponding to the propID of the variables WFD
3. Combined Survey: These are observations combining DEEP and WFD in the DDF. Note that this leads to duplicate observations which must be subsequently dropped.
End of explanation
OpSim_combined = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is ' + str(DD) + ' or ' + str(WFD), engine, index_col='obsHistID')
OpSim_Deep = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is ' + str(DD), engine, index_col='obsHistID')
Explanation: If we knew ahead of time the proposal ID, then we could have done this quicker using
End of explanation
OpSimDeepSummary = so.SummaryOpsim(OpSim_Deep)
OpSimCombinedSummary = so.SummaryOpsim(OpSim_combined)
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection='mollweide');
fig = OpSimDeepSummary.showFields(ax=fig.axes[0], marker='o', s=40)
OpSimCombinedSummary.showFields(ax=ax, marker='o', color='r', s=8)
#fieldList = EnigmaDeepSummary.fieldIds
Explanation: We can also sub-select this from the all-encompassing Summay Table. This can be done in two way:
Some properties of the OpSim Outputs
Construct our Summary
End of explanation
firstSeasonDeep = OpSimDeepSummary.cadence_plot(fieldID=1427, observedOnly=False, sql_query='night < 366')
firstSeasonCombined = OpSimCombinedSummary.cadence_plot(fieldID=1427, observedOnly=False, sql_query='night < 366')
OpSimCombinedSummary.mjdvalfornight(300)
firstSeasonCombined[0].savefig('minion_1427.pdf')
firstSeason = OpSimCombinedSummary.cadence_plot(fieldID=1430, observedOnly=False, sql_query='night <366', 
                                             nightMin=0, nightMax=366)
Explanation: First Season
We can visualize the cadence during the first season using the cadence plot for a particular field: The following plot shows how many visits we have in different filters on a particular night:
End of explanation
SN = OpSimCombinedSummary.cadence_plot(fieldID=1427, #racol='fieldRA', deccol='fieldDec',
                                                  observedOnly=False, mjd_center=59580 + 270., mjd_range=[-30., 50.])
# ax = plt.gca()
# ax.axvline(49540, color='r', lw=2.)
# ax.xaxis.get_major_formatter().set_useOffset(False)
SN[0].savefig("Minion_1427_59850.pdf")
lsst_bp = BandpassDict.loadTotalBandpassesFromFiles()
# sncosmo Bandpasses required for fitting
throughputsdir = os.getenv('THROUGHPUTS_DIR')
from astropy.units import Unit
bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
for band in bandPassList:
    # setup sncosmo bandpasses
    bandfname = banddir + "/total_" + band + '.dat'
    # register the LSST bands to the SNCosmo registry
    # Not needed for LSST, but useful to compare independent codes
    # Usually the next two lines can be merged,
    # but there is an astropy bug currently which affects only OSX.
    numpyband = np.loadtxt(bandfname)
    print(band)
    sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
                                   trans=numpyband[:, 1],
                                   wav\e_unit=Unit('nm'),
                                   name=band)
    sncosmo.registry.register(sncosmoband, force=True)
snDaily = PerSNMetric(summarydf=OpSimCombinedSummary.simlib(1427),t0=59580 + 270., 
                      raCol='fieldRA', decCol='fieldDec', lsst_bp=lsst_bp)
snDaily.lightcurve
notCoadded = snDaily.lcplot(scattered=True)
CoaddedNights = snDaily.lcplot(nightlyCoadd=True, scattered=True)
notCoadded.savefig('SingleVisitNights.pdf')
CoaddedNights.savefig('CoaddedNights.pdf')
!pwd
OpSimCombinedSummary.df.night.min()
SN[0].savefig('SN_observaton.pdf')
Explanation: Suppose we have a supernova with a peak around a particular MJD of 49540, and we want to see what the observations happened around it:
End of explanation
SN_matrix.sum(axis=1).sum()
EnigmaDeep.query('fieldID == 744 and expMJD < 49590 and expMJD > 49510').expMJD.size
SN_matrix[SN_matrix > 0.5] = 1
SN_matrix.sum().sum()
len(SN_matrix.sum(axis=1).dropna())
nightlySN_matrix = SN_matrix.copy(deep=True)
nightlySN_matrix[SN_matrix > 0.5] =1
nightlySN_matrix.sum(axis=1).dropna().sum()
nightlySN_matrix.sum(axis=1).dropna().size
nightlySN_matrix.sum(ax)
Explanation: Scratch
End of explanation
<END_TASK>