Allow for an artificial horizon when generating an *Observer sky map.

pygdsm/base_observer.py

@@ -40,10 +40,11 @@ def _setup(self):
 
         self._pix0 = None
         self._mask = None
+        self._horizon_elevation = 0.0
         self._observed_ra = None
         self._observed_dec = None
 
-    def generate(self, freq=None, obstime=None):
+    def generate(self, freq=None, obstime=None, horizon_elevation=None):
         """ Generate the observed sky for the observer, based on the GSM.
 
         Parameters
@@ -52,6 +53,8 @@ def generate(self, freq=None, obstime=None):
             Frequency of map to generate, in units of MHz (default).
         obstime: astropy.time.Time
             Time of observation to generate
+        horizon_elevation: float
+            Elevation of the artificial horizon (default 0.0)
 
         Returns
         -------
@@ -75,8 +78,20 @@ def generate(self, freq=None, obstime=None):
             self._time = Time(obstime)  # This will catch datetimes, but Time() object should be passed
             self.date  = obstime.to_datetime()
 
+        # Match pyephem convertion -- string is degrees, int/float is rad
+        horizon_elevation = ephem.degrees(horizon_elevation or 0.0)
+        if self._horizon_elevation == horizon_elevation:
+            horizon_has_changed = False
+        else:
+            self._horizon_elevation = horizon_elevation
+            horizon_has_changed = True
+
+        # Checking this separately encase someone tries to be smart and sets both the attribute and kwarg
+        if self._horizon_elevation < 0:
+                raise ValueError(f"Horizon elevation must be greater or equal to 0 degrees ({np.rad2deg(horizon_elevation)=}).")
+
         # Rotation is quite slow, only recompute if time or frequency has changed, or it has never been run
-        if time_has_changed or self.observed_sky is None:
+        if time_has_changed or self.observed_sky is None or horizon_has_changed:
             # 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'))
@@ -89,7 +104,7 @@ def generate(self, freq=None, obstime=None):
 
             # 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)
+            pix_visible = hp.query_disc(self._n_side, vec=vec_zen, radius=np.pi/2 - self._horizon_elevation)
             mask[pix_visible] = 0
             self._mask = mask
 

tests/test_gsm2016.py

@@ -45,6 +45,29 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+    ov = GSMObserver16()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    horizon_elevation = 85.0
+    ov.generate(1400, horizon_elevation=str(horizon_elevation))
+    d_85deg_horizon = ov.view(logged=True)
+
+    ov = GSMObserver16()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    ov.generate(1400, horizon_elevation=np.deg2rad(horizon_elevation))
+    d_85deg2rad_horizon = ov.view(logged=True)
+
+    assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match."
+    assert np.ma.count_masked(d_85deg_horizon).sum() == 12558953
+    plt.show()
+
 def test_interp():
     f = np.arange(40, 80, 5)
     for interp in ('pchip', 'cubic'):

tests/test_gsm_observer.py

@@ -5,6 +5,7 @@
 import numpy as np
 import os
 from astropy.time import Time
+import pytest
 
 def test_gsm_observer(show=False):
     """ Test GSMObserver() is working
@@ -42,6 +43,10 @@ def test_gsm_observer(show=False):
     ov.generate(50)
     ov.view(logged=True)
     ov.view_observed_gsm(logged=True)
+
+    with pytest.raises(ValueError) as e:
+        ov.generate(horizon_elevation=-1e-3)
+
     if show:
         plt.show()
 
@@ -57,9 +62,11 @@ def test_observed_mollview():
     obs = []
     if not os.path.exists('generated_sky'):
         os.mkdir('generated_sky')
+    freq = 50
+    horizon_elevation = '85.0'
     for ii in range(0, 24, 4):
         ov.date = datetime(2000, 1, 1, ii, 0)
-        ov.generate(50)
+        ov.generate(freq)
         sky = ov.view_observed_gsm(logged=True, show=False, min=9, max=20)
         plt.savefig('generated_sky/galactic-%02d.png' % ii)
         plt.close()
@@ -76,9 +83,15 @@ def test_observed_mollview():
         plt.savefig('generated_sky/ortho-%02d.png' % ii)
         plt.close()
 
+        ov.generate(freq=freq, horizon_elevation=horizon_elevation)
+        ov.view(logged=True, show=False, min=9, max=20)
+        plt.savefig('generated_sky/ortho_85deg_horizon-%02d.png' % ii)
+        plt.close()
+
         print(ii)
 
     os.system('convert -delay 20 generated_sky/ortho-*.png ortho.gif')
+    os.system('convert -delay 20 generated_sky/ortho_85deg_horizon-*.png ortho_85deg_horizon.gif')
     os.system('convert -delay 20 generated_sky/galactic-*.png galactic.gif')
     os.system('convert -delay 20 generated_sky/ecliptic-*.png ecliptic.gif')
     os.system('convert -delay 20 generated_sky/equatorial-*.png equatorial.gif')
@@ -106,6 +119,10 @@ def test_generate_with_and_without_args():
     ov.generate(obstime=now)
     ov.generate(obstime=now, freq=53)
     ov.generate(obstime=now, freq=52)
+    ov.generate(obstime=now, freq=53, horizon_elevation=0.0)
+    ov.generate(obstime=now, freq=52, horizon_elevation='0.0')
+    ov.generate(obstime=now, freq=53, horizon_elevation=np.deg2rad(85.0))
+    ov.generate(obstime=now, freq=52, horizon_elevation='85.0')
 
 if __name__ == "__main__":
     test_gsm_observer(show=True)

tests/test_haslam.py

@@ -42,6 +42,29 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+    ov = GSMObserver()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    horizon_elevation = 85.0
+    ov.generate(1400, horizon_elevation=str(horizon_elevation))
+    d_85deg_horizon = ov.view(logged=True)
+
+    ov = GSMObserver()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    ov.generate(1400, horizon_elevation=np.deg2rad(horizon_elevation))
+    d_85deg2rad_horizon = ov.view(logged=True)
+
+    assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match."
+    assert np.ma.count_masked(d_85deg_horizon).sum() == 3139749
+    plt.show()
+
 def test_cmb_removal():
     g = HaslamSkyModel(freq_unit='MHz', include_cmb=False)
     sky_no_cmb = g.generate(400)

tests/test_lfsm.py

@@ -44,6 +44,29 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+    ov = LFSMObserver()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    horizon_elevation = 85.0
+    ov.generate(200, horizon_elevation=str(horizon_elevation))
+    d_85deg_horizon = ov.view(logged=True)
+
+    ov = LFSMObserver()
+    ov.lon = longitude
+    ov.lat = latitude
+    ov.elev = elevation
+    ov.date = datetime(2000, 1, 1, 23, 0)
+
+    ov.generate(200, horizon_elevation=np.deg2rad(horizon_elevation))
+    d_85deg2rad_horizon = ov.view(logged=True)
+
+    assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match."
+    assert np.ma.count_masked(d_85deg_horizon).sum() == 784951
+    plt.show()
+
 def test_cmb_removal():
     g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=False)
     sky_no_cmb = g.generate(400)