Coverage for /opt/anaconda/lib/python3.7/site-packages/llama/files/healpix/skymap.py : 46%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# (c) Stefan Countryman 2016-2018
3"""
4Define class containing a HEALPix skymap and convenience methods for
5manipulating, creating, and reformating skymaps.
6"""
8from abc import ABC, abstractproperty, abstractmethod
9from copy import deepcopy
10from deprecated import deprecated
11from llama.utils import memoize
12from llama.filehandler import TriggerList
13from llama.files.healpix.utils import (
14 dangle_rad,
15 sky_area,
16)
19@deprecated
20def psf_gaussian(delta_angle, sigma, degrees=True):
21 """Return a 2-D Gaussian PDF function to be used as a PSF (point spread
22 function) on the sphere. Note that this Gaussian PDF *does not* take
23 account of curvature and is accurate only for small angular spreads.
25 Parameters
26 ----------
27 delta_angle : array_like
28 A scalar or array of angular distances from the center of the Gaussian.
29 Must be in same units (degrees or radians) as sigma.
30 sigma : float
31 The standard deviation of the Gaussian. Must be in same units (degrees
32 or radians) as delta_angle.
33 degrees : bool, optional
34 Whether angles are given in degrees. If False, units are assumed to be
35 in radians.
37 Returns
38 -------
39 densities : array_like
40 The values of the Gaussian at the distances given by delta_angle in
41 units of [1/steradian] (even if inputs are given in degrees!).
43 >>> import numpy as np
44 >>> psf_gaussian(1, 1, degrees=False) == 1/(2*np.pi) * np.exp(-1/2)
45 True
46 """
47 import numpy as np
48 two_sigma2 = 2*sigma**2
49 densities = np.exp(-delta_angle**2/two_sigma2) / (np.pi*two_sigma2)
50 if degrees:
51 from astropy.units import Quantity
52 return Quantity(densities, "1/deg^2").to("1/steradian").value
53 return densities
56@deprecated
57def psf_gaussian_scales(sigma):
58 """Get the characteristic maximum and minimum scales of a Gaussian PSF for
59 a given angle sigma.
61 Parameters
62 ----------
63 sigma : int or float
64 The angular standard deviation of the Gaussian PSF in degrees or
65 radians.
67 Returns
68 -------
69 min_scale : float
70 The scale of inter-pixel distances required to capture the features of
71 a Gaussian with the given standard deviation. In same units as
72 ``sigma``.
73 max_scale : float
74 The diameter of a pixel mesh centered on the Gaussian mean that is
75 large enough to capture almost all of the probability in the
76 distribution. In same units as ``sigma``.
77 """
78 return (sigma / 10., sigma * 10)
81class SkyMap(ABC):
82 """An array representing a skymap."""
84 # pylint: disable=invalid-name
85 def __init__(self, pixels, ra, dec, degrees=True, area_per_pixel=None):
86 """Take an input array or list of pixels. ``pixels``, ``ra``, and
87 ``dec`` must have the same dimensions.
89 Parameters
90 ----------
91 pixels : array-like
92 Pixel values.
93 ra : array-like
94 Right Ascensions of the pixels.
95 dec : array-like
96 Declinations of the pixels.
97 degrees : bool, optional
98 Whether input angular positions are given in degrees. If ``False``,
99 assume radians.
100 area_per_pixel : astropy.Quantity, optional
101 Area per pixel. Can be a scalar value (if all pixels are the same)
102 or an array of pixel areas corresponding to the pixels given in
103 ``pixels``. If none is provided (DEFAULT), just divide the total
104 area of the sky by the number of pixels.
106 CURRENTLY ASSUMES EQUATORIAL COORDINATE SYSTEM.
107 """
108 import numpy as np
109 if not (len(pixels) == len(ra) and len(pixels) == len(dec)):
110 raise ValueError("ra, dec, and pixels must have same length.")
111 self.pixels = np.array(pixels)
112 if degrees:
113 self.ra = np.array(ra)
114 self.dec = np.array(dec)
115 else:
116 # these pylint disables are for spurious pylint errors from numpy;
117 # E1111=assignment-from-no-return
118 self.ra = np.degrees(ra) # pylint: disable=E1111
119 # pylint: disable=assignment-from-no-return
120 self.dec = np.degrees(dec)
121 if area_per_pixel is not None:
122 self.area_per_pixel = area_per_pixel
124 @classmethod
125 def _from_dict(cls, dictionary):
126 """Create a new instance from a dictionary of the sort returned by
127 ``to_dict``."""
128 pixels = dictionary.pop('pixels')
129 return cls(pixels, **dictionary)
131 def _to_dict(self):
132 """Return a dictionary whose keys are the attributes of this instance.
133 Can be used to initialize a new identical skymap instance."""
134 return {'pixels': self.pixels, 'ra': self.ra, 'dec': self.dec,
135 'area_per_pixel': self.area_per_pixel, 'degrees': True}
137 @property
138 def npix(self):
139 return len(self.pixels)
141 @abstractmethod
142 def dangle(self, ra, dec, degrees=True):
143 """Get a new skymap with the same pixelization where each pixel value
144 is the distance from that pixel to the point with right ascension
145 ``ra`` and declination ``dec``. By default, ``ra`` and ``dec`` are in
146 degrees, though this can be overridden by setting ``degrees=False``,
147 in which case they will be interpretted as radians. The measured angle
148 differences in the returned skymap will always be in degrees. Dot
149 product formula per-pixel is:
151 sin(DEC1)*sin(DEC2) + cos(DEC1)*cos(DEC2)*cos(RA2-RA1)
152 """
153 raise NotImplementedError()
154 # TODO Implement
156 def map(self, function):
157 """Map a function to the pixel locations of this skymap and return an
158 identically sized and ordered skymap with the output values of that
159 function as its pixel values. ``function`` must take vector arguments
160 ``ra`` and ``dec`` and return a list of pixel values of the same
161 length. No normalization takes place."""
162 result = self._to_dict()
163 result['pixels'] = function(self.ra, self.dec)
164 if len(result['pixels']) != self.npix:
165 raise ValueError('mapped function produced non-matching output '
166 'length.')
167 return self._from_dict(result)
169 @property
170 def ra(self):
171 """Right Ascensions in degrees."""
172 return self._ra
174 @ra.setter
175 def ra(self, ra):
176 """Set the Right Ascensions in degrees."""
177 self._ra = ra
179 @property
180 def dec(self):
181 """Declinations in degrees."""
182 return self._dec
184 @dec.setter
185 def dec(self, dec):
186 """Set the Declinations in degrees."""
187 self._dec = dec
189 @property
190 def area_per_pixel(self):
191 """Get the solid angle per pixel. Can be set to a
192 scalar value or a vector with the same length as the list of pixels; in
193 either case, the pixel areas are returned as a vector for each pixel.
194 If not set, assumes that all pixel areas are identical and that the
195 pixel areas sum to a full sphere."""
196 try:
197 return self._area_per_pixel
198 except AttributeError:
199 self.area_per_pixel = sky_area() / self.npix
200 return self.area_per_pixel
202 @area_per_pixel.setter
203 def area_per_pixel(self, area_per_pixel):
204 import numpy as np
205 from astropy.units import Quantity
206 if not (isinstance(area_per_pixel, Quantity) and
207 area_per_pixel.unit.physical_type == 'solid angle'):
208 raise ValueError("area_per_pixel must be an astropy.units.Quantity"
209 "of physical_type solid angle.")
210 # assign all pixels same value if given a scalar quantity
211 if area_per_pixel.shape == ():
212 self._area_per_pixel = np.ones_like(self.pixels)*area_per_pixel
213 elif area_per_pixel.shape == (self.npix,):
214 self._area_per_pixel = area_per_pixel
215 else:
216 raise ValueError("area_per_pixel must be scalar or a vector of the"
217 " same length as pixels.")
219 @memoize(method=True)
220 def integrate(self):
221 """If ``pixels`` is given in units of inverse solid angle, return the
222 integral of the pixel values over the whole sky. If ``pixels`` is
223 dimensionless, assume the pixel values are probability per-pixel and
224 simply return the sum of pixel values. Either way, the integral is
225 dimensionless and equals the total probability of the SkyMap (given the
226 stated assumptions)."""
227 from astropy.units import Unit
228 if self.unit == Unit(""):
229 return sum(self.pixels)
230 if (self.unit**-1).physical_type == 'solid angle':
231 return sum(self.pixels*self.area_per_pixel)
232 raise ValueError("Cannot integrate ``pixels`` with units of: "
233 "``{}``".format(self.unit))
235 @property
236 def unit(self):
237 """Get the units of the pixels of this skymap. If ``pixels`` is already
238 an ``astropy.Quantity``, simply return its ``unit``. Otherwise, return
239 ``astropy.units.Unit("")`` (i.e. dimensionless)."""
240 try:
241 return self.pixels.units
242 except AttributeError:
243 from astropy.units import Unit
244 return Unit("")
246 @memoize(method=True)
247 def normalize(self, norm_to=1):
248 """Scale the pixel values of this skymap so that they integrate to the
249 value specified by ``norm_to`` (DEFAULT: 1). Returns a new skymap with
250 the same units as ``self.pixels`` (assuming ``self.integrate`` is
251 defined for ``self.units`` used in ``self.pixels``."""
252 normalized = deepcopy(self)
253 normalized.pixels = (normalized.pixels/self.integrate())*norm_to
254 return normalized
256 @property
257 @memoize(method=True)
258 def as_pdf(self):
259 """Return a skymap whose pixels represent proability/pixel
260 (rather than probability per solid angle, i.e. a proper PDF) and
261 putting them into units of probability per solid angle. Returns a new
262 skymap with units of inverse solid angle. If this skymap is already a
263 PDF in units of inverse solid angle, this operation returns the
264 original skymap."""
265 from astropy.units import Unit
266 if self.unit == Unit(""):
267 pdf = deepcopy(self)
268 pdf.pixels = pdf.pixels / pdf.area_per_pixel
269 return pdf
270 if (self.unit**-1).physical_type == 'solid angle':
271 return self
272 raise ValueError("Cannot convert ``pixels`` with units of: "
273 "``{}`` to a PDF.".format(self.unit))
275 @memoize(method=True)
276 def confidence_region_size(self, percent=90):
277 """Calculate the size of the confidence region in square degrees. By
278 default, finds the 90% confidence region, though the percentage can be
279 manually specified. Does NOT normalize pixels before summing them (in
280 case the total probability does not add to 1, e.g. if we are already
281 dealing with a partial skymap)."""
282 import numpy as np
283 from astropy.units import Unit
284 pdf = self.as_pdf
285 sorted_inds = pdf.pixels.argsort()[::-1]
286 cumsum = (pdf.area_per_pixel*pdf.pixels)[sorted_inds].cumsum()
287 assert cumsum.unit == Unit("")
288 final_ind = np.nonzero(cumsum.value <= (percent/100.))[0][-1]
289 return sum(pdf.area_per_pixel[sorted_inds][:final_ind+1])
292class HEALPixRepresentableFileHandler(TriggerList):
293 """
294 A ``FileHandler`` class containing triggers that can be represented in
295 HEALPix. Provides ``get_healpix`` and ``get_healpix_descriptions`` methods
296 for this purpose. Implementations need not store data in a HEALPix format
297 natively, as long as they provide a ``get_healpix`` method for returning
298 valid ``HEALPixSkyMap`` instances.
299 """
301 @abstractmethod
302 def get_healpix(self, nest=None, template_skymap=None):
303 """Return a list of HEALPix skymaps representing the triggers in this
304 ``FileHandler``."""
306 def get_healpix_descriptions(self, skymap_indices=None):
307 """Get a list of descriptive strings for each HEALPix skymap contained
308 in this file.
310 Parameters
311 ----------
312 skymap_indices : list
313 A list of tuples describing the index of each skymap used to
314 produce this output skymap, with each entry in the tuple
315 corresponding to a detector. In this case, the descriptive strings
316 will be generated based on this ``FileHandler`` instance's
317 ``detector`` property and the indices of the skymap_indices in
318 ``skymap_indices``. The ordering of the returned descriptive
319 strings will match the ordering of skymap_indices. Note that for
320 single-detector skymaps, the elements of the ``skymap_indices``
321 list should just be single-element tuples corresponding to the
322 index of each skymap.
324 Returns
325 -------
326 descriptions : list
327 Descriptions of each HEALPix skymap in the order in which they
328 appear in the file.
329 """
330 if skymap_indices is None:
331 raise ValueError('Must provide skymap indices.')
332 descriptions = list()
333 for ituple in skymap_indices:
334 fmt = ('{}_' + '-{}_'.join(['{:03d}'.format(s) for s in ituple]))
335 descriptions.append(fmt.format(*self.DETECTORS))
336 return descriptions
339class HEALPixPSF(HEALPixRepresentableFileHandler):
340 """
341 A class with utility functions for generating skymaps in-memory that
342 follow the HEALPix parameters (NSIDE, ORDERING) of the file described by a
343 template ``HEALPixSkyMapFileHandler``. Used to dynamically generate skymaps
344 that are compatible with the template skymap based on some PSF (Gaussian by
345 default). Implementations of this class should contain point-like events
346 with some sort of point-spread function that can be used to generate
347 skymaps.
348 """
350 @abstractproperty
351 def template_skymap_filehandler(self):
352 """The ``HEALPixSkyMapFileHandler`` whose HEALPix parameters should be
353 used as a template for the ``HEALPixSkyMap`` output by this
354 ``FileHandler``. It is assumed that this ``HEALPixSkyMapFileHandler``
355 only returns one skymap in its list; consequently, only the first
356 skymap loaded by this file handler will be used. Note that the
357 template skymap file handler is likely not a dependency of this file
358 handler; it is only used as a default for formatting generated skymaps
359 from this file handler's data."""
361 @abstractmethod
362 def source_locations(self):
363 """Returns a list of tuples of ``(RA, Dec, sigma)`` for all point-like
364 sources, where ``(RA, Dec)`` is the central sky position and ``sigma``
365 is the standard deviation of the Gaussian PSF. Since this depends on
366 the structure of data in the relevant file handler, it must be
367 specified for each subclass."""
369 def get_healpix(self, nest=None, template_skymap=None):
370 """Generate a list of ``HEALPixSkyMap`` instances from the events
371 contained in this FileHandler.
373 Parameters
374 ----------
375 nest : bool, optional
376 If ``True``, return skymaps with NESTED ordering. If False, return
377 skymaps with RING ordering. Overrides the ordering (but not NSIDE)
378 defined in ``template_skymap_filehandler``.
379 template_skymap : HEALPixSkyMap, optional
380 Provide a ``HEALPixSkyMap`` as a template to override the ones
381 provided by ``template_skymap_filehandler``. Also overrides
382 ``nest`` (if provided).
384 Returns
385 -------
386 skymaps : list
387 A list of ``HEALPixSkyMap`` instances of the triggers in this
388 ``FileHandler`` that conform to the ``HEALPixSkyMap`` parameters in
389 this class's ``template_skymap_filehandler``
390 """
391 if template_skymap is None:
392 templatefh = self.template_skymap_filehandler(self)
393 template_skymap = templatefh.get_healpix()[0]
394 if nest is not None:
395 if nest:
396 template_skymap = template_skymap.to_nest()
397 else:
398 template_skymap = template_skymap.to_ring()
399 return [template_skymap.psf_gaussian(*p)
400 for p in self.source_locations()]
403class MollviewMixin(ABC):
404 """Defines an interface for plotting Mollweide projections."""
406 @abstractmethod
407 def mollview(self, **kwargs):
408 """Return a Mollweide projection of this skymap as a ``matplotlib``
409 figure."""
412class HEALPixSkyMap(SkyMap, MollviewMixin):
413 """
414 A HEALPix skymap object residing in memory, including methods for
415 working with the skymap and combining skymaps with each other.
416 """
418 def __init__(self, pixels, nest=True, title=None):
419 """Take an input array of pixel values and the HEALPix ORDER parameter,
420 "NEST" if ``nest=True`` (DEFAULT), "RING" if ``nest=False``.
422 Optionally provide a ``title`` for plots made from this skymap
423 describing what it shows.
425 CURRENTLY ASSUMES EQUATORIAL COORDINATE SYSTEM.
426 """
427 import numpy as np
428 self.nest = nest
429 self.pixels = np.array(pixels)
430 self.title = title
432 # pylint: disable=C0103
433 @classmethod
434 @memoize(method=True)
435 def nside2ang(cls, nside, nest=True, degrees=True):
436 """Return (RA, Dec) in degrees for the pixels corresponding to a
437 HEALPix skymap with ``NSIDE=nside`` where ``nest=True`` (DEFAULT)
438 indicates NEST ordering and ``nest=False`` indicates RING ordering. If
439 ``degrees=True`` (DEFAULT), the return value will be in degrees;
440 otherwise, radians.
442 Examples
443 --------
444 >>> import numpy as np
445 >>> ra, dec = HEALPixSkyMap.nside2ang(1, nest=True, degrees=True)
446 >>> all(1e-13 > ra - np.array([
447 ... 45., 135., 225., 315., 0., 90., 180., 270., 45., 135.,
448 ... 225., 315.
449 ... ]))
450 True
451 """
452 import numpy as np
453 if not degrees:
454 ra, dec = cls.nside2ang(nside, nest=nest, degrees=True)
455 return (np.radians(ra), np.radians(dec))
456 import healpy as hp
457 return hp.pix2ang(nside, range(hp.pixelfunc.nside2npix(nside)),
458 nest=nest, lonlat=True)
460 def ang2pix(self, ra, dec, degrees=True):
461 """Return the pixel index on this skymap at right ascensions ``ra`` and
462 declinations ``dec``."""
463 import numpy as np
464 import healpy as hp
465 if not degrees:
466 ra = dec * (180./np.pi)
467 dec = dec * (180./np.pi)
468 return hp.ang2pix(self.nside, ra, dec, nest=self.nest, lonlat=True)
470 def _to_dict(self):
471 """Return a dictionary whose keys are the attributes of this instance.
472 Can be used to initialize a new identical skymap instance."""
473 return dict(nest=self.nest, pixels=self.pixels)
475 @classmethod
476 def _read_from_hdf5_group(cls, group):
477 """Read a HEALPix skymap from the provided HDF5 group and return a new
478 instance containing the data in that HDF5 group.
480 >>> import tempfile, os, h5py
481 >>> f = tempfile.NamedTemporaryFile(suffix='.hdf5', delete=False)
482 >>> fname = f.name
483 >>> f.close()
484 >>> f = h5py.File(fname, 'w')
485 >>> g = f.create_group('foo')
486 >>> skymap = HEALPixSkyMap(range(12))
487 >>> skymap == HEALPixSkyMap._read_from_hdf5_group(
488 ... skymap._write_to_hdf5_group(g)
489 ... )
490 True
491 >>> os.remove(fname)
492 """
493 skymap_dict = dict()
494 for key in group:
495 skymap_dict[key] = group[key][()]
496 return cls._from_dict(skymap_dict)
498 def _write_to_hdf5_group(self, group, compressed=True):
499 """Write this HEALPix skymap to the provided HDF5 group. Optionally,
500 require compression by setting ``compressed`` to ``True`` (DEFAULT:
501 ``False``)."""
502 skymap_dict = self._to_dict()
503 for key, val in skymap_dict.items():
504 if compressed and key in ['pixels']:
505 group.create_dataset(key, data=val, compression='gzip')
506 else:
507 group.create_dataset(key, data=val)
508 return group
510 @property
511 def ra(self): # pylint: disable=C0103
512 """Right Ascensions in degrees."""
513 return self.nside2ang(self.nside, nest=self.nest)[0]
515 @property
516 def dec(self):
517 """Declinations in degrees."""
518 return self.nside2ang(self.nside, nest=self.nest)[1]
520 @property
521 def nside(self):
522 """The NSIDE HEALPix parameter indicating how many times each base
523 pixel has been subdivided in each of the 2 dimensions on the sphere."""
524 import healpy.pixelfunc as pxlfnc
525 return pxlfnc.npix2nside(self.npix)
527 def to_ring(self):
528 """Return a skymap with the same pixel values and RING ordering."""
529 if not self.nest:
530 return self
531 import healpy.pixelfunc as pxlfnc
532 return type(self)(
533 self.pixels[pxlfnc.nest2ring(self.nside, range(self.npix))],
534 nest=False
535 )
537 def to_nest(self):
538 """Return a skymap with the same pixel values and NESTED ordering."""
539 if self.nest:
540 return self
541 import healpy.pixelfunc as pxlfnc
542 return type(self)(
543 self.pixels[pxlfnc.ring2nest(self.nside, range(self.npix))],
544 nest=True
545 )
547 def conform_skymaps(self, *args):
548 """Make sure all skymaps contained in ``*args`` have the same ordering
549 (NESTED or RING) as the current skymap (``self``), reordering the
550 pixels if necessary, and returning a tuple of the (possibly reordered)
551 input skymaps."""
552 # make sure all sizes are the same
553 if len(set([skymap.npix for skymap in args])) != 1:
554 raise ValueError('Cannot conform skymaps of different sizes.')
555 if self.nest:
556 return tuple([skymap.to_nest() for skymap in args])
557 return tuple([skymap.to_ring() for skymap in args])
559 def dangle(self, ra, dec, degrees=True): # pylint: disable=C0103
560 """Get a new skymap with the same pixelization where each pixel value
561 is the distance from that pixel to the point with right ascension
562 ``ra`` and declination ``dec``. By default, ``ra`` and ``dec`` are in
563 degrees, though this can be overridden by setting ``degrees=False``,
564 in which case they will be interpretted as radians. The measured angle
565 differences in the returned skymap will always be in degrees. Dot
566 product formula per-pixel is:
568 sin(DEC1)*sin(DEC2) + cos(DEC1)*cos(DEC2)*cos(RA2-RA1)
570 So we get the angle by taking the arccos of this function.
572 Examples
573 --------
574 We should find that the distance from any pixel to the North pole is
575 equal to 90 minus the declination (within some small error):
577 >>> skymap = HEALPixSkyMap([1]*12)
578 >>> skymap.dangle(32, 90).pixels + skymap.dec - 90 < 1e-12
579 array([ True, True, True, True, True, True, True, True, True,
580 True, True, True])
582 Likewise, the distance from any pixel to the South pole should be
583 equal to 90 plus the declination:
585 >>> skymap.dangle(359, -90).pixels - skymap.dec - 90. < 1e-12
586 array([ True, True, True, True, True, True, True, True, True,
587 True, True, True])
588 """
589 import numpy as np
590 if degrees:
591 ra = np.radians(ra) # pylint: disable=assignment-from-no-return
592 dec = np.radians(dec) # pylint: disable=assignment-from-no-return
593 mapra = self.nside2ang(self.nside, nest=self.nest, degrees=False)[0]
594 mapdec = self.nside2ang(self.nside, nest=self.nest, degrees=False)[1]
595 return type(self)(np.degrees(dangle_rad(ra, dec, mapra, mapdec)),
596 nest=self.nest)
598 # pylint: disable=C0103,R0913
599 def psf_gaussian(self, ra, dec, sigma, degrees=True, normalize=True):
600 """Calculate a Gaussian PSF at all pixels in this skymap.
602 Parameters
603 ----------
604 ra : float
605 Right Ascension of the Gaussian mean.
606 dec : float
607 Declination of the Gaussian mean.
608 sigma : float
609 Standard deviation of the Gaussian PSF.
610 degrees : bool, optional
611 If true, input arguments (``ra``, ``dec``, ``sigma``) are assumed
612 to be in degrees. Otherwise, they are assumed to be in radians.
613 normalize : bool, optional
614 If true, normalize the returned skymap so that the integral over
615 area is 1 (useful for when the Gaussian is meant to represent a
616 probability density function).
618 Returns
619 -------
620 skymap : HEALPixSkyMap
621 A new skymap with the same pixelization where a 2-D Gaussian PSF
622 (point spread function) centered at ``(ra, dec)`` with standard
623 deviation ``sigma`` has been calculated for each pixel in this
624 skymap.
625 """
626 import numpy as np
627 if not degrees:
628 # pylint: disable=assignment-from-no-return
629 sigma = np.degrees(sigma)
630 min_scale, max_scale = psf_gaussian_scales(sigma)
631 skymap = type(self)(
632 psf_gaussian(self.dangle(ra, dec, degrees=degrees).pixels, sigma,
633 degrees=degrees),
634 nest=self.nest
635 )
636 if normalize:
637 skymap = skymap.normalize()
638 return skymap
640 def __eq__(self, other):
641 """Check whether two skymaps are the same WITHOUT conforming them
642 first. This means that two instances representing the SAME skymap with
643 different ordering will not be equal."""
644 if self.npix != other.npix:
645 return False
646 if self.nest != other.nest:
647 return False
648 return all(self.pixels == other.pixels)
650 def __add__(self, other):
651 return type(self)(self.pixels + self.conform_skymaps(other)[0].pixels,
652 nest=self.nest)
654 def __sub__(self, other):
655 return type(self)(self.pixels - self.conform_skymaps(other)[0].pixels,
656 nest=self.nest)
658 def __mul__(self, other):
659 return type(self)(self.pixels * self.conform_skymaps(other)[0].pixels,
660 nest=self.nest)
662 def __div__(self, other):
663 return type(self)(self.pixels / self.conform_skymaps(other)[0].pixels,
664 nest=self.nest)
666 def __pow__(self, power):
667 return type(self)(self.pixels ** power, nest=self.nest)
669 def __len__(self):
670 return self.npix
672 def __getitem__(self, key):
673 return self.pixels[key]
675 def __setitem__(self, key, value):
676 self.pixels[key] = value
678 def __repr__(self):
679 fmt = '{}({}, nest={})'
680 dic = self._to_dict()
681 return fmt.format(type(self).__name__, dic['pixels'], dic['nest'])
683 def mollview(self, **kwargs):
684 """Return a Mollweide projection of this skymap as a ``matplotlib``
685 figure with the correct nesting. Remaining ``**kwargs`` are passed to
686 ``healpy.mollview``. Use the LIGO standard rotation where 180 degrees
687 is at the center of the map."""
688 import healpy as hp
689 import matplotlib.pyplot as plt
690 if self.title is not None and "title" not in kwargs:
691 kwargs['title'] = self.title
692 kwargs['rot'] = kwargs.get('rot', (180, 0, 0))
693 hp.mollview(self.pixels, nest=self.nest, **kwargs)
694 hp.graticule(dpar=15, dmer=30)
695 return plt.gcf()