In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
# 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 simlibs.summarize_opsim as so
from sqlalchemy import create_engine
import pandas as pd
print so.__file__
In [3]:
# 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
In [4]:
# DD = 366
# WFD = 364
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
In [4]:
dbname = 'sqlite:////Users/rbiswas/data/LSST/OpSimData/enigma_1189_sqlite.db'
opsdb = db.OpsimDatabase(dbname)
propID, propTags = opsdb.fetchPropInfo()
DD = propTags['DD'][0]
WFD = propTags['DD'][0]
In [5]:
engine = create_engine(dbname)
In [6]:
dbname_2168 = 'sqlite:////Users/rbiswas/data/LSST/OpSimData/opsim2_168_sqlite.db'
engine_2168 = create_engine(dbname_2168)
dbAddressDict = {'dbAddress': dbname_2168, 'Summary':'summary'}
print dbAddressDict
opsimDb2168 = opsimUtils.connectOpsimDb(dbAddressDict)
In [7]:
propID, propTags = opsimDb2168.fetchPropInfo()
In [9]:
propTags
Out[9]:
This does not work for the old format. Read the SSTar document: http://opsimcvs.tuc.noao.edu/runs/opsim2.168/SSTAR-opsim2-168.pdf and find the deep drilling propIDs from Table 1.
In [10]:
DDold = 359, 360
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.
In [6]:
# Load to a dataframe
# Summary = pd.read_hdf('storage.h5', 'table')
# Summary = pd.read_sql_table('Summary', engine)
EnigmaDeep = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is 366', engine)
In [7]:
EnigmaDeep.fieldID.unique()
Out[7]:
In [8]:
# EnigmaDeepFields = pd.read_sql_table('Summary', engine)
# EnigmaDeepFields = EnigmaDeepFields.drop_duplicates()
In [8]:
# List columns
EnigmaDeep.columns
Out[8]:
In [16]:
# Alternatively, you can use Summary to read the entire table, and then slice using pandas
# DeepFields = Summary[Summary['propID'] == 366]
# WFD = Summary[Summary['propID'] == 364]
# season1Deep = DeepFields[DeepFields['night'] < 366]
We can read information about the Enigma Deep field summary table downloaded from Opsim with the command below to summarize the table. There are 110114 rows, and a total of 45 columns. The total memory on disk used is 38.6 MB
In [17]:
EnigmaDeep.fieldID.unique()
Out[17]:
In [18]:
OpSim_2168_Deep = pd.read_sql_query('SELECT * FROM SUMMARY WHERE PROPID is 359 or PROPID is 360', engine_2168)
As shown from a similar command as above this OpSim database has 197332 rows (close to about twice as the enigma case) and 46 columns. The memory on disc is also larger by a similar factor and is 70 MB
In [19]:
OpSim_2168_Deep.info()
In [20]:
OpSim_2168_Deep['propID'].unique()
Out[20]:
We can select a season, for example the first season by the following method. It would be trivial to add a column giving a season by dividing the night variable by 365.
In [21]:
Season0_2_168 = OpSim_2168_Deep[OpSim_2168_Deep['night'] < 366]
In [22]:
OpSim_2168_Deep.fieldID.unique()
Out[22]:
This is the point, where we start using this package. The stuff before was simply interacting with OpSim
In [23]:
OpSim_2_168_DeepFull = so.SummaryOpsim(OpSim_2168_Deep)
We can view the summary for a particular field
In [24]:
OpSim_2_168_DeepFull.simlib(2082)[['expMJD','obsHistID','filter', 'night']].head()
Out[24]:
In [25]:
OpSim_2_168_DeepFull.writeSimlib('OpSim_2.168.DEEP.Full')
In [26]:
OpSim_2_168_DeepFull.fieldIds
Out[26]:
In [27]:
fig = OpSim_2_168_DeepFull.cadence_plot(519, sql_query='night <365', nightMin=0, nightMax=365)
In [28]:
fig = map(lambda x: OpSim_2_168_DeepFull.cadence_plot(x, sql_query='night <365', nightMin=0, nightMax=365), OpSim_2_168_DeepFull.fieldIds)
In [29]:
fig = map(lambda x: OpSim_2_168_DeepFull.cadence_plot(x, sql_query='night <3650', nightMin=0, nightMax=3650), OpSim_2_168_DeepFull.fieldIds)
In [7]:
EnigmaDeepSummary = so.SummaryOpsim(EnigmaDeep)
In [12]:
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection='mollweide');
fig = OpSim_2_168_DeepFull.showFields(ax=fig.axes[0], marker='o', s=80)
In [11]:
fig = EnigmaDeepSummary.showFields(ax=fig.axes[0], marker='s', color='r')
In [33]:
ax = fig.gca()
ax.set_title('Deep Drilling Fields for 2.168, Enigma_1189')
Out[33]:
In [34]:
fig
Out[34]:
In [9]:
EnigmaDeepSummary.simlib(290).head()
Out[9]:
In [31]:
mat = EnigmaDeepSummary.cadence_Matrix(290, observedOnly=True)
In [33]:
mat[180:]
Out[33]:
In [25]:
gg = pd.DataFrame({'gg':mat.g.diff().values, 'rr':mat.r.diff().values})
In [30]:
gg['gg'].values[170:190]
Out[30]:
In [9]:
fig = EnigmaDeepSummary.cadence_plot(290, observedOnly=True)
plt.savefig('cadence.pdf')
In [20]:
matrix = EnigmaDeepSummary.cadence_plot(290, observedOnly=True)
In [106]:
axesImage = plt.matshow(matrix.transpose(), aspect='auto', cmap=plt.cm.gray_r, vmin=0., vmax=1.)
# plt.colorbar(orientation='horizontal')
In [107]:
ax = axesImage.axes
In [56]:
fig = ax.figure
In [103]:
EnigmaDeepSummary.dec(290)
Out[103]:
In [120]:
title = 'fieldID: {:0>2d} (ra: {:+3f} dec: {:+3f})'
title = title.format(290, EnigmaDeepSummary.ra(290), EnigmaDeepSummary.dec(290))
print title
In [67]:
ax.set_title(title)
Out[67]:
In [113]:
for i, item in enumerate(ax.get_yticklabels()[1:-1]):
print i, item.get_text(), ax.set_text(Filters[i]), item.get_text()
In [68]:
fig.colorbar(cax=axesImage,)
In [87]:
Filters = [u'u', u'g', u'r', u'i', u'z', u'Y']
In [92]:
len(yticklabels)
Out[92]:
In [91]:
for i, item in enumerate(yticklabels):
item.set_text(Filters[i])
In [89]:
yticklabels = ax.get_yticklabels()
In [93]:
[item.get_text() for item in yticklabels]
Out[93]:
In [81]:
ax.xaxis.tick_bottom()
In [114]:
fig = ax.figure
fig
Out[114]:
In [78]:
xticks = ax.get_xticks()
In [80]:
ax.set_xticks(xticks)
In [118]:
ax.set_yticklabels(['0'] +Filters, minor=False)
Out[118]:
In [119]:
fig
Out[119]:
In [36]:
fig = map(lambda x: EnigmaDeepSummary.cadence_plot(x, sql_query='night < 3650', nightMin=0, nightMax=3650), EnigmaDeepSummary.fieldIds)
In [37]:
EnigmaDeepFieldSummary = so.SummaryOpsim(EnigmaDeepFields)
In [38]:
EnigmaDeepFieldSummary.fieldIds[10]
Out[38]:
In [39]:
fig = EnigmaDeepFieldSummary.cadence_plot(, sql_query='night < 365', nightMax=365)
In [41]:
fig = EnigmaDeepSummary.cadence_plot(744, sql_query='night < 730', nightMin=0, nightMax=730)
In [40]:
fig = map(lambda x: EnigmaDeepFieldSummary.cadence_plot(x, sql_query='night < 730', nightMin=0, nightMax=730), EnigmaDeepFieldSummary.fieldIds)
In [ ]:
season = season1Deep.copy(deep=True)
In [ ]:
oldseason = Season0_2_168.copy(deep=True)
In [ ]:
simlib_opsim_2_168 = so.SummaryOpsim(oldseason)
In [ ]:
simlib_opsim_2_168.fieldIds
In [ ]:
simlib_opsim_2_168.writeSimlib('Opsim.2.168.Simlib')
In [ ]:
simlib = so.SummaryOpsim(season)
In [ ]:
simlib.df['filter'].unique()
In [ ]:
simlib.fieldIds
In [ ]:
simlib.writeSimlib('test_out.simlib')
In [ ]:
filtergroups = EnigmaDeepSummary.simlib(290).query('night < 366').groupby('filter')
In [ ]:
filters = filtergroups.groups.keys()
Filters = [u'u', u'g', u'r', u'i', u'z', u'Y']
In [ ]:
filtered = map(lambda x: filtergroups.get_group(x).copy(deep=True), filters)
In [ ]:
def cadence_plot(dd, fieldID, sql_query='night < 366', Filters=[u'u', u'g', u'r', u'i', u'z', u'Y'], nightMin=0, nightMax=365):
filtergroups = dd.simlib(fieldID).query(sql_query).groupby('filter')
times = dict()
numExps = dict()
numDays = nightMax - nightMin
H = np.zeros(shape=(len(Filters), numDays))
for i, filt in enumerate(Filters):
expVals = np.zeros(numDays, dtype=int)
filtered = map(lambda x: filtergroups.get_group(x).copy(deep=True), Filters)
# filtered[i]['timeIndex'] = filtered[i].expMJD // deltaT
timeBinned = filtered[i].groupby('night')
timeKeys = timeBinned.groups.keys()
times = map(int, timeBinned.night.apply(np.mean))
times = np.array(times) - nightMin
# print times
# times = map(lambda x: timeBinned.get_group(x).expMJD.mean(), timeKeys)
numExps = timeBinned.apply(len)
expVals[times] = numExps
H[i, :] = expVals
#fig = plt.figure()
# ax = fig.add_subplot(111)
ax = plt.matshow(H, aspect='auto')
plt.colorbar(orientation='horizontal', cmap=cm.gray_r)
# ax.axhline(1,0, 365)
fig = ax.figure
# ax.grid(True)
# fig = ax.figure
return fig
In [ ]:
OpSim_2_168_DeepFull.fieldIds
In [ ]:
EnigmaDeepSummary.fieldIds
In [ ]:
In [ ]:
fig = EnigmaDeepSummary.cadence_plot(290, sql_query='night <366', nightMin=0, nightMax=365)
In [ ]:
fig = EnigmaDeepSummary.cadence_plot(1427, sql_query='night <365', nightMin=0, nightMax=365)
In [ ]:
fig = EnigmaDeepSummary.cadence_plot(2412, sql_query='night <365', nightMin=0, nightMax=365)
In [ ]:
fig = EnigmaDeepSummary.cadence_plot(2786, sql_query='night <365', nightMin=0, nightMax=365)
In [ ]:
fig = plt.figure()#figsize=(10*10, 4))
ax = fig.add_subplot(1,3,9)
In [ ]:
help(fig.add_subplot(10,1,1))
In [ ]:
help(plt.matshow)
In [ ]:
np.shape(H)
In [ ]:
plt.matshow(H, aspect='auto')
plt.colorbar(orientation='horizontal', cmap = cm.gray_r)
In [ ]:
help(plt.matshow)
In [ ]:
import matplotlib.cm as cm
In [ ]:
plt.colormaps()
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(times[u'g'], numExps[u'g'])
In [ ]:
expVals = np.zeros(365)
In [ ]:
expVals
In [ ]:
times[u'g']
In [ ]:
t = np.zeros(len(times[u'g'].values), dtype=int)
In [ ]:
t = times[u'g'].values
In [ ]:
t = map(int, times[u'g'].values)
In [ ]:
expVals[t] = numExps[u'g']
In [ ]:
expVals
In [ ]:
numExps[u'g'].dtype
In [ ]:
cadence = filtergroups.get_group(u'g')
In [ ]:
cad = cadence.copy(deep=True)
In [ ]:
deltaT = 1
cad['timeIndex'] = cad.expMJD // deltaT
In [ ]:
gg = cad.groupby('timeIndex').apply(len)
In [ ]:
filters
In [ ]:
filtered[0]
In [ ]:
map(len,gg.values())
In [ ]:
obs = cadence.groupby('night')
In [ ]:
In [ ]:
map(lambda x: EnigmaDeepSummary.ra(x), EnigmaDeepSummary.fieldIds)
In [ ]:
EnigmaDeepSummary.ra(290)
In [ ]:
fieldSummary = EnigmaDeepSummary.simlib(290)
In [ ]:
dd = fieldSummary[fieldSummary['night'] < 366 ]
In [ ]:
xx = fieldSummary.query('night < 366')
In [ ]:
all(xx == dd)
In [ ]:
def nightsObserved(dd, band, color='k', ax=None):
xx = dd.groupby(dd['filter']).get_group(band).groupby('night')
yy = np.sort(xx.groups.keys())
if ax is None:
fig, ax = plt.subplots()
ax.bar(yy, np.ones(len(yy)), color=color)
In [ ]:
def plotnightsObserved(dd):
nightsObserved(dd, u'u', color='y')
nightsObserved(dd, u'g', color='g')
nightsObserved(dd, u'r', color='b')
nightsObserved(dd, u'i', color ='k')
nightsObserved(dd, u'z', color='red')
nightsObserved(dd, u'Y', color='magenta')
In [ ]:
def plotnightsDeep()
In [ ]:
nightsObserved(dd, u'g')
In [ ]:
plotnightsObserved(EnigmaDeepSummary.simlib(744)[EnigmaDeepSummary.simlib(744)['night']<730])
In [ ]:
dd.night.hist(by=dd['filter'], bins=365, histtype='step', sharex=True)
In [ ]:
help(dd.night.hist)
In [ ]:
fd = dd.groupby('filter')
In [ ]:
dp=fd.get_group(u'g')
In [ ]:
dd.groups.keys()
In [ ]:
dp.plot(dp['night'], bins=365, histtype='step')
In [ ]:
dd.get_group(u'g').night.plot(kind='bar')
In [ ]:
xd = dd.get_group(u'g').('expMJD', kind='bar')
In [ ]:
np.sort(xd.groups.keys())
In [ ]:
(map(lambda x: xd.get_group(x).obsHistID.size, xd.groups.keys()))
In [ ]:
def filterDiff(summary, fieldID, filter):
x = summary.simlib(fieldID).groupby('filter').get_group(filter).night.unique()
return np.diff(x)
In [ ]:
plt.hist(filterDiff(OpSim_2_168_DeepFull, 2082, 'g'), histtype='step', bins=20, color='g')
plt.hist(filterDiff(OpSim_2_168_DeepFull, 2082, 'u'), histtype='step', bins=20, color='r')
plt.hist(filterDiff(OpSim_2_168_DeepFull, 2082, 'Y'), histtype='step', bins=20, color='k')
In [ ]:
def drawfields(summ, ax=None, marker=None):
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111, projection='mollweide')
if marker is None:
marker = 'o'
ax.plot(summ.coords()[0], summ.coords()[1], marker);
ax.grid(True)
fig = ax.figure;
return fig
#ax.plot(np.radians([30., 45.]), np.radians([75, 45.]),'o')
In [ ]:
import seaborn as sns
sns.set()
# Load the example flights dataset and conver to long-form
flights_long = sns.load_dataset("flights")
flights = flights_long.pivot("month", "year", "passengers")
# Draw a heatmap with the numeric values in each cell
sns.heatmap(flights, annot=True, fmt="d", linewidths=.5)
In [ ]:
type(flights_long)
In [ ]:
flights_long
In [ ]:
type(flights)
In [ ]:
flights
In [ ]:
gg = EnigmaDeepSummary.simlib(744).groupby(['filter', 'night'])
In [ ]:
f1, nt = zip(*gg.groups.keys())
In [ ]:
gg.apply(len).values
In [ ]:
filts, nights = zip( *gg.groups.keys())
nums = gg.apply(len).values
df = pd.DataFrame({'filts': list(filts), 'nights': list(nights), 'nums': nums})
In [ ]:
df.columns
In [ ]:
dd = df.pivot('nights','filts','nums' )
In [ ]:
dd.columns
In [ ]:
ddd = dd.fillna(0).transpose()
In [ ]:
ss = pd.Series(np.arange(0,365))
In [ ]:
ddd.index
In [ ]:
ddd = ddd.reindex(ss, fill_value=0)
In [ ]:
xx = ddd[[u'u', u'g', u'r', u'i', u'z', u'Y']]
In [ ]:
xx.iloc[180:185]
In [ ]:
ddd.iloc[180:185]
In [ ]:
gg.apply(len).values
In [ ]:
fig = sns.heatmap(ddd.transpose(), annot=True, fmt='.2g', linewidth=0.5)
In [ ]:
def grayify_cmap(cmap):
"""Return a grayscale version of the colormap"""
cmap = plt.cm.get_cmap(cmap)
colors = cmap(np.arange(cmap.N))
# convert RGBA to perceived greyscale luminance
# cf. http://alienryderflex.com/hsp.html
RGB_weight = [0.299, 0.587, 0.114]
luminance = np.sqrt(np.dot(colors[:, :3] ** 2, RGB_weight))
colors[:, :3] = luminance[:, np.newaxis]
return cmap.from_list(cmap.name + "_grayscale", colors, cmap.N)
def show_colormap(cmap):
im = np.outer(np.ones(10), np.arange(100))
fig, ax = plt.subplots(2, figsize=(6, 1.5),
subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(hspace=0.1)
ax[0].imshow(im, cmap=cmap)
ax[1].imshow(im, cmap=grayify_cmap(cmap))
show_colormap('gray')
In [ ]:
tips = sns.load_dataset("tips")
In [ ]:
tips
In [ ]:
sns.FacetGrid(tips, col='time', col_wrap=1)
In [ ]:
fig = plt.figure()
In [ ]:
import matplotlib.gridspec as gridspec
In [ ]:
gs = gridspec.GridSpec(10,1)
In [ ]:
ax0 = fig.add_subplot(gs[0])
In [ ]:
fig = EnigmaDeepSummary.cadence_plot(744, sql_query='night <366', nightMin=0, nightMax=365)
In [ ]:
plt.matshow(xx.transpose(), aspect='auto', cmap=plt.cm.gray_r)
plt.yticks(np.arange(6),('u', 'g', 'r', 'i', 'z', 'Y'))
plt.grid(True)
In [ ]:
np.arange(5, 12)
In [ ]:
(180 / np.pi)**2
In [17]:
l = np.array([1,0, 2, 3, 0, 0, 3])
In [19]:
np.where(l, 1, 0)
Out[19]:
In [25]:
fig, ax = plt.subplots()
#ax = fig.axes[0]
ax.plot(l,l)
ax.set_title('tit;e')
Out[25]:
In [ ]: