LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
multiindex.py
Go to the documentation of this file.
1 import os
2 import numpy
3 import pyfits
4 import lsst.utils
5 from lsst.pex.logging import getDefaultLog
6 from .astrometry_net import multiindex_new, multiindex_add_index, INDEX_ONLY_LOAD_METADATA, healpixDistance
7 from .astrometryNetDataConfig import AstrometryNetDataConfig
8 
9 __all__ = ["getIndexPath", "getConfigFromEnvironment", "AstrometryNetCatalog", "generateCache"]
10 
11 def getIndexPath(fn):
12  """!Get the path to the specified astrometry.net index file
13 
14  No effort is made to confirm that the file exists, so it may be used to locate the
15  path to a non-existent file (e.g., to write).
16 
17  @param[in] fn path to index file; if relative, then relative to astrometry_net_data
18  if that product is setup, else relative to the current working directory
19  @return the absolute path to the index file
20  """
21  if os.path.isabs(fn):
22  return fn
23  try:
24  andir = lsst.utils.getPackageDir('astrometry_net_data')
25  except:
26  # Relative to cwd
27  return os.path.abspath(fn)
28  return os.path.join(andir, fn)
29 
30 
32  """Find the config file from the environment
33 
34  The andConfig.py file is in the astrometry_net_data directory.
35  """
36  try:
37  anDir = lsst.utils.getPackageDir('astrometry_net_data')
38  except:
39  anDir = os.getcwd()
40  andConfigPath = "andConfig.py"
41  if not os.path.exists(andConfigPath):
42  raise RuntimeError("Unable to find andConfig.py in the current directory. "
43  "Did you forget to setup astrometry_net_data?")
44  else:
45  andConfigPath = os.path.join(anDir, "andConfig.py")
46  if not os.path.exists(andConfigPath):
47  raise RuntimeError("Unable to find andConfig.py in astrometry_net_data directory %s" % (anDir,))
48 
49  andConfig = AstrometryNetDataConfig()
50  andConfig.load(andConfigPath)
51  return andConfig
52 
53 class MultiIndexCache(object):
54  """A wrapper for the multiindex_t, which only reads the data when it needs to
55 
56  The MultiIndexCache may be instantiated directly, or via the 'fromFilenameList'
57  class method, which loads it from a list of filenames.
58  """
59  def __init__(self, filenameList, healpix, nside):
60  """!Constructor
61 
62  @param filenameList List of filenames; first is the multiindex, then
63  follows the individual index files
64  @param healpix Healpix number
65  @param nside Healpix nside
66  """
67  if len(filenameList) < 2:
68  raise RuntimeError("Insufficient filenames provided for multiindex (%s): expected >= 2" %
69  (filenameList,))
70  self._filenameList = filenameList
71  self._healpix = int(healpix)
72  self._nside = int(nside)
73  self._mi = None
74  self._loaded = False
75  self.log = getDefaultLog()
76 
77  @classmethod
78  def fromFilenameList(cls, filenameList):
79  """Construct from a list of filenames
80 
81  The list of filenames should contain the multiindex filename first,
82  then the individual index filenames. The healpix and nside are
83  determined by reading the indices, so this is not very efficient.
84  """
85  self = cls(filenameList, 0, 0)
86  self.reload()
87  healpixes = set(self[i].healpix for i in range(len(self)))
88  nsides = set(self[i].hpnside for i in range(len(self)))
89  assert len(healpixes) == 1
90  assert len(nsides) == 1
91  self._healpix = healpixes.pop()
92  self._nside = nsides.pop()
93  return self
94 
95  def read(self):
96  """Read the indices"""
97  if self._mi is not None:
98  return
99  fn = getIndexPath(self._filenameList[0])
100  if not os.path.exists(fn):
101  raise RuntimeError("Unable to get filename for astrometry star file %s" % (self._filenameList[0],))
102  self._mi = multiindex_new(fn)
103  if self._mi is None:
104  # Can't proceed at all without stars
105  raise RuntimeError('Failed to read stars from astrometry multiindex filename "%s"' % fn)
106  for i, fn in enumerate(self._filenameList[1:]):
107  if fn is None:
108  self.log.logdebug('Unable to find index part of multiindex %s' % fn)
109  continue
110  fn = getIndexPath(fn)
111  if not os.path.exists(fn):
112  self.log.warn("Unable to get filename for astrometry index %s" % (fn,))
113  continue
114  self.log.logdebug('Reading index from multiindex file "%s"' % fn)
115  if multiindex_add_index(self._mi, fn, INDEX_ONLY_LOAD_METADATA):
116  raise RuntimeError('Failed to read index from multiindex filename "%s"' % fn)
117  ind = self._mi[i]
118  self.log.logdebug(' index %i, hp %i (nside %i), nstars %i, nquads %i' %
119  (ind.indexid, ind.healpix, ind.hpnside, ind.nstars, ind.nquads))
120 
121  def reload(self):
122  """Reload the indices."""
123  if self._loaded:
124  return
125  if self._mi is None:
126  self.read()
127  else:
128  self._mi.reload()
129  self._loaded = True
130 
131  def unload(self):
132  """Unload the indices"""
133  if not self._loaded:
134  return
135  self._mi.unload()
136  self._loaded = False
137 
138  def isWithinRange(self, coord, distance):
139  """!Is the index within range of the provided coordinates?
140 
141  @param coord Coordinate to check (lsst.afw.coord.Coord)
142  @param distance Angular distance (lsst.afw.geom.Angle)
143  """
144  return (self._nside == 0 or healpixDistance(self._healpix, self._nside, coord) <= distance)
145 
146  def __getitem__(self, i):
147  self.reload()
148  return self._mi[i]
149 
150  def __len__(self):
151  return len(self._filenameList) - 1 # The first is the multiindex; the rest are the indices
152 
153  def __iter__(self):
154  self.reload()
155  return iter(self._mi)
156 
157 
158 class AstrometryNetCatalog(object):
159  """An interface to an astrometry.net catalog
160 
161  Behaves like a list of MultiIndexCache (or multiindex_t).
162 
163  These should usually be constructed using the 'fromEnvironment'
164  class method, which wraps the 'fromIndexFiles' and 'fromCache'
165  alternative class methods.
166  """
167  _cacheFilename = "andCache.fits"
168 
169  def __init__(self, andConfig):
170  """!Constructor
171 
172  @param andConfig Configuration (an AstrometryNetDataConfig)
173  """
174  self.config = andConfig
175  cacheName = getIndexPath(self._cacheFilename)
176  if self.config.allowCache and os.path.exists(cacheName):
177  self._initFromCache(cacheName)
178  else:
179  self._initFromIndexFiles(self.config)
180 
181  def _initFromIndexFiles(self, andConfig):
182  """Initialise from the index files in an AstrometryNetDataConfig"""
183  indexFiles = zip(andConfig.indexFiles, andConfig.indexFiles) + andConfig.multiIndexFiles
184  self._multiInds = [MultiIndexCache.fromFilenameList(fnList) for fnList in indexFiles]
185 
186  def writeCache(self):
187  """Write a cache file
188 
189  The cache file is a FITS file with all the required information to build the
190  AstrometryNetCatalog quickly. The first table extension contains a row for each multiindex,
191  storing the healpix and nside values. The second table extension contains a row
192  for each filename in all the multiindexes. The two may be JOINed through the
193  'id' column.
194  """
195  outName = getIndexPath(self._cacheFilename)
196  numFilenames = sum(len(ind._filenameList) for ind in self._multiInds)
197  maxLength = max(len(fn) for ind in self._multiInds for fn in ind._filenameList) + 1
198 
199  # First table
200  first = pyfits.new_table([pyfits.Column(name="id", format="K"),
201  pyfits.Column(name="healpix", format="K"),
202  pyfits.Column(name="nside", format="K"),
203  ], nrows=len(self._multiInds))
204  first.data.field("id")[:] = numpy.arange(len(self._multiInds), dtype=int)
205  first.data.field("healpix")[:] = numpy.array([ind._healpix for ind in self._multiInds])
206  first.data.field("nside")[:] = numpy.array([ind._nside for ind in self._multiInds])
207 
208  # Second table
209  second = pyfits.new_table([pyfits.Column(name="id", format="K"),
210  pyfits.Column(name="filename", format="%dA" % (maxLength)),
211  ], nrows=numFilenames)
212  ident = second.data.field("id")
213  filenames = second.data.field("filename")
214  i = 0
215  for j, ind in enumerate(self._multiInds):
216  for fn in ind._filenameList:
217  ident[i] = j
218  filenames[i] = fn
219  i += 1
220 
221  pyfits.HDUList([pyfits.PrimaryHDU(), first, second]).writeto(outName, clobber=True)
222 
223  def _initFromCache(self, filename):
224  """Initialise from a cache file
225 
226  Ingest the cache file written by the 'writeCache' method and
227  use that to quickly instantiate the AstrometryNetCatalog.
228  """
229  with pyfits.open(filename) as hduList:
230  first = hduList[1].data
231  second = hduList[2].data
232 
233  # first JOIN second USING(id)
234  filenames = {i: [] for i in first.field("id")}
235  for id2, fn in zip(second.field("id"), second.field("filename")):
236  filenames[id2].append(fn)
237  self._multiInds = [MultiIndexCache(filenames[i], hp, nside) for i, hp, nside in
238  zip(first.field("id"), first.field("healpix"), first.field("nside"))]
239 
240  # Check for consistency
241  cacheFiles = set(second.field("filename"))
242  configFiles = set(sum(self.config.multiIndexFiles, []) + self.config.indexFiles)
243  assert(cacheFiles == configFiles)
244 
245  def __getitem__(self, ii):
246  return self._multiInds[ii]
247 
248  def __iter__(self):
249  return iter(self._multiInds)
250 
251  def __len__(self):
252  return len(self._multiInds)
253 
254 
255 def generateCache(andConfig=None):
256  """Generate a cache file"""
257  if andConfig is None:
258  andConfig = getConfigFromEnvironment()
259  catalog = AstrometryNetCatalog(andConfig)
260  try:
261  for index in catalog:
262  index.reload()
263  catalog.writeCache()
264  finally:
265  for index in catalog:
266  index.unload()
int iter
def getIndexPath
Get the path to the specified astrometry.net index file.
Definition: multiindex.py:11
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
std::string getPackageDir(std::string const &packageName)
return the root directory of a setup package
Definition: Utils.cc:66
def isWithinRange
Is the index within range of the provided coordinates?
Definition: multiindex.py:138