Bump to 1.6.0

CHANGELOG.md

@@ -1,3 +1,6 @@
+* 1.6.0 (2024.12.21) - Moved data to datacentral.org.au, now downloads maps only as needed.
+                       BaseObserver.generate() now allow for the horizon to be set (thanks D. McKenna).
+                       Removed case statements to support older Python versions (thanks @ sjoerd-bouma)
 * 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).

pygdsm/__init__.py

@@ -28,17 +28,20 @@
         sky_model (various): Corresponding sky model
     """
     gsm_name = gsm_name.lower().strip()
-    match gsm_name:
-        case "gsm":  # Shorthand for GSM08
-            return GlobalSkyModel()
-        case "gsm08":
-            return GlobalSkyModel()
-        case "gsm16":
-            return GlobalSkyModel16()
-        case "lfsm":
-            return LowFrequencySkyModel()
-        case "haslam":
-            return HaslamSkyModel()
+
+    if gsm_name == 'gsm': # Shorthand for GSM08
+        return GlobalSkyModel()
+    elif gsm_name == 'gsm08':
+        return GlobalSkyModel()
+    elif gsm_name == 'gsm16':
+        return GlobalSkyModel16()
+    elif gsm_name == 'lfsm':
+        return LowFrequencySkyModel()
+    elif gsm_name == 'haslam':
+        return HaslamSkyModel()
+    else:
+        raise ValueError(f'Invalid model specification "{gsm_name}"')
+
 
 
 def init_observer(gsm_name: str = "gsm08"):
@@ -61,14 +64,16 @@
         observer (various): Corresponding sky model observer
     """
     gsm_name = gsm_name.lower().strip()
-    match gsm_name:
-        case "gsm":  # Shorthand for GSM08
-            return GSMObserver()
-        case "gsm08":
-            return GSMObserver()
-        case "gsm16":
-            return GSMObserver16()
-        case "lfsm":
-            return LFSMObserver()
-        case "haslam":
-            return HaslamObserver()
+
+    if gsm_name == 'gsm': # Shorthand for GSM08
+        return GSMObserver()
+    elif gsm_name == 'gsm08':
+        return GSMObserver()
+    elif gsm_name == 'gsm16':
+        return GSMObserver16()
+    elif gsm_name == 'lfsm':
+        return LFSMObserver()
+    elif gsm_name == 'haslam':
+        return HaslamObserver()
+    else:
+        raise ValueError(f'Invalid model specification "{gsm_name}"')

pygdsm/base_observer.py

@@ -42,18 +42,21 @@
 
         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):
-        """Generate the observed sky for the observer, based on the GSM.
+    def generate(self, freq=None, obstime=None, horizon_elevation=None):
+        """ Generate the observed sky for the observer, based on the GSM.
 
         Parameters
         ----------
         freq: float
             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
         -------
@@ -74,13 +77,35 @@
             time_has_changed = False
         else:
             time_has_changed = True
-            self._time = Time(
-                obstime
-            )  # This will catch datetimes, but Time() object should be passed
+            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)=}).")
+
+        # 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"))
@@ -92,8 +117,8 @@
             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 = np.ones(shape=self._n_pix, dtype='bool')
+            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
 

setup.cfg

@@ -7,6 +7,9 @@ universal=1
 [metadata]
 description_file=README.md
 
+[options]
+python_requires = >=3.6
+
 [aliases]
 test=pytest
 

setup.py

@@ -21,7 +21,7 @@
 
 setup(
     name="pygdsm",
-    version="1.5.5",
+    version="1.6.0",
     description="Python Global Sky Model of diffuse Galactic radio emission",
     long_description=long_description,
     long_description_content_type="text/markdown",
@@ -39,8 +39,12 @@
         "Topic :: Scientific/Engineering :: Astronomy",
         "License :: OSI Approved :: MIT License",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.10",
-        "Programming Language :: Python :: 3.11",
+        'Programming Language :: Python :: 3.6',
+        'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
+        'Programming Language :: Python :: 3.11',
     ],
     # What does your project relate to?
     keywords="radio astronomy sky model galactic diffuse emission",

tests/test_gsm2016.py

@@ -44,6 +44,28 @@ 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)

tests/test_gsm_observer.py

@@ -5,6 +5,7 @@
 import numpy as np
 import pylab as plt
 from astropy.time import Time
+import pytest
 
 from pygdsm import GSMObserver, GSMObserver16, LFSMObserver
 
@@ -44,6 +45,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()
 
@@ -58,11 +63,13 @@ def test_observed_mollview():
     ov.elev = elevation
 
     obs = []
-    if not os.path.exists("generated_sky"):
-        os.mkdir("generated_sky")
+    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()
@@ -79,12 +86,23 @@ 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()
+
+        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/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")
+    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')
+
 
 
 def test_generate_with_and_without_args():
@@ -109,7 +127,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

@@ -44,6 +44,28 @@ 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)

tests/test_lfsm.py

@@ -51,6 +51,28 @@ 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)