Merge pull request #23 from telegraphic/hpix_to_sky_utils

Adding observed_gsm property and hpix2sky, sky2hpix

CHANGELOG.md

@@ -1,5 +1,9 @@
+* 1.5.4 (2024.04.15) - Added observed_gsm property to `BaseObserver`,
+                       Added `hpix2sky` and `sky2hpix` helpers,
+                       Now using `query_disc` to find pixels below horizon (neater code).
 * 1.5.3 (2024.04.15) - Added `init_gsm()` and `init_observer()` helper functions
 * 1.5.2 (2024.03.29) - Fix PyPi distribution files (thanks to D. McKenna)
+
 * 1.5.0 (2024.03.29) - Fix to `GlobalSkyModel16`, correct T_CMB monopole addition (thanks to R. Francis)
 * 1.4.1 (2023.04.28) - Added `include_cmb` option to sky models.
 * 1.4.0 (2023.04.27) - Added new interpolation method (courtesy R. Braun)

pygdsm/base_observer.py

@@ -3,7 +3,8 @@
 import numpy as np
 from astropy.time import Time
 from pygdsm.plot_utils import show_plt
-
+from pygdsm.utils import hpix2sky, sky2hpix
+from astropy.coordinates import SkyCoord
 
 class BaseObserver(ephem.Observer):
     """ Observer of the Global Sky Model.
@@ -59,13 +60,13 @@ def generate(self, freq=None, obstime=None):
             with below the horizon masked.
         """
         # Check to see if frequency has changed.
-        if freq is not None: 
+        if freq is not None:
             if not np.isclose(freq, self._freq):
                 self.gsm.generate(freq)
                 self._freq = freq
-        
+
         sky = self.gsm.generated_map_data
-        
+
         # Check if time has changed -- astropy allows None == Time() comparison
         if obstime == self._time or obstime is None:
             time_has_changed = False
@@ -78,9 +79,20 @@ def generate(self, freq=None, obstime=None):
         if time_has_changed or self.observed_sky is None:
             # Get RA and DEC of zenith
             ra_zen, dec_zen = self.radec_of(0, np.pi/2)
+            sc_zen = SkyCoord(ra_zen, dec_zen, unit=('rad', 'rad'))
+            pix_zen = sky2hpix(self._n_side, sc_zen)
+            vec_zen = hp.pix2vec(self._n_side, pix_zen)
+
+            # Convert to degrees
             ra_zen  *= (180 / np.pi)
             dec_zen *= (180 / np.pi)
 
+            # Generate below-horizon mask using query_disc
+            mask = np.ones(shape=self._n_pix, dtype='bool')
+            pix_visible = hp.query_disc(self._n_side, vec=vec_zen, radius=np.pi/2)
+            mask[pix_visible] = 0
+            self._mask = mask
+
             # Transform from Galactic coordinates to Equatorial
             rot = hp.Rotator(coord=['G', 'C'])
             eq_theta, eq_phi = rot(self._theta, self._phi)
@@ -89,12 +101,6 @@ def generate(self, freq=None, obstime=None):
             dec = 90.0 - np.abs(eq_theta*(180/np.pi))
             ra = ( (eq_phi + 2*np.pi) % (2*np.pi) )*(180/np.pi)
 
-            # Generate a mask for below horizon (distance from zenith > 90 deg)
-            dist = 2 * np.arcsin( np.sqrt( np.sin((dec_zen-dec)*(np.pi/180)/2)**2 + \
-                    np.cos(dec_zen*(np.pi/180))*np.cos(dec*(np.pi/180))*np.sin((ra_zen-ra)*(np.pi/180)/2)**2 ) )
-            mask = dist > (np.pi/2)
-            self._mask = mask
-
             # Apply rotation to convert from Galactic to Equatorial and center on zenith
             hrot = hp.Rotator(rot=[ra_zen, dec_zen], coord=['G', 'C'], inv=True)
             g0, g1 = hrot(self._theta, self._phi)
@@ -105,7 +111,7 @@ def generate(self, freq=None, obstime=None):
             ra_rotated = ra[self._pix0]
             self._observed_ra = ra_rotated
             self._observed_dec = dec_rotated
- 
+
         sky_rotated = sky[self._pix0]
         mask_rotated = self._mask[self._pix0]
 
@@ -132,11 +138,10 @@ def view(self, logged=False, show=False, **kwargs):
             show_plt()
         return sky
 
-    def view_observed_gsm(self, logged=False, show=False, **kwargs):
-        """ View the GSM (Mollweide), with below-horizon area masked. """
+    @property
+    def observed_gsm(self):
+        """ Return the GSM (Mollweide), with below-horizon area masked. """
         sky = self.observed_sky
-        if logged:
-            sky = np.log2(sky)
 
         # Get RA and DEC of zenith
         ra_rad, dec_rad = self.radec_of(0, np.pi / 2)
@@ -153,10 +158,24 @@ def view_observed_gsm(self, logged=False, show=False, **kwargs):
         g0, g1 = coordrotate(self._theta, self._phi)
         pix0 = hp.ang2pix(self._n_side, g0, g1)
         sky = sky[pix0]
+        return sky
 
-        hp.mollview(sky, coord='G', **kwargs)
+    def view_observed_gsm(self, logged=False, show=False, **kwargs):
+        """ View the GSM (Mollweide), with below-horizon area masked.
+
+        Args:
+            logged (bool): Apply log2 to data (default False)
+            show (bool): Call plt.show() (default False)
+
+        Returns:
+            sky (np.array): Healpix map of observed GSM.
+        """
+        sky = self.observed_gsm
 
+        if logged:
+            sky = np.log2(sky)
+
+        hp.mollview(sky, coord='G', **kwargs)
         if show:
             show_plt()
-
         return sky

pygdsm/plot_utils.py

@@ -1,4 +1,5 @@
 def show_plt(): # pragma: no cover
+    """ Call plt.show() """
     try:
         plt.show() # noqa: F821
     except: