In [1]:
import os
import numpy as np
import pandas as pd
Note For this to work, you will need the lsst.sims stack to be installed.
healpy which is installed with the sims stack, but also available from conda
In [2]:
import opsimsummary as oss
from opsimsummary import Tiling, HealpixTiles
# import snsims
import healpy as hp
In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
In [4]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
noTile = snsims.Tiling()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-9-5f6f8a94508e> in <module>()
----> 1 noTile = snsims.Tiling()
TypeError: Can't instantiate abstract class Tiling with abstract methods __init__, area, pointingSequenceForTile, tileIDSequence, tileIDsForSNThe class snsims.Tiling is an abstract Base class. Therefore, this cannot be instantiated. It must be subclassed, and the set of methods outlined have to be implemented for this to work.
In [5]:
class NoTile(Tiling):
pass
In [6]:
noTile = NoTile()
"""noTile = NoTile()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-4-8ddedac7fb97> in <module>()
----> 1 noTile = NoTile()
TypeError: Can't instantiate abstract class NoTile with abstract methods __init__, area, pointingSequenceForTile, positions, tileIDSequence, tileIDsForSN
"""The above fails because the methods are not implemented. Below is a stupid (ie. not useful) but minimalist class that would work:
In [7]:
class MyTile(Tiling):
def __init__(self):
pass
@property
def tileIDSequence(self):
return np.arange(100)
def tileIDsForSN(self, ra, dec):
x = ra + dec
y = np.remainder(x, 100.)
return np.floor(y)
def area(self, tileID):
return 1.
def pointingSequenceForTile(self, tileID, pointings):
return None
def positions(self):
pass
In [8]:
myTile = MyTile()
Currently there is only concrete tiling class that has been implemented. This is the snsims.HealpixTiles class.
This shows how to use the HealpixTiles Class
In [9]:
issubclass(HealpixTiles, Tiling)
Out[9]:
In [10]:
help(HealpixTiles)
In [11]:
datadir = os.path.join(oss.__path__[0], 'example_data')
opsimdb = os.path.join(datadir, 'enigma_1189_micro.db')
healpixelizedDB = os.path.join(datadir, 'healpixels_micro.db')
#opsimdb = '/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db'
In [12]:
NSIDE = 256
In [13]:
hpOpSim = oss.HealPixelizedOpSim.fromOpSimDB(opsimdb, NSIDE=NSIDE)
In [18]:
opsout = oss.OpSimOutput.fromOpSimDB(opsimdb)
In [14]:
#x = np.ones(hp.nside2npix(256)) * hp.UNSEEN
#x[1] = -1
x = np.arange(hp.nside2npix(256))
hp.mollview(x, nest=True)
In [15]:
hptiles = HealpixTiles(healpixelizedOpSim=hpOpSim, nside=NSIDE)
In [16]:
from opsimsummary import HPTileVis
In [20]:
hptvis = HPTileVis(hptiles, opsout)
In [21]:
from opsimsummary import pixelsForAng, plot_south_steradian_view
In [22]:
pixelsForAng(54., -27.5, 4)
Out[22]:
In [23]:
ra, dec = hptvis.pointingCenters(0)
In [24]:
figoverlaps = hptvis.plotTilePointings(135, projection='cyl', paddingFactors=1,
query='night < 365',
**dict(fill=False, color='g', alpha=0.1))
In [25]:
hpchips = HealpixTiles(nside=128, preComputedMap=healpixelizedDB)
In [29]:
hpchips.pointingSequenceForTile(0)
Out[29]:
In [60]:
np.unique(hpchips.pointingSequenceForTile(0)).size
Out[60]:
In [32]:
from opsimsummary import convertToSphericalCoordinates
In [31]:
angs = opsout.summary[['fieldRA', 'fieldDec']]
ra, dec = angs.fieldRA.values, angs.fieldDec.values
In [33]:
convertToSphericalCoordinates??
In [34]:
theta, phi = convertToSphericalCoordinates(ra, dec, unit='radians')
In [35]:
vecs = hp.ang2vec(theta, phi)
In [37]:
vec = hp.pix2vec(128, 0, nest=True)
In [49]:
opsout.summary.iloc[np.abs(np.dot(vecs, vec)).argmin()].fieldID
Out[49]:
In [50]:
obsHistIDsStatic = opsout.summary.query('fieldID==1859').index.values
In [52]:
def staticPoints(tileID):
return obsHistIDsStatic
In [53]:
staticPoints(0)
Out[53]:
In [54]:
hpchipsStatic = HealpixTiles(nside=128, preComputedMap=staticPoints)
In [62]:
np.unique(hpchipsStatic.pointingSequenceForTile(0))
Out[62]:
In [24]:
hpvis = HPTileVis(hpchips, opsout)
In [25]:
fig2, tc, _ = hpvis.plotTilePointings(140785, projection='cea', paddingFactors=4,
drawPointings=True,
**dict(fill=False, color='g', alpha=1.0))
In [47]:
figtile.savefig('HealpixTile_Pointings.pdf')
In [276]:
from opsimsummary import convertToCelestialCoordinates, healpix_boundaries, pixelsForAng
from matplotlib.patches import Polygon
In [371]:
class HPTileVis(HealpixTiles):
def __init__(self, hpTile, opsout):
"""
"""
self.hpTile = hpTile
self.opsout = opsout
self.nside = self.hpTile.nside
def tileIDfromCelestialCoordinates(self, ra, dec, opsout, units='degrees'):
"""
Parameters
-----------
ra :
dec :
units: {'degrees', 'radians'}
"""
return pixelsForAng(lon=ra,lat=dec, unit=units )
def tileCenter(self, tileID):
theta, phi = hp.pix2ang(self.nside, tileID, nest=True)
ra, dec = convertToCelestialCoordinates(theta, phi, input_unit='radians',
output_unit='degrees')
return ra, dec
def pointingSummary(self, tileID, columns=('ditheredRA', 'ditheredDec'),
allPointings=None):
#if allPointings is None:
# allPointings = opsout.summary
# if query is not None:
# allPointings = allPointings.query(query)
# allPointings = allPointings.index.values
obsHistIDs = self.hpTile.pointingSequenceForTile(tileID=tileID,
allPointings=allPointings)
return self.opsout.summary.ix[obsHistIDs]#[list(columns)]
def pointingCenters(self,
tileID,
raCol='ditheredRA',
decCol='ditheredDec',
query=None):
summary = self.pointingSummary(tileID)#, columns=[raCol, decCol])
if query is not None:
summary = summary.query(query)
ra = summary[raCol].apply(np.degrees).values
dec = summary[decCol].apply(np.degrees).values
return ra, dec
def plotTilePointings(self, tileID, raCol='ditheredRA', decCol='ditheredDec', radius=1.75,
paddingFactors=1, query=None, ax=None, projection='cyl',**kwargs):
"""
Parameters
----------
"""
if ax is None:
fig, ax = plt.subplots()
padding = np.degrees(hp.max_pixrad(self.nside)) + radius
ra_tile, dec_tile = self.tileCenter(tileID)
llcrnrlat = dec_tile - padding * paddingFactors
urcrnrlat = dec_tile + padding * paddingFactors
llcrnrlon = ra_tile - padding * paddingFactors
urcrnrlon = ra_tile + padding * paddingFactors
m = Basemap(llcrnrlat=llcrnrlat, llcrnrlon=llcrnrlon,
urcrnrlat=urcrnrlat, urcrnrlon=urcrnrlon,
projection=projection, lon_0=ra_tile, lat_0=dec_tile,
ax=ax)
parallels = np.linspace(llcrnrlat, urcrnrlat, 3)
meridians = np.linspace(llcrnrlon, urcrnrlon, 3)
m.drawparallels(parallels, labels=(1, 0, 0, 0)) #np.ones(len(parallels), dtype=bool))
m.drawmeridians(meridians, labels=(0, 1, 1, 1)) #np.ones(len(meridians), dtype=bool))
ra, dec = self.pointingCenters(tileID, raCol=raCol, decCol=decCol, query=query)
lon, lat = healpix_boundaries(tileID, nside=self.nside, units='degrees',
convention='celestial', step=10,
nest=True)
x, y = m(lon, lat)
xy = zip(x, y)
healpixels = Polygon(xy, facecolor='w',fill=False, alpha=1., edgecolor='k', lw=2)
for ra, dec in zip(ra, dec):
m.tissot(ra, dec, radius, 100, **kwargs)
ax.add_patch(healpixels)
return fig
In [365]:
tileID = pixelsForAng(54., -27.5, 4, unit='degrees')
In [366]:
tileID
Out[366]:
In [367]:
theta, phi = hp.pix2ang(hptiles.nside, 0, nest=True)
ra, dec = oss.convertToCelestialCoordinates(theta, phi)
In [368]:
Out[368]:
In [372]:
hptvis = HPTileVis(hptiles, opsout)
In [388]:
from palettable.colorbrewer import sequential
In [389]:
import palettable
In [ ]:
palettable
In [ ]:
from matplotlib.colorbar import
In [81]:
hptvis.plotTilePointings(, raCol='ditheredRA', decCol='ditheredDec', projection='gnom',
query=None,#'night <3650',
**dict(fill=False, color='g', alpha=1., lw=1., ls='solid'))
Out[81]:
In [376]:
hptvis.centers(0)
In [115]:
opsout.summary.ix[hpTileshpOpSim.pointingSequenceForTile(0, allPointings=None)][['ditheredRA', 'ditheredDec']]
Out[115]:
In [ ]:
In [16]:
phi, theta = hpTileshpOpSim.positions(1, 10000)
In [ ]:
In [17]:
mapvals = np.ones(hp.nside2npix(NSIDE)) * hp.UNSEEN
In [18]:
mapvals[1] = 100
In [ ]:
In [19]:
hp.ang2pix(NSIDE, np.radians(theta), np.radians(phi), nest=True)
Out[19]:
In [20]:
theta_c, phi_c = hp.pix2ang(4, 1, nest=True)
In [21]:
hp.mollview(mapvals, nest=True)
hp.projscatter(np.radians(theta), np.radians(phi), **dict(s=0.0002))
hp.projscatter(theta_c, phi_c, **dict(s=8., c='r'))
Out[21]:
In [22]:
%timeit hpTileshpOpSim.pointingSequenceForTile(33, allPointings=None)
In [23]:
preCompMap = os.path.join(oss.__path__[0], 'example_data', 'healpixels_micro.db')
In [26]:
import os
os.path.exists(preCompMap)
Out[26]:
In [24]:
hpTilesMap = HealpixTiles(nside=1, preComputedMap=preCompMap)
In [25]:
obsHistIDs = hpTilesMap.pointingSequenceForTile(10, allPointings=None)
In [28]:
from sqlalchemy import create_engine
In [36]:
engine = create_engine('sqlite:////' + preCompMap)
In [37]:
import pandas as pd
In [38]:
df = pd.read_sql_query('SELECT * FROM simlib Limit 5', engine)
In [35]:
preCompMap
Out[35]:
In [35]:
%timeit hpOpSim.obsHistIdsForTile(34)
In [36]:
hpTiles = HealpixTiles(healpixelizedOpSim=hpOpSim)
In [37]:
hpTiles.pointingSequenceForTile(34, allPointings=None)
Out[37]:
In [390]:
df = opsout.summary.copy()
In [393]:
df.fieldID.unique().size
Out[393]:
In [398]:
df.query('fieldID == 744')[['fieldRA', 'fieldDec', 'ditheredRA', 'ditheredDec', 'expMJD', 'filter']].head()
Out[398]:
In [399]:
df.query('fieldID == 744')['ditheredDec'] = df.query('fieldID == 744')['fieldDec']
In [405]:
def fixdec(mydf, fieldIDs):
return mydf.fieldID.apply(lambda x: x in fieldIDs)
In [407]:
fixdec(df, (744, 1427, 290)).astype(np.int) * df.fieldRA + (1.0 - fixdec(df, (744, 1427, 290)).astype(np.int))
Out[407]:
In [400]:
df.query('fieldID == 744')[['fieldRA', 'fieldDec', 'ditheredRA', 'ditheredDec', 'expMJD', 'filter']].head()
Out[400]:
In [54]:
df = opsout.summary.copy()
In [55]:
grouped = df.groupby('propID')
In [59]:
grouped.groups.keys()
Out[59]:
In [62]:
grouped.get_group(366).fieldRA
Out[62]:
In [ ]: