diff --git a/CHANGELOG.md b/CHANGELOG.md index 2036982..99ee077 100644 --- a/CHANGELOG.md +++ b/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). diff --git a/pygdsm/__init__.py b/pygdsm/__init__.py index a2c37bc..e82a091 100644 --- a/pygdsm/__init__.py +++ b/pygdsm/__init__.py @@ -28,17 +28,20 @@ def init_gsm(gsm_name: str = "gsm08"): 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 @@ def init_observer(gsm_name: str = "gsm08"): 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}"') diff --git a/pygdsm/base_observer.py b/pygdsm/base_observer.py index 1386f82..56ad67a 100644 --- a/pygdsm/base_observer.py +++ b/pygdsm/base_observer.py @@ -42,11 +42,12 @@ 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): - """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 ---------- @@ -54,6 +55,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 ------- @@ -74,13 +77,35 @@ def generate(self, freq=None, obstime=None): 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 @@ def generate(self, freq=None, obstime=None): 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 diff --git a/setup.cfg b/setup.cfg index 8b09a37..09c4c7c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -7,6 +7,9 @@ universal=1 [metadata] description_file=README.md +[options] +python_requires = >=3.6 + [aliases] test=pytest diff --git a/setup.py b/setup.py index 5980f2e..34017ff 100644 --- a/setup.py +++ b/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", diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index 6d7f7bd..c4a4e4f 100644 --- a/tests/test_gsm2016.py +++ b/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) diff --git a/tests/test_gsm_observer.py b/tests/test_gsm_observer.py index 5ed762b..7b48c48 100644 --- a/tests/test_gsm_observer.py +++ b/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) diff --git a/tests/test_haslam.py b/tests/test_haslam.py index 12a28d8..2d75251 100644 --- a/tests/test_haslam.py +++ b/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) diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index bda238d..74c738f 100644 --- a/tests/test_lfsm.py +++ b/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)