From bca0bcd8d7aa01e3b6a22a6b28ec8ed3c0d45827 Mon Sep 17 00:00:00 2001 From: David McKenna Date: Tue, 30 Apr 2024 16:50:09 +0800 Subject: [PATCH 1/4] Allow for an artificial horizon when generating an *Observer sky map. --- pygdsm/base_observer.py | 21 ++++++++++++++++++--- tests/test_gsm2016.py | 22 ++++++++++++++++++++++ tests/test_gsm_observer.py | 19 ++++++++++++++++++- tests/test_lfsm.py | 22 ++++++++++++++++++++++ 4 files changed, 80 insertions(+), 4 deletions(-) diff --git a/pygdsm/base_observer.py b/pygdsm/base_observer.py index 8013aa4..a5bfb37 100644 --- a/pygdsm/base_observer.py +++ b/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 horionz (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 diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index ca9cf70..caf65ef 100644 --- a/tests/test_gsm2016.py +++ b/tests/test_gsm2016.py @@ -45,6 +45,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) + + ov.generate(1400, horizon_elevation = '20.0') + d_20deg_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(20.0)) + d_20deg2rad_horizon = ov.view(logged=True) + + assert np.all(d_20deg_horizon == d_20deg2rad_horizon), "The two methods for calculating the artifial horizon do not match." + assert np.ma.count_masked(d_20deg_horizon).sum() == 8443259 # This is higher than LFSM and Haslam due to the nside=1024 sampling for the 'hi' resolution default + plt.show() + def test_interp(): f = np.arange(40, 80, 5) for interp in ('pchip', 'cubic'): diff --git a/tests/test_gsm_observer.py b/tests/test_gsm_observer.py index 4382d42..6e98270 100644 --- a/tests/test_gsm_observer.py +++ b/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 + elevation = '20.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=elevation) + ov.view(logged=True, show=False, min=9, max=20) + plt.savefig('generated_sky/ortho_20deg_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_20deg_horizon-*.png ortho_20deg_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(10)) + ov.generate(obstime=now, freq=52, horizon_elevation='10.0') if __name__ == "__main__": test_gsm_observer(show=True) diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index 1735c27..2b7488f 100644 --- a/tests/test_lfsm.py +++ b/tests/test_lfsm.py @@ -44,6 +44,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) + + ov.generate(200, horizon_elevation = '20.0') + d_20deg_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(20.0)) + d_20deg2rad_horizon = ov.view(logged=True) + + assert np.all(d_20deg_horizon == d_20deg2rad_horizon), "The two methods for calculating the artifial horizon do not match." + assert np.ma.count_masked(d_20deg_horizon).sum() == 527705 + plt.show() + def test_cmb_removal(): g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=False) sky_no_cmb = g.generate(400) From 0def15e7d8c31c14c53f6868eea8e4b79050eeb0 Mon Sep 17 00:00:00 2001 From: David McKenna Date: Mon, 6 May 2024 14:03:22 +0800 Subject: [PATCH 2/4] Fix "horionz" typo. --- pygdsm/base_observer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygdsm/base_observer.py b/pygdsm/base_observer.py index a5bfb37..fa422b1 100644 --- a/pygdsm/base_observer.py +++ b/pygdsm/base_observer.py @@ -54,7 +54,7 @@ def generate(self, freq=None, obstime=None, horizon_elevation=None): obstime: astropy.time.Time Time of observation to generate horizon_elevation: float - Elevation of the artificial horionz (default 0.0) + Elevation of the artificial horizon (default 0.0) Returns ------- From 5e9d51539f6d2b3111b2cb14a8e8d4cb5c6471a6 Mon Sep 17 00:00:00 2001 From: David McKenna Date: Thu, 9 May 2024 11:59:51 +0800 Subject: [PATCH 3/4] Increase the artificial horizon elevation in the tests to 85 deg. I seem to have run into a case of -ffast-math optimizations differing by platform, resulting in the tests returning different values on my ARM64 laptop, vs the x86_64 runners. Increasing the elevation reduces the effective number of boundary pixels, reducing the chance of the optimizations resulting in a variable number of pixels being masked by platform. --- tests/test_gsm2016.py | 13 +++++++------ tests/test_gsm_observer.py | 12 ++++++------ tests/test_haslam.py | 23 +++++++++++++++++++++++ tests/test_lfsm.py | 13 +++++++------ 4 files changed, 43 insertions(+), 18 deletions(-) diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index caf65ef..abf7330 100644 --- a/tests/test_gsm2016.py +++ b/tests/test_gsm2016.py @@ -51,8 +51,9 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(1400, horizon_elevation = '20.0') - d_20deg_horizon = ov.view(logged=True) + horizon_elevation = 85.0 + ov.generate(1400, horizon_elevation = str(horizon_elevation)) + d_85deg_horizon = ov.view(logged=True) ov = GSMObserver16() ov.lon = longitude @@ -60,11 +61,11 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(1400, horizon_elevation = np.deg2rad(20.0)) - d_20deg2rad_horizon = ov.view(logged=True) + ov.generate(1400, horizon_elevation = np.deg2rad(horizon_elevation)) + d_85deg2rad_horizon = ov.view(logged=True) - assert np.all(d_20deg_horizon == d_20deg2rad_horizon), "The two methods for calculating the artifial horizon do not match." - assert np.ma.count_masked(d_20deg_horizon).sum() == 8443259 # This is higher than LFSM and Haslam due to the nside=1024 sampling for the 'hi' resolution default + 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(): diff --git a/tests/test_gsm_observer.py b/tests/test_gsm_observer.py index 6e98270..482d9a6 100644 --- a/tests/test_gsm_observer.py +++ b/tests/test_gsm_observer.py @@ -63,7 +63,7 @@ def test_observed_mollview(): if not os.path.exists('generated_sky'): os.mkdir('generated_sky') freq = 50 - elevation = '20.0' + horizon_elevation = '85.0' for ii in range(0, 24, 4): ov.date = datetime(2000, 1, 1, ii, 0) ov.generate(freq) @@ -83,15 +83,15 @@ def test_observed_mollview(): plt.savefig('generated_sky/ortho-%02d.png' % ii) plt.close() - ov.generate(freq=freq, horizon_elevation=elevation) + ov.generate(freq=freq, horizon_elevation=horizon_elevation) ov.view(logged=True, show=False, min=9, max=20) - plt.savefig('generated_sky/ortho_20deg_horizon-%02d.png' % ii) + 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_20deg_horizon-*.png ortho_20deg_horizon.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') @@ -121,8 +121,8 @@ def test_generate_with_and_without_args(): 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(10)) - ov.generate(obstime=now, freq=52, horizon_elevation='10.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 a8aaf26..33d1987 100644 --- a/tests/test_haslam.py +++ b/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) diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index 2b7488f..4fcfc3b 100644 --- a/tests/test_lfsm.py +++ b/tests/test_lfsm.py @@ -50,8 +50,9 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(200, horizon_elevation = '20.0') - d_20deg_horizon = ov.view(logged=True) + horizon_elevation = 85.0 + ov.generate(200, horizon_elevation = str(horizon_elevation)) + d_85deg_horizon = ov.view(logged=True) ov = LFSMObserver() ov.lon = longitude @@ -59,11 +60,11 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(200, horizon_elevation = np.deg2rad(20.0)) - d_20deg2rad_horizon = ov.view(logged=True) + ov.generate(200, horizon_elevation = np.deg2rad(horizon_elevation)) + d_85deg2rad_horizon = ov.view(logged=True) - assert np.all(d_20deg_horizon == d_20deg2rad_horizon), "The two methods for calculating the artifial horizon do not match." - assert np.ma.count_masked(d_20deg_horizon).sum() == 527705 + 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(): From afaab59f5289c9a4eb62afc61962da34e08505e1 Mon Sep 17 00:00:00 2001 From: David McKenna Date: Thu, 9 May 2024 12:03:10 +0800 Subject: [PATCH 4/4] Match existing coding style a bit closer. --- tests/test_gsm2016.py | 4 ++-- tests/test_haslam.py | 4 ++-- tests/test_lfsm.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index abf7330..ec314e6 100644 --- a/tests/test_gsm2016.py +++ b/tests/test_gsm2016.py @@ -52,7 +52,7 @@ def test_observer_test(): ov.date = datetime(2000, 1, 1, 23, 0) horizon_elevation = 85.0 - ov.generate(1400, horizon_elevation = str(horizon_elevation)) + ov.generate(1400, horizon_elevation=str(horizon_elevation)) d_85deg_horizon = ov.view(logged=True) ov = GSMObserver16() @@ -61,7 +61,7 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(1400, horizon_elevation = np.deg2rad(horizon_elevation)) + 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." diff --git a/tests/test_haslam.py b/tests/test_haslam.py index 33d1987..23c00c6 100644 --- a/tests/test_haslam.py +++ b/tests/test_haslam.py @@ -49,7 +49,7 @@ def test_observer_test(): ov.date = datetime(2000, 1, 1, 23, 0) horizon_elevation = 85.0 - ov.generate(1400, horizon_elevation = str(horizon_elevation)) + ov.generate(1400, horizon_elevation=str(horizon_elevation)) d_85deg_horizon = ov.view(logged=True) ov = GSMObserver() @@ -58,7 +58,7 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(1400, horizon_elevation = np.deg2rad(horizon_elevation)) + 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." diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index 4fcfc3b..cd9e0a1 100644 --- a/tests/test_lfsm.py +++ b/tests/test_lfsm.py @@ -51,7 +51,7 @@ def test_observer_test(): ov.date = datetime(2000, 1, 1, 23, 0) horizon_elevation = 85.0 - ov.generate(200, horizon_elevation = str(horizon_elevation)) + ov.generate(200, horizon_elevation=str(horizon_elevation)) d_85deg_horizon = ov.view(logged=True) ov = LFSMObserver() @@ -60,7 +60,7 @@ def test_observer_test(): ov.elev = elevation ov.date = datetime(2000, 1, 1, 23, 0) - ov.generate(200, horizon_elevation = np.deg2rad(horizon_elevation)) + 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."