diff --git a/README.md b/README.md index 6073306..3798789 100644 --- a/README.md +++ b/README.md @@ -1,29 +1,30 @@ [![ascl:1603.013](https://img.shields.io/badge/ascl-1603.013-blue.svg)](https://ascl.net/1603.013) [![codecov](https://codecov.io/gh/telegraphic/pygdsm/branch/master/graph/badge.svg)](https://codecov.io/gh/telegraphic/pygdsm) -[![astropy](http://img.shields.io/badge/powered%20by-AstroPy-orange.svg?style=flat)](http://www.astropy.org/) - +[![astropy](http://img.shields.io/badge/powered%20by-AstroPy-orange.svg?style=flat)](http://www.astropy.org/) + PyGDSM ===== ![skymodels.jpg](https://github.com/telegraphic/pygdsm/raw/master/docs/skymodels.jpg) - `PyGDSM` is a Python interface for global diffuse sky models: all-sky maps in Healpix format of diffuse Galactic radio emission. This package includes interfaces to: - * **GSM2008:** A model of diffuse Galactic radio emission from 10 MHz to 100 GHz, [Oliveira-Costa et. al., (2008)](https://ui.adsabs.harvard.edu/abs/2008MNRAS.388..247D/abstract). + * **GSM2008:** A model of diffuse Galactic radio emission from 10 MHz to 100 GHz, [Oliveira-Costa et. al., (2008)](https://ui.adsabs.harvard.edu/abs/2008MNRAS.388..247D/abstract). * **GSM2016:** An improved model of diffuse galactic radio emission from 10 MHz to 5 THz, [Zheng et. al., (2016)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.464.3486Z/abstract). - * **LFSS:** The LWA1 Low Frequency Sky Survey (10-408 MHz) [Dowell et. al. (2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.469.4537D/abstract). + * **LFSS:** The LWA1 Low Frequency Sky Survey (10-408 MHz) [Dowell et. al. (2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.469.4537D/abstract). * **Haslam:** A frequency-scaled model using a spectral index, based on the [Haslam 408 MHz](https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_info.cfm) all-sky map. -In general, these are *not* wrappers of the original code (GSM2008 was written in Fortan and GSM2016 in C); instead they provides a uniform API with some additional features and advantages, such as healpy integration for imaging, and sky rotation for observed skies. +In general, these are *not* wrappers of the original code (GSM2008 was written in Fortan and GSM2016 in C); instead they provides a uniform API with some additional features and advantages, such as healpy integration for imaging, and sky rotation for observed skies. +> [!TIP] +> We thank [datacentral.org.au](https://datacentral.org.au) for hosting the component map data. Quickstart ---------- -The first thing to do will be to make sure you've got the dependencies: +The first thing to do will be to make sure you've got the dependencies: * [numpy](http://www.numpy.org/) * [scipy](http://www.scipy.org/install.html) @@ -38,13 +39,13 @@ Then you should be able to install with: Alternatively, clone the directory: git clone https://github.com/telegraphic/pygdsm - -An run `pip install .`. On first run, the sky model data will be downloaded from Zenodo / LAMBDA. These are about 500 MB total, and will be downloaded into your astropy cache (`~/.astropy/`). The data are hosted on [Zenodo](https://zenodo.org/record/3479985#.XaASx79S-AY). + +An run `pip install .`. On first run, the sky model data will be downloaded from Zenodo / LAMBDA. These are about 500 MB total, and will be downloaded into your astropy cache (`~/.astropy/`). The data are hosted by [datacentral.org.au/](https://datacentral.org.au), and are also available on [Zenodo](https://zenodo.org/record/3479985#.XaASx79S-AY). Examples --------- -To get a quick feel of what `PyGDSM` does, have a look at the +To get a quick feel of what `PyGDSM` does, have a look at the [GSM2008 quickstart guide](http://nbviewer.ipython.org/github/telegraphic/PyGDSM/blob/master/docs/pygdsm_quickstart.ipynb), and the new [GSM2016 quickstart guide](http://nbviewer.ipython.org/github/telegraphic/PyGDSM/blob/master/docs/pygdsm2016_quickstart.ipynb). @@ -54,27 +55,27 @@ Q & A **Q. What's the difference between this and the `gsm.f` from the main GSM2008 website?** The `gsm.f` is a very basic Fortran code, which reads and writes values to and from ASCII files, and uses a command line interface for input. If you want to run this code - on an ancient computer with nothing but Fortran installed, then `gsm.f` is the way to go. - In contrast, `PyGDSM` is a Python code that leverages a lot of other Packages so that you - can do more stuff more efficiently. For example: you can view a sky model in a healpy - image; you can write a sky model to a Healpix FITS file; and believe it or not, the - Python implementation is *much faster*. Have a look at the + on an ancient computer with nothing but Fortran installed, then `gsm.f` is the way to go. + In contrast, `PyGDSM` is a Python code that leverages a lot of other Packages so that you + can do more stuff more efficiently. For example: you can view a sky model in a healpy + image; you can write a sky model to a Healpix FITS file; and believe it or not, the + Python implementation is *much faster*. Have a look at the [quickstart guide](http://nbviewer.ipython.org/github/telegraphic/PyGDSM/blob/master/docs/pygdsm_quickstart.ipynb) to get a feel for what `PyGDSM` does. -**Q. Are the outputs of `gsm.f` and `pygdsm` identical?**. - **NO**. The cubic spline interpolation implementation differs, so values will differ by as +**Q. Are the outputs of `gsm.f` and `pygdsm` identical?**. + **NO**. The cubic spline interpolation implementation differs, so values will differ by as much as a few percent. The interpolation code used in `gsm.f` does not have an open-source - license (it's from [Numerical Recipes](http://www.nr.com/licenses/) ), so we haven't + license (it's from [Numerical Recipes](http://www.nr.com/licenses/) ), so we haven't implemented it (one could probably come up with an equivalent that didn't infringe). Nevertheless, the underlying PCA data are identical, and I've run tests to check that - the two outputs are indeed comparable. + the two outputs are indeed comparable. **Q. What's the difference between this and the [Zheng et. al. github repo](https://github.com/jeffzhen/gsm2016)?** `pygdsm` provides two classes: `GlobalSkyModel16()` and `GSMObserver16()`, which once instantiated - provide methods for programatically generating sky models. The Zheng et. al. github repo is a + provide methods for programatically generating sky models. The Zheng et. al. github repo is a simple, low-dependency, command line tool. As of PyGDSM 1.4.0, we have implemented improved interpolation - via cubic spline or PCHIP, which avoids discontinuities identified using the 2016 method. Have a look at the + via cubic spline or PCHIP, which avoids discontinuities identified using the 2016 method. Have a look at the [GSM2016 quickstart guide](http://nbviewer.ipython.org/github/telegraphic/PyGDSM/blob/master/docs/pygdsm2016_quickstart.ipynb) to get a feel for what `PyGDSM` does. @@ -84,7 +85,7 @@ Q & A files that come in the original `gsm.tar.gz` file. The next biggest thing is test data, so that the output can be compared against precomputed output from `gsm.f`. The package now also includes the Zheng et. al. data, which is another ~300 MB. - + References ---------- @@ -112,11 +113,11 @@ J. Dowell, G. B. Taylor, F. Schinzel, N. E. Kassim, K. Stovall MNRAS, 469, 4, 4537-4550 (2017) https://ui.adsabs.harvard.edu/abs/2017MNRAS.469.4537D/abstract -An improved source-subtracted and destriped 408-MHz all-sky map +An improved source-subtracted and destriped 408-MHz all-sky map M. Remazeilles, C. Dickinson,A.J. Banday, M. Bigot-Sazy, T. Ghosh MNRAS 451, 4, 4311-4327 (2014) https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.4311R/abstract - + ``` PyGSDM has an [ascl.net entry](https://ascl.net/1603.013): @@ -128,6 +129,5 @@ D. C. Price, 2016, 2.0.0, Astrophysics Source Code Library, 1603.013 License ------- -All *code* in PyGDSM is licensed under the MIT license (not the underlying *data*). +All *code* in PyGDSM is licensed under the MIT license (not the underlying *data*). The PCA data, by Zheng et. al. is licensed under MIT also (see https://github.com/jeffzhen/gsm2016). - diff --git a/pygdsm/__init__.py b/pygdsm/__init__.py index 7b53832..a2c37bc 100644 --- a/pygdsm/__init__.py +++ b/pygdsm/__init__.py @@ -1,22 +1,15 @@ - -from .gsm08 import GlobalSkyModel -from .gsm08 import GSMObserver - -from .gsm16 import GlobalSkyModel16 -from .gsm16 import GSMObserver16 - -from .lfsm import LowFrequencySkyModel -from .lfsm import LFSMObserver - -from .haslam import HaslamSkyModel -from .haslam import HaslamObserver +from .gsm08 import GlobalSkyModel, GSMObserver +from .gsm16 import GlobalSkyModel16, GSMObserver16 +from .haslam import HaslamObserver, HaslamSkyModel +from .lfsm import LFSMObserver, LowFrequencySkyModel +from .component_data import download_map_data as download_map_data # Add aliases GlobalSkyModel08 = GlobalSkyModel -GSMObserver08 = GSMObserver +GSMObserver08 = GSMObserver -def init_gsm(gsm_name: str='gsm08'): - """ Initialize a GDSM object by ID/name +def init_gsm(gsm_name: str = "gsm08"): + """Initialize a GDSM object by ID/name Returns a diffuse sky model (subclass of BaseSkyModel), based on one of: * **GSM08:** A model of diffuse Galactic radio emission from 10 MHz to 100 GHz, @@ -36,20 +29,20 @@ def init_gsm(gsm_name: str='gsm08'): """ gsm_name = gsm_name.lower().strip() match gsm_name: - case 'gsm': # Shorthand for GSM08 + case "gsm": # Shorthand for GSM08 return GlobalSkyModel() - case 'gsm08': + case "gsm08": return GlobalSkyModel() - case 'gsm16': + case "gsm16": return GlobalSkyModel16() - case 'lfsm': + case "lfsm": return LowFrequencySkyModel() - case 'haslam': + case "haslam": return HaslamSkyModel() -def init_observer(gsm_name: str='gsm08'): - """ Initialize a GDSM Observer object by ID/name +def init_observer(gsm_name: str = "gsm08"): + """Initialize a GDSM Observer object by ID/name Returns an observer (subclass of BaseObserver), where the diffuse sky is created from one of: * **GSM08:** A model of diffuse Galactic radio emission from 10 MHz to 100 GHz, @@ -69,13 +62,13 @@ def init_observer(gsm_name: str='gsm08'): """ gsm_name = gsm_name.lower().strip() match gsm_name: - case 'gsm': # Shorthand for GSM08 + case "gsm": # Shorthand for GSM08 return GSMObserver() - case 'gsm08': + case "gsm08": return GSMObserver() - case 'gsm16': + case "gsm16": return GSMObserver16() - case 'lfsm': + case "lfsm": return LFSMObserver() - case 'haslam': - return HaslamObserver() \ No newline at end of file + case "haslam": + return HaslamObserver() diff --git a/pygdsm/base_observer.py b/pygdsm/base_observer.py index 8013aa4..1386f82 100644 --- a/pygdsm/base_observer.py +++ b/pygdsm/base_observer.py @@ -1,13 +1,15 @@ import ephem import healpy as hp import numpy as np +from astropy.coordinates import SkyCoord 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. + """Observer of the Global Sky Model. Generates the Observed sky, for a given point on Earth. Applies the necessary rotations and coordinate transformations @@ -20,7 +22,7 @@ class BaseObserver(ephem.Observer): """ def __init__(self, gsm): - """ Initialize the Observer object. + """Initialize the Observer object. Calls ephem.Observer.__init__ function and adds on gsm """ @@ -34,7 +36,7 @@ def _setup(self): self._time = Time(self.date.datetime()) # Generate mapping from pix <-> angles self.gsm.generate(self._freq) - self._n_pix = hp.get_map_size(self.gsm.generated_map_data) + self._n_pix = hp.get_map_size(self.gsm.generated_map_data) self._n_side = hp.npix2nside(self._n_pix) self._theta, self._phi = hp.pix2ang(self._n_side, np.arange(self._n_pix)) @@ -44,7 +46,7 @@ def _setup(self): self._observed_dec = None def generate(self, freq=None, obstime=None): - """ Generate the observed sky for the observer, based on the GSM. + """Generate the observed sky for the observer, based on the GSM. Parameters ---------- @@ -72,37 +74,39 @@ 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.date = obstime.to_datetime() + self._time = Time( + obstime + ) # This will catch datetimes, but Time() object should be passed + self.date = obstime.to_datetime() # 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: # 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')) + 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) + 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 = 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']) + rot = hp.Rotator(coord=["G", "C"]) eq_theta, eq_phi = rot(self._theta, self._phi) # Convert from Equatorial colatitude and longitude to normal RA and DEC - dec = 90.0 - np.abs(eq_theta*(180/np.pi)) - ra = ( (eq_phi + 2*np.pi) % (2*np.pi) )*(180/np.pi) + dec = 90.0 - np.abs(eq_theta * (180 / np.pi)) + ra = ((eq_phi + 2 * np.pi) % (2 * np.pi)) * (180 / np.pi) # 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) + hrot = hp.Rotator(rot=[ra_zen, dec_zen], coord=["G", "C"], inv=True) g0, g1 = hrot(self._theta, self._phi) pix0 = hp.ang2pix(self._n_side, g0, g1) self._pix0 = pix0 @@ -121,7 +125,7 @@ def generate(self, freq=None, obstime=None): return self.observed_sky def view(self, logged=False, show=False, **kwargs): - """ View the local sky, in orthographic projection. + """View the local sky, in orthographic projection. Parameters ---------- @@ -140,12 +144,12 @@ def view(self, logged=False, show=False, **kwargs): @property def observed_gsm(self): - """ Return the GSM (Mollweide), with below-horizon area masked. """ + """Return the GSM (Mollweide), with below-horizon area masked.""" sky = self.observed_sky # Get RA and DEC of zenith ra_rad, dec_rad = self.radec_of(0, np.pi / 2) - ra_deg = ra_rad / np.pi * 180 + ra_deg = ra_rad / np.pi * 180 dec_deg = dec_rad / np.pi * 180 # Apply rotation @@ -154,14 +158,14 @@ def observed_gsm(self): pix0 = hp.ang2pix(self._n_side, g0, g1) sky = sky[pix0] - coordrotate = hp.Rotator(coord=['C', 'G'], inv=True) + coordrotate = hp.Rotator(coord=["C", "G"], inv=True) g0, g1 = coordrotate(self._theta, self._phi) pix0 = hp.ang2pix(self._n_side, g0, g1) sky = sky[pix0] return sky def view_observed_gsm(self, logged=False, show=False, **kwargs): - """ View the GSM (Mollweide), with below-horizon area masked. + """View the GSM (Mollweide), with below-horizon area masked. Args: logged (bool): Apply log2 to data (default False) @@ -175,7 +179,7 @@ def view_observed_gsm(self, logged=False, show=False, **kwargs): if logged: sky = np.log2(sky) - hp.mollview(sky, coord='G', **kwargs) + hp.mollview(sky, coord="G", **kwargs) if show: show_plt() return sky diff --git a/pygdsm/base_skymodel.py b/pygdsm/base_skymodel.py index 58e2523..08b55fa 100644 --- a/pygdsm/base_skymodel.py +++ b/pygdsm/base_skymodel.py @@ -1,10 +1,13 @@ import os -import numpy as np + import h5py import healpy as hp +import numpy as np from astropy.io import fits + from .plot_utils import show_plt + def is_fits(filepath): """ Check if file is a FITS file @@ -15,11 +18,13 @@ def is_fits(filepath): filepath: str Path to file """ - FITS_SIGNATURE = (b"\x53\x49\x4d\x50\x4c\x45\x20\x20\x3d\x20\x20\x20\x20\x20" - b"\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20" - b"\x20\x54") + FITS_SIGNATURE = ( + b"\x53\x49\x4d\x50\x4c\x45\x20\x20\x3d\x20\x20\x20\x20\x20" + b"\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20\x20" + b"\x20\x54" + ) try: - with open(str(filepath),'rb') as f: + with open(str(filepath), "rb") as f: return f.read(30) == FITS_SIGNATURE except FileNotFoundError as e: print(e) @@ -27,10 +32,10 @@ def is_fits(filepath): class BaseSkyModel(object): - """ Global sky model (GSM) class for generating sky models. - """ + """Global sky model (GSM) class for generating sky models.""" + def __init__(self, name, filepath, freq_unit, data_unit, basemap): - """ Initialise basic sky model class + """Initialise basic sky model class Parameters ---------- @@ -65,7 +70,7 @@ def generate(self, freqs): raise NotImplementedError def view(self, idx=0, logged=False, show=False): - """ View generated map using healpy's mollweide projection. + """View generated map using healpy's mollweide projection. Parameters ---------- @@ -78,7 +83,9 @@ def view(self, idx=0, logged=False, show=False): """ if self.generated_map_data is None: - raise RuntimeError("No GSM map has been generated yet. Run generate() first.") + raise RuntimeError( + "No GSM map has been generated yet. Run generate() first." + ) if self.generated_map_data.ndim == 2: gmap = self.generated_map_data[idx] @@ -90,13 +97,15 @@ def view(self, idx=0, logged=False, show=False): if logged: gmap = np.log2(gmap) - hp.mollview(gmap, coord='G', title='%s %s, %s' % (self.name, str(freq), self.basemap)) + hp.mollview( + gmap, coord="G", title="%s %s, %s" % (self.name, str(freq), self.basemap) + ) if show: show_plt() def get_sky_temperature(self, coords, freqs=None, include_cmb=True): - """ Get sky temperature at given coordinates. + """Get sky temperature at given coordinates. Returns sky temperature at a given SkyCoord (e.g. Ra/Dec or galactic l/b). Useful for estimating sky contribution to system temperature. @@ -117,19 +126,20 @@ def get_sky_temperature(self, coords, freqs=None, include_cmb=True): if freqs is not None: self.generate(freqs) - pix = hp.ang2pix(self.nside, coords.galactic.l.deg, coords.galactic.b.deg, lonlat=True) + pix = hp.ang2pix( + self.nside, coords.galactic.l.deg, coords.galactic.b.deg, lonlat=True + ) if self.generated_map_data.ndim == 2: return self.generated_map_data[:, pix] + T_cmb else: return self.generated_map_data[pix] + T_cmb - def write_fits(self, filename): - """ Write out map data as FITS file. + """Write out map data as FITS file. Parameters ---------- filename: str file name for output FITS file """ - hp.write_map(filename, self.generated_map_data, column_units=self.data_unit) \ No newline at end of file + hp.write_map(filename, self.generated_map_data, column_units=self.data_unit) diff --git a/pygdsm/component_data.py b/pygdsm/component_data.py index c56f4bf..7b8b50f 100644 --- a/pygdsm/component_data.py +++ b/pygdsm/component_data.py @@ -1,19 +1,25 @@ -from astropy.utils.data import conf, download_file +from astropy.utils.data import download_file +from astropy.config import get_cache_dir -# Zenodo can be slow to respond, so set timeout to 30 seconds -conf.remote_timeout = 30 +FILE_HOST = "datacentral.org.au" +GSM_DATA_URL = "https://apps.datacentral.org.au/pygdsm/data/gsm_components.h5" +GSM2016_DATA_URL = "https://apps.datacentral.org.au/pygdsm/data/gsm2016_components.h5" +GSM2008_TEST_DATA_URL = "https://apps.datacentral.org.au/pygdsm/data/gsm_fortran_test_data.h5" +LFSM_DATA_URL = "https://apps.datacentral.org.au/pygdsm/data/lfsm.h5" +HASLAM_DATA_URL = "https://lambda.gsfc.nasa.gov/data/foregrounds/haslam_2014/haslam408_dsds_Remazeilles2014.fits" -GSM_FILEPATH = download_file('https://zenodo.org/record/3835582/files/gsm_components.h5?download=1', - cache=True, show_progress=True) +DATA_URLS = [ + GSM_DATA_URL, + GSM2016_DATA_URL, + GSM2008_TEST_DATA_URL, + LFSM_DATA_URL, + HASLAM_DATA_URL + ] -GSM2016_FILEPATH = download_file('https://zenodo.org/record/3835582/files/gsm2016_components.h5?download=1', - cache=True, show_progress=True) - -GSM2008_TEST_DATA = download_file('https://zenodo.org/record/3835582/files/gsm_fortran_test_data.h5?download=1', - cache=True, show_progress=True) - -LFSM_FILEPATH = download_file('https://zenodo.org/record/3835582/files/lfsm.h5?download=1', - cache=True, show_progress=True) - -HASLAM_FILEPATH = download_file('https://lambda.gsfc.nasa.gov/data/foregrounds/haslam_2014/haslam408_dsds_Remazeilles2014.fits', - cache=True, show_progress=True) +def download_map_data(): + """Download component data.""" + print("Checking for component data files...") + for URL in DATA_URLS: + download_file(URL, cache=True, show_progress=True) + CACHE_PATH = get_cache_dir() + print(f"Data saved in {CACHE_PATH}") \ No newline at end of file diff --git a/pygdsm/gsm08.py b/pygdsm/gsm08.py index 83c0cab..4ed3ce6 100644 --- a/pygdsm/gsm08.py +++ b/pygdsm/gsm08.py @@ -19,25 +19,28 @@ """ import numpy as np -from scipy.interpolate import interp1d, pchip -import h5py from astropy import units -import healpy as hp - +from astropy.utils.data import download_file +from scipy.interpolate import interp1d, pchip -from .component_data import GSM_FILEPATH -from .plot_utils import show_plt from .base_observer import BaseObserver from .base_skymodel import BaseSkyModel +from .component_data import GSM_DATA_URL +from .plot_utils import show_plt T_CMB = 2.725 class GlobalSkyModel(BaseSkyModel): - """ Global sky model (GSM) class for generating sky models. - """ + """Global sky model (GSM) class for generating sky models.""" - def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip', include_cmb=False): - """ Global sky model (GSM) class for generating sky models. + def __init__( + self, + freq_unit="MHz", + basemap="haslam", + interpolation="pchip", + include_cmb=False, + ): + """Global sky model (GSM) class for generating sky models. Upon initialization, the map PCA data are loaded into memory and interpolation functions are pre-computed. @@ -72,23 +75,33 @@ def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip', inc """ + # download component data as needed using astropy cache + GSM_FILEPATH = download_file( + GSM_DATA_URL, + cache=True, + show_progress=True, + ) + try: - assert basemap in {'5deg', 'wmap', 'haslam'} + assert basemap in {"5deg", "wmap", "haslam"} except AssertionError: - raise RuntimeError("GSM basemap unknown: %s. Choose '5deg', 'haslam' or 'wmap'" % basemap) + raise RuntimeError( + "GSM basemap unknown: %s. Choose '5deg', 'haslam' or 'wmap'" % basemap + ) try: - assert interpolation in {'cubic', 'pchip'} + assert interpolation in {"cubic", "pchip"} except AssertionError: raise RuntimeError("Interpolation must be set to either 'cubic' or 'pchip'") - data_unit = 'K' - super(GlobalSkyModel, self).__init__('GSM2008', GSM_FILEPATH, freq_unit, data_unit, basemap) + data_unit = "K" + super(GlobalSkyModel, self).__init__( + "GSM2008", GSM_FILEPATH, freq_unit, data_unit, basemap + ) self.interpolation_method = interpolation self.include_cmb = include_cmb - self.pca_map_data = None self.interp_comps = None self.update_interpolants() @@ -98,37 +111,39 @@ def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip', inc self.nside = 512 def update_interpolants(self): - # Choose the PCA map to load from the HDF5 file - pca_map_dict = {"5deg": "component_maps_5deg", - "haslam": "component_maps_408locked", - "wmap": "component_maps_23klocked"} + # Choose the PCA map to load from the HDF5 file + pca_map_dict = { + "5deg": "component_maps_5deg", + "haslam": "component_maps_408locked", + "wmap": "component_maps_23klocked", + } pca_map_key = pca_map_dict[self.basemap] self.pca_map_data = self.h5[pca_map_key][:] # Now, load the PCA eigenvalues pca_table = self.h5["components"][:] pca_freqs_mhz = pca_table[:, 0] - pca_scaling = pca_table[:, 1] - pca_comps = pca_table[:, 2:].T + pca_scaling = pca_table[:, 1] + pca_comps = pca_table[:, 2:].T # Interpolate to the desired frequency values ln_pca_freqs = np.log(pca_freqs_mhz) - if self.interpolation_method == 'cubic': - spl_scaling = interp1d(ln_pca_freqs, np.log(pca_scaling), kind='cubic') - spl1 = interp1d(ln_pca_freqs, pca_comps[0], kind='cubic') - spl2 = interp1d(ln_pca_freqs, pca_comps[1], kind='cubic') - spl3 = interp1d(ln_pca_freqs, pca_comps[2], kind='cubic') + if self.interpolation_method == "cubic": + spl_scaling = interp1d(ln_pca_freqs, np.log(pca_scaling), kind="cubic") + spl1 = interp1d(ln_pca_freqs, pca_comps[0], kind="cubic") + spl2 = interp1d(ln_pca_freqs, pca_comps[1], kind="cubic") + spl3 = interp1d(ln_pca_freqs, pca_comps[2], kind="cubic") else: spl_scaling = pchip(ln_pca_freqs, np.log(pca_scaling)) - spl1 = pchip(ln_pca_freqs, pca_comps[0]) - spl2 = pchip(ln_pca_freqs, pca_comps[1]) - spl3 = pchip(ln_pca_freqs, pca_comps[2]) + spl1 = pchip(ln_pca_freqs, pca_comps[0]) + spl2 = pchip(ln_pca_freqs, pca_comps[1]) + spl3 = pchip(ln_pca_freqs, pca_comps[2]) self.interp_comps = (spl_scaling, spl1, spl2, spl3) def generate(self, freqs): - """ Generate a global sky model at a given frequency or frequencies + """Generate a global sky model at a given frequency or frequencies Parameters ---------- @@ -144,7 +159,7 @@ def generate(self, freqs): """ # convert frequency values into Hz freqs = np.array(freqs) * units.Unit(self.freq_unit) - freqs_mhz = freqs.to('MHz').value + freqs_mhz = freqs.to("MHz").value if isinstance(freqs_mhz, float): freqs_mhz = np.array([freqs_mhz]) @@ -156,16 +171,16 @@ def generate(self, freqs): raise RuntimeError("Frequency values lie outside 10 MHz < f < 94 GHz") # Load interpolators and do interpolation - ln_freqs = np.log(freqs_mhz) + ln_freqs = np.log(freqs_mhz) spl_scaling, spl1, spl2, spl3 = self.interp_comps comps = np.row_stack((spl1(ln_freqs), spl2(ln_freqs), spl3(ln_freqs))) scaling = np.exp(spl_scaling(ln_freqs)) # Finally, compute the dot product via einsum (awesome function) # c=comp, f=freq, p=pixel. We want to dot product over c for each freq - #print comps.shape, self.pca_map_data.shape, scaling.shape - map_out = np.einsum('cf,pc,f->fp', comps, self.pca_map_data, scaling) - + # print comps.shape, self.pca_map_data.shape, scaling.shape + map_out = np.einsum("cf,pc,f->fp", comps, self.pca_map_data, scaling) + if self.include_cmb: map_out += T_CMB @@ -174,8 +189,6 @@ def generate(self, freqs): self.generated_map_data = map_out self.generated_map_freqs = freqs - - return map_out def set_basemap(self, new_basemap): @@ -199,8 +212,8 @@ def set_interpolation_method(self, new_method): class GSMObserver(BaseObserver): def __init__(self): - """ Initialize the Observer object. + """Initialize the Observer object. Calls ephem.Observer.__init__ function and adds on gsm """ - super(GSMObserver, self).__init__(gsm=GlobalSkyModel) \ No newline at end of file + super(GSMObserver, self).__init__(gsm=GlobalSkyModel) diff --git a/pygdsm/gsm16.py b/pygdsm/gsm16.py index 9d9f13d..323c5c1 100644 --- a/pygdsm/gsm16.py +++ b/pygdsm/gsm16.py @@ -1,14 +1,12 @@ +import healpy as hp import numpy as np -from scipy.interpolate import interp1d, pchip -import h5py from astropy import units -import healpy as hp -import ephem +from astropy.utils.data import download_file +from scipy.interpolate import interp1d, pchip -from .component_data import GSM2016_FILEPATH -from .plot_utils import show_plt from .base_observer import BaseObserver from .base_skymodel import BaseSkyModel +from .component_data import GSM2016_DATA_URL kB = 1.38065e-23 C = 2.99792e8 @@ -17,21 +15,22 @@ hoverk = h / kB -def K_CMB2MJysr(K_CMB, nu):#in Kelvin and Hz - B_nu = 2 * (h * nu)* (nu / C)**2 / (np.exp(hoverk * nu / T) - 1) - conversion_factor = (B_nu * C / nu / T)**2 / 2 * np.exp(hoverk * nu / T) / kB - return K_CMB * conversion_factor * 1e20#1e-26 for Jy and 1e6 for MJy +def K_CMB2MJysr(K_CMB, nu): # in Kelvin and Hz + B_nu = 2 * (h * nu) * (nu / C) ** 2 / (np.exp(hoverk * nu / T) - 1) + conversion_factor = (B_nu * C / nu / T) ** 2 / 2 * np.exp(hoverk * nu / T) / kB + return K_CMB * conversion_factor * 1e20 # 1e-26 for Jy and 1e6 for MJy -def K_RJ2MJysr(K_RJ, nu):#in Kelvin and Hz - conversion_factor = 2 * (nu / C)**2 * kB - return K_RJ * conversion_factor * 1e20#1e-26 for Jy and 1e6 for MJy + +def K_RJ2MJysr(K_RJ, nu): # in Kelvin and Hz + conversion_factor = 2 * (nu / C) ** 2 * kB + return K_RJ * conversion_factor * 1e20 # 1e-26 for Jy and 1e6 for MJy def rotate_map(hmap, rot_theta, rot_phi, nest=True): nside = hp.npix2nside(len(hmap)) # Get theta, phi for non-rotated map - t, p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside)), nest= nest) # theta, phi + t, p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside)), nest=nest) # theta, phi # Define a rotator r = hp.Rotator(deg=False, rot=[rot_phi, rot_theta]) @@ -40,17 +39,25 @@ def rotate_map(hmap, rot_theta, rot_phi, nest=True): trot, prot = r(t, p) # Inerpolate map onto these co-ordinates - rot_map = hp.get_interp_val(hmap, trot, prot, nest= nest) + rot_map = hp.get_interp_val(hmap, trot, prot, nest=nest) return rot_map class GlobalSkyModel16(BaseSkyModel): - """ Global sky model (GSM) class for generating sky models. - """ - - def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot=0, phi_rot=0, interpolation='pchip', include_cmb=False): - """ Global sky model (GSM) class for generating sky models. + """Global sky model (GSM) class for generating sky models.""" + + def __init__( + self, + freq_unit="MHz", + data_unit="TCMB", + resolution="hi", + theta_rot=0, + phi_rot=0, + interpolation="pchip", + include_cmb=False, + ): + """Global sky model (GSM) class for generating sky models. Upon initialization, the map PCA data are loaded into memory and interpolation functions are pre-computed. @@ -83,41 +90,55 @@ def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot """ - if data_unit not in ['MJysr', 'TCMB', 'TRJ']: - raise RuntimeError("UNIT ERROR: %s not supported. Only MJysr, TCMB, TRJ are allowed." % data_unit) - - if resolution.lower() in ('hi', 'high', 'h'): - resolution = 'hi' - elif resolution.lower() in ('low', 'lo', 'l'): - resolution = 'low' + # download component data as needed using astropy cache + GSM2016_FILEPATH = download_file( + GSM2016_DATA_URL, + cache=True, + show_progress=True, + ) + + if data_unit not in ["MJysr", "TCMB", "TRJ"]: + raise RuntimeError( + "UNIT ERROR: %s not supported. Only MJysr, TCMB, TRJ are allowed." + % data_unit + ) + + if resolution.lower() in ("hi", "high", "h"): + resolution = "hi" + elif resolution.lower() in ("low", "lo", "l"): + resolution = "low" else: - raise RuntimeError("RESOLUTION ERROR: Must be either hi or low, not %s" % resolution) + raise RuntimeError( + "RESOLUTION ERROR: Must be either hi or low, not %s" % resolution + ) - super(GlobalSkyModel16, self).__init__('GSM2016', GSM2016_FILEPATH, freq_unit, data_unit, basemap='') + super(GlobalSkyModel16, self).__init__( + "GSM2016", GSM2016_FILEPATH, freq_unit, data_unit, basemap="" + ) self.interpolation_method = interpolation self.resolution = resolution self.include_cmb = include_cmb # Map data to load - labels = ['Synchrotron', 'CMB', 'HI', 'Dust1', 'Dust2', 'Free-Free'] + labels = ["Synchrotron", "CMB", "HI", "Dust1", "Dust2", "Free-Free"] self.n_comp = len(labels) - if resolution=='hi': + if resolution == "hi": self.nside = 1024 - self.map_ni = np.array([self.h5['highres_%s_map'%lb][:] for lb in labels]) + self.map_ni = np.array([self.h5["highres_%s_map" % lb][:] for lb in labels]) else: self.nside = 64 - self.map_ni = np.array(self.h5['lowres_maps']) + self.map_ni = np.array(self.h5["lowres_maps"]) - self.spec_nf = self.h5['spectra'][:] + self.spec_nf = self.h5["spectra"][:] if theta_rot or phi_rot: - for i,map in enumerate(self.map_ni): + for i, map in enumerate(self.map_ni): self.map_ni[i] = rotate_map(map, theta_rot, phi_rot, nest=True) def generate(self, freqs): - """ Generate a global sky model at a given frequency or frequencies + """Generate a global sky model at a given frequency or frequencies Parameters ---------- @@ -134,7 +155,7 @@ def generate(self, freqs): # convert frequency values into Hz freqs = np.array(freqs) * units.Unit(self.freq_unit) - freqs_ghz = freqs.to('GHz').value + freqs_ghz = freqs.to("GHz").value if isinstance(freqs_ghz, float): freqs_ghz = np.array([freqs_ghz]) @@ -153,65 +174,69 @@ def generate(self, freqs): spec_nf = self.spec_nf nfreq = spec_nf.shape[1] - + # Now borrow code from the orignal GSM2008 model to do a sensible interpolation pca_freqs_ghz = spec_nf[0] - pca_scaling = spec_nf[1] - pca_comps = spec_nf[2:] - # Interpolate to the desired frequency values + pca_scaling = spec_nf[1] + pca_comps = spec_nf[2:] + # Interpolate to the desired frequency values ln_pca_freqs = np.log(pca_freqs_ghz) - if self.interpolation_method == 'cubic': - spl_scaling = interp1d(ln_pca_freqs, np.log(pca_scaling), kind='cubic') - spl1 = interp1d(ln_pca_freqs, pca_comps[0], kind='cubic') - spl2 = interp1d(ln_pca_freqs, pca_comps[1], kind='cubic') - spl3 = interp1d(ln_pca_freqs, pca_comps[2], kind='cubic') - spl4 = interp1d(ln_pca_freqs, pca_comps[3], kind='cubic') - spl5 = interp1d(ln_pca_freqs, pca_comps[4], kind='cubic') - spl6 = interp1d(ln_pca_freqs, pca_comps[5], kind='cubic') + if self.interpolation_method == "cubic": + spl_scaling = interp1d(ln_pca_freqs, np.log(pca_scaling), kind="cubic") + spl1 = interp1d(ln_pca_freqs, pca_comps[0], kind="cubic") + spl2 = interp1d(ln_pca_freqs, pca_comps[1], kind="cubic") + spl3 = interp1d(ln_pca_freqs, pca_comps[2], kind="cubic") + spl4 = interp1d(ln_pca_freqs, pca_comps[3], kind="cubic") + spl5 = interp1d(ln_pca_freqs, pca_comps[4], kind="cubic") + spl6 = interp1d(ln_pca_freqs, pca_comps[5], kind="cubic") else: spl_scaling = pchip(ln_pca_freqs, np.log(pca_scaling)) - spl1 = pchip(ln_pca_freqs, pca_comps[0]) - spl2 = pchip(ln_pca_freqs, pca_comps[1]) - spl3 = pchip(ln_pca_freqs, pca_comps[2]) - spl4 = pchip(ln_pca_freqs, pca_comps[3]) - spl5 = pchip(ln_pca_freqs, pca_comps[4]) - spl6 = pchip(ln_pca_freqs, pca_comps[5]) - + spl1 = pchip(ln_pca_freqs, pca_comps[0]) + spl2 = pchip(ln_pca_freqs, pca_comps[1]) + spl3 = pchip(ln_pca_freqs, pca_comps[2]) + spl4 = pchip(ln_pca_freqs, pca_comps[3]) + spl5 = pchip(ln_pca_freqs, pca_comps[4]) + spl6 = pchip(ln_pca_freqs, pca_comps[5]) + self.interp_comps = (spl_scaling, spl1, spl2, spl3, spl4, spl5, spl6) - + ln_freqs = np.log(freqs_ghz) - comps = np.row_stack((spl1(ln_freqs), spl2(ln_freqs), spl3(ln_freqs), spl4(ln_freqs), spl5(ln_freqs), spl6(ln_freqs))) + comps = np.row_stack(( + spl1(ln_freqs), + spl2(ln_freqs), + spl3(ln_freqs), + spl4(ln_freqs), + spl5(ln_freqs), + spl6(ln_freqs), + )) scaling = np.exp(spl_scaling(ln_freqs)) - + # Finally, compute the dot product via einsum (awesome function) # c=comp, f=freq, p=pixel. We want to dot product over c for each freq - #print comps.shape, self.pca_map_data.shape, scaling.shape - - output = np.single(np.einsum('cf,pc,f->fp', comps, map_ni.T, scaling)) - + # print comps.shape, self.pca_map_data.shape, scaling.shape + + output = np.single(np.einsum("cf,pc,f->fp", comps, map_ni.T, scaling)) for ifreq, freq in enumerate(freqs_ghz): - output[ifreq] = hp.pixelfunc.reorder(output[ifreq], n2r=True) # DCP 2024.03.29 - Add CMB if requested if self.include_cmb: - output[ifreq] += K_CMB2MJysr(T, 1e9 * freq) - - if self.data_unit == 'TCMB': - conversion = 1. / K_CMB2MJysr(1., 1e9 * freq) - elif self.data_unit == 'TRJ': - conversion = 1. / K_RJ2MJysr(1., 1e9 * freq) + output[ifreq] += K_CMB2MJysr(T, 1e9 * freq) + + if self.data_unit == "TCMB": + conversion = 1.0 / K_CMB2MJysr(1.0, 1e9 * freq) + elif self.data_unit == "TRJ": + conversion = 1.0 / K_RJ2MJysr(1.0, 1e9 * freq) else: - conversion = 1. + conversion = 1.0 output[ifreq] *= conversion - if len(output) == 1: output = output[0] - #else: + # else: # map_data = np.row_stack(output) self.generated_map_freqs = freqs @@ -222,7 +247,7 @@ def generate(self, freqs): class GSMObserver16(BaseObserver): def __init__(self): - """ Initialize the Observer object. + """Initialize the Observer object. Calls ephem.Observer.__init__ function and adds on gsm """ diff --git a/pygdsm/haslam.py b/pygdsm/haslam.py index fd41613..a521ed2 100644 --- a/pygdsm/haslam.py +++ b/pygdsm/haslam.py @@ -1,19 +1,20 @@ +import healpy as hp import numpy as np from astropy import units -import healpy as hp -from scipy.interpolate import interp1d +from astropy.utils.data import download_file -from .component_data import HASLAM_FILEPATH -from .base_skymodel import BaseSkyModel from .base_observer import BaseObserver +from .base_skymodel import BaseSkyModel +from .component_data import HASLAM_DATA_URL T_CMB = 2.725 + class HaslamSkyModel(BaseSkyModel): - """ Haslam destriped, desourced sky model """ + """Haslam destriped, desourced sky model""" - def __init__(self, freq_unit='MHz', spectral_index=-2.6, include_cmb=False): - """ Generate a sky model at a given frequency based off Haslam 408 MHz map + def __init__(self, freq_unit="MHz", spectral_index=-2.6, include_cmb=False): + """Generate a sky model at a given frequency based off Haslam 408 MHz map Parameters ---------- @@ -24,7 +25,7 @@ def __init__(self, freq_unit='MHz', spectral_index=-2.6, include_cmb=False): Notes ----- The basemap is the destriped and desoruced Haslam map from Remazeilles (2015) - A T_CMB value of 2.725 K is subtracted to account for the microwave background component. + A T_CMB value of 2.725 K is subtracted to account for the microwave background component. This is a crude model; the sky's spectral index changes with sidereal time and a singular spectral index is not realistic. Here are some measured values @@ -42,15 +43,24 @@ def __init__(self, freq_unit='MHz', spectral_index=-2.6, include_cmb=False): References ---------- - Remazeilles et al (2015), improved Haslam map, https://doi.org/10.1093/mnras/stv1274 + Remazeilles et al (2015), improved Haslam map, https://doi.org/10.1093/mnras/stv1274 Mozdzen et al (2016), EDGES high-band, https://doi.org/10.1093/mnras/stw2696 Mozdzen et al (2019), EDGES low-band, https://doi.org/10.1093/mnras/sty3410 Dickinson et al (2019), C-BASS experiment, https://doi.org/10.1093/mnras/stz522 """ - data_unit = 'K' - basemap = 'Haslam' - - super(HaslamSkyModel, self).__init__('Haslam', HASLAM_FILEPATH, freq_unit, data_unit, basemap) + data_unit = "K" + basemap = "Haslam" + + # download component data as needed using astropy cache + HASLAM_FILEPATH = download_file( + HASLAM_DATA_URL, + cache=True, + show_progress=True, + ) + + super(HaslamSkyModel, self).__init__( + "Haslam", HASLAM_FILEPATH, freq_unit, data_unit, basemap + ) self.spectral_index = spectral_index self.data = hp.read_map(self.fits, dtype=np.float64) - T_CMB self.nside = 512 @@ -58,7 +68,7 @@ def __init__(self, freq_unit='MHz', spectral_index=-2.6, include_cmb=False): self.include_cmb = include_cmb def generate(self, freqs): - """ Generate a global sky model at a given frequency or frequencies + """Generate a global sky model at a given frequency or frequencies Parameters ---------- @@ -74,12 +84,14 @@ def generate(self, freqs): """ # convert frequency values into Hz freqs = np.array(freqs) * units.Unit(self.freq_unit) - freqs_mhz = freqs.to('MHz').value + freqs_mhz = freqs.to("MHz").value if isinstance(freqs_mhz, float): freqs_mhz = np.array([freqs_mhz]) - map_out = np.outer((freqs_mhz / 408.0) ** (self.spectral_index), self.data).squeeze() + map_out = np.outer( + (freqs_mhz / 408.0) ** (self.spectral_index), self.data + ).squeeze() if self.include_cmb: map_out += T_CMB @@ -90,8 +102,8 @@ def generate(self, freqs): class HaslamObserver(BaseObserver): def __init__(self): - """ Initialize the Observer object. + """Initialize the Observer object. Calls ephem.Observer.__init__ function and adds on gsm """ - super(HaslamObserver, self).__init__(gsm=HaslamSkyModel) \ No newline at end of file + super(HaslamObserver, self).__init__(gsm=HaslamSkyModel) diff --git a/pygdsm/lfsm.py b/pygdsm/lfsm.py index 4bf1818..fa9587d 100644 --- a/pygdsm/lfsm.py +++ b/pygdsm/lfsm.py @@ -1,19 +1,21 @@ +import healpy as hp import numpy as np from astropy import units -import healpy as hp +from astropy.utils.data import download_file from scipy.interpolate import interp1d -from .component_data import LFSM_FILEPATH -from .base_skymodel import BaseSkyModel from .base_observer import BaseObserver +from .base_skymodel import BaseSkyModel +from .component_data import LFSM_DATA_URL T_CMB = 2.725 + def rotate_equatorial_to_galactic(map): """ Given a map in equatorial coordinates, convert it to Galactic coordinates. """ - rotCG = hp.rotator.Rotator(coord=('C', 'G')) + rotCG = hp.rotator.Rotator(coord=("C", "G")) nSides = hp.pixelfunc.npix2nside(map.size) theta, phi = hp.pixelfunc.pix2ang(nSides, range(map.size)) theta_new, phi_new = rotCG(theta, phi, inv=True) @@ -22,11 +24,10 @@ def rotate_equatorial_to_galactic(map): class LowFrequencySkyModel(BaseSkyModel): - """ LWA1 Low Frequency Sky Model - """ + """LWA1 Low Frequency Sky Model""" - def __init__(self, freq_unit='MHz', include_cmb=False): - """ Global sky model (GSM) class for generating sky models. + def __init__(self, freq_unit="MHz", include_cmb=False): + """Global sky model (GSM) class for generating sky models. Parameters ---------- @@ -34,28 +35,38 @@ def __init__(self, freq_unit='MHz', include_cmb=False): include_cmb (bool): Choose whether to include the CMB. Defaults to False. A value of T_CMB = 2.725 K is used if True. """ - data_unit = 'K' - basemap = 'LFSS' - super(LowFrequencySkyModel, self).__init__('LFSM', LFSM_FILEPATH, freq_unit, data_unit, basemap) - - self.pca_map = self.h5['lfsm_component_maps_3.0deg.dat'][:] - self.pca_components = self.h5['lfsm_components.dat'][:] + data_unit = "K" + basemap = "LFSS" + + # download component data as needed using astropy cache + LFSM_FILEPATH = download_file( + LFSM_DATA_URL, + cache=True, + show_progress=True, + ) + + super(LowFrequencySkyModel, self).__init__( + "LFSM", LFSM_FILEPATH, freq_unit, data_unit, basemap + ) + + self.pca_map = self.h5["lfsm_component_maps_3.0deg.dat"][:] + self.pca_components = self.h5["lfsm_components.dat"][:] self.nside = 256 self.include_cmb = include_cmb - freqs = self.pca_components[:, 0] + freqs = self.pca_components[:, 0] sigmas = self.pca_components[:, 1] - comps = self.pca_components[:, 2:] + comps = self.pca_components[:, 2:] - self.scaleFunc = interp1d(np.log(freqs), np.log(sigmas), kind='slinear') + self.scaleFunc = interp1d(np.log(freqs), np.log(sigmas), kind="slinear") self.compFuncs = [] for i in range(comps.shape[1]): - self.compFuncs.append(interp1d(np.log(freqs), comps[:, i], kind='cubic')) + self.compFuncs.append(interp1d(np.log(freqs), comps[:, i], kind="cubic")) def generate(self, freqs): - """ Generate a global sky model at a given frequency or frequencies + """Generate a global sky model at a given frequency or frequencies Parameters ---------- @@ -71,7 +82,7 @@ def generate(self, freqs): """ # convert frequency values into Hz freqs = np.array(freqs) * units.Unit(self.freq_unit) - freqs_mhz = freqs.to('MHz').value + freqs_mhz = freqs.to("MHz").value if isinstance(freqs_mhz, float): freqs_mhz = np.array([freqs_mhz]) @@ -87,7 +98,7 @@ def generate(self, freqs): if freqs.ndim > 0: map_out = np.zeros(shape=(freqs.shape[0], hp.nside2npix(self.nside))) else: - map_out = np.zeros(shape=(1, hp.nside2npix(self.nside))) + map_out = np.zeros(shape=(1, hp.nside2npix(self.nside))) else: map_out = np.zeros(shape=(1, hp.nside2npix(self.nside))) @@ -96,7 +107,7 @@ def generate(self, freqs): map_out[ff] += compFunc(np.log(freqs_mhz[ff])) * self.pca_map[:, i] map_out[ff] *= np.exp(self.scaleFunc(np.log(freqs_mhz[ff]))) - map_out[ff] = rotate_equatorial_to_galactic(map_out[ff]) + map_out[ff] = rotate_equatorial_to_galactic(map_out[ff]) map_out = map_out.squeeze() @@ -110,8 +121,8 @@ def generate(self, freqs): class LFSMObserver(BaseObserver): def __init__(self): - """ Initialize the Observer object. + """Initialize the Observer object. Calls ephem.Observer.__init__ function and adds on gsm """ - super(LFSMObserver, self).__init__(gsm=LowFrequencySkyModel) \ No newline at end of file + super(LFSMObserver, self).__init__(gsm=LowFrequencySkyModel) diff --git a/pygdsm/plot_utils.py b/pygdsm/plot_utils.py index 51d2fb8..bc139d5 100644 --- a/pygdsm/plot_utils.py +++ b/pygdsm/plot_utils.py @@ -1,7 +1,8 @@ -def show_plt(): # pragma: no cover - """ Call plt.show() """ +def show_plt(): # pragma: no cover + """Call plt.show()""" try: - plt.show() # noqa: F821 + plt.show() # noqa: F821 except: import pylab as plt - plt.show() \ No newline at end of file + + plt.show() diff --git a/pygdsm/utils.py b/pygdsm/utils.py index 1e6fb85..f1f4a6f 100644 --- a/pygdsm/utils.py +++ b/pygdsm/utils.py @@ -1,10 +1,10 @@ -import numpy as np import healpy as hp +import numpy as np from astropy.coordinates import SkyCoord def hpix2sky(nside: int, pix_ids: np.ndarray) -> SkyCoord: - """ Convert a healpix pixel_id into a SkyCoord + """Convert a healpix pixel_id into a SkyCoord Args: nside (int): Healpix NSIDE parameter @@ -14,12 +14,12 @@ def hpix2sky(nside: int, pix_ids: np.ndarray) -> SkyCoord: sc (SkyCoord): Corresponding SkyCoordinates """ gl, gb = hp.pix2ang(nside, pix_ids, lonlat=True) - sc = SkyCoord(gl, gb, frame='galactic', unit=('deg', 'deg')) + sc = SkyCoord(gl, gb, frame="galactic", unit=("deg", "deg")) return sc def sky2hpix(nside: int, sc: SkyCoord) -> np.ndarray: - """ Convert a SkyCoord into a healpix pixel_id + """Convert a SkyCoord into a healpix pixel_id Args: nside (int): Healpix NSIDE parameter @@ -28,6 +28,6 @@ def sky2hpix(nside: int, sc: SkyCoord) -> np.ndarray: Returns: pix (np.array): Array of healpix pixel IDs """ - gl, gb = sc.galactic.l.to('deg').value, sc.galactic.b.to('deg').value + gl, gb = sc.galactic.l.to("deg").value, sc.galactic.b.to("deg").value pix = hp.ang2pix(nside, gl, gb, lonlat=True) return pix diff --git a/setup.py b/setup.py index 9a89064..5980f2e 100644 --- a/setup.py +++ b/setup.py @@ -1,80 +1,56 @@ -"""A setuptools based setup module. -""" +"""A setuptools based setup module.""" # Always prefer setuptools over distutils -from setuptools import setup, find_packages # To use a consistent encoding from codecs import open from os import path +from setuptools import find_packages, setup + here = path.abspath(path.dirname(__file__)) # Get the long description from the relevant file -with open(path.join(here, 'README.md'), encoding='utf-8') as f: +with open(path.join(here, "README.md"), encoding="utf-8") as f: long_description = f.read() -with open("requirements.txt", 'r') as fh: +with open("requirements.txt", "r") as fh: requirements = fh.read().splitlines() -with open("requirements_test.txt", 'r') as fh: +with open("requirements_test.txt", "r") as fh: test_requirements = fh.read().splitlines() setup( - name='pygdsm', - - # Versions should comply with PEP440. For a discussion on single-sourcing - # the version across setup.py and the project code, see - # https://packaging.python.org/en/latest/single_source_version.html - version='1.5.4', - - description='Python Global Sky Model of diffuse Galactic radio emission', + name="pygdsm", + version="1.5.5", + description="Python Global Sky Model of diffuse Galactic radio emission", long_description=long_description, - long_description_content_type='text/markdown', - + long_description_content_type="text/markdown", # The project's main homepage. - url='https://github.com/telegraphic/pygdsm', - + url="https://github.com/telegraphic/pygdsm", # Author details - author='Danny C. Price', - author_email='daniel.price@skao.int', - + author="Danny C. Price", + author_email="daniel.price@skao.int", # Choose your license - license='MIT', - + license="MIT", # See https://pypi.python.org/pypi?%3Aaction=list_classifiers classifiers=[ - # How mature is this project? Common values are - # 3 - Alpha - # 4 - Beta - # 5 - Production/Stable - 'Development Status :: 5 - Production/Stable', - - # Indicate who your project is intended for - 'Intended Audience :: Science/Research', - 'Topic :: Scientific/Engineering :: Astronomy', - - # Pick your license as you wish (should match "license" above) - 'License :: OSI Approved :: MIT License', - - # Specify the Python versions you support here. In particular, ensure - # that you indicate whether you support Python 2, Python 3 or both. - 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.10', - 'Programming Language :: Python :: 3.11', + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Astronomy", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", ], - # What does your project relate to? - keywords='radio astronomy sky model galactic diffuse emission', - + keywords="radio astronomy sky model galactic diffuse emission", # You can just specify the packages manually here if your project is # simple. Or you can use find_packages(). - packages=find_packages(exclude=['contrib', 'docs', 'tests*']), - + packages=find_packages(exclude=["contrib", "docs", "tests*"]), # List run-time dependencies here. These will be installed by pip when # your project is installed. For an analysis of "install_requires" vs pip's # requirements files see: # https://packaging.python.org/en/latest/requirements.html install_requires=requirements, tests_require=test_requirements, - ) diff --git a/tests/test_base_skymodel.py b/tests/test_base_skymodel.py index 524e1ea..0915a76 100644 --- a/tests/test_base_skymodel.py +++ b/tests/test_base_skymodel.py @@ -5,25 +5,45 @@ Tests for GSM base skymodel """ +from pathlib import Path + import pytest -from pygdsm.base_skymodel import BaseSkyModel -from pygdsm.component_data import GSM_FILEPATH -from pathlib import Path +from astropy.utils.data import download_file +from pygdsm.base_skymodel import BaseSkyModel +from pygdsm.component_data import GSM_DATA_URL -def test_base_skymodel_init(): +# download component data as needed using astropy cache +GSM_FILEPATH = download_file( + GSM_DATA_URL, + cache=True, + show_progress=True, +) - gsm = BaseSkyModel('TEST_GSM', GSM_FILEPATH, - freq_unit='MHz', basemap='haslam', data_unit='K') +def test_base_skymodel_init(): + gsm = BaseSkyModel( + "TEST_GSM", GSM_FILEPATH, freq_unit="MHz", basemap="haslam", data_unit="K" + ) with pytest.raises(RuntimeError): - gsm = BaseSkyModel('TEST_GSM', 'not_a_file.madeup', - freq_unit='MHz', basemap='haslam', data_unit='K') + gsm = BaseSkyModel( + "TEST_GSM", + "not_a_file.madeup", + freq_unit="MHz", + basemap="haslam", + data_unit="K", + ) with pytest.raises(RuntimeError): - gsm = BaseSkyModel('TEST_GSM', Path(__file__).absolute(), - freq_unit='MHz', basemap='haslam', data_unit='K') + gsm = BaseSkyModel( + "TEST_GSM", + Path(__file__).absolute(), + freq_unit="MHz", + basemap="haslam", + data_unit="K", + ) + if __name__ == "__main__": - test_base_skymodel_init() \ No newline at end of file + test_base_skymodel_init() diff --git a/tests/test_download_data.py b/tests/test_download_data.py new file mode 100644 index 0000000..01c8d6a --- /dev/null +++ b/tests/test_download_data.py @@ -0,0 +1,7 @@ +from pygdsm import download_map_data + +def test_download_map_data(): + download_map_data() + +if __name__ == "__main__": + download_map_data() \ No newline at end of file diff --git a/tests/test_gsm.py b/tests/test_gsm.py index aeb2142..c854cb8 100644 --- a/tests/test_gsm.py +++ b/tests/test_gsm.py @@ -5,20 +5,25 @@ Tests for GSM module. """ -from pygdsm import GlobalSkyModel -from pygdsm.component_data import GSM2008_TEST_DATA - import os import time -import pytest -import numpy as np + import h5py import healpy as hp +import numpy as np +import pytest from astropy.coordinates import SkyCoord +from astropy.utils.data import download_file + +from pygdsm import GlobalSkyModel +from pygdsm.component_data import GSM2008_TEST_DATA_URL + +# Download test data as needed (stored in astropy cache) +GSM2008_TEST_DATA = download_file(GSM2008_TEST_DATA_URL, cache=True, show_progress=True) def test_gsm_generate(): - """ Test maps generate successfully, and that output shapes are correct. """ + """Test maps generate successfully, and that output shapes are correct.""" freq = 40 gsm = GlobalSkyModel() @@ -26,20 +31,20 @@ def test_gsm_generate(): assert map.shape == (3145728,) freqs = [40, 80, 120, 160] - gsm = GlobalSkyModel(interpolation='cubic') + gsm = GlobalSkyModel(interpolation="cubic") map = gsm.generate(freqs) assert map.shape == (4, 3145728) freqs = np.linspace(1, 90, 10) - freq_unit = 'GHz' - for map_id in ['5deg', 'haslam', 'wmap']: - gsm = GlobalSkyModel(basemap=map_id, interpolation='pchip', freq_unit=freq_unit) + freq_unit = "GHz" + for map_id in ["5deg", "haslam", "wmap"]: + gsm = GlobalSkyModel(basemap=map_id, interpolation="pchip", freq_unit=freq_unit) map = gsm.generate(freqs) assert map.shape == (10, 3145728) def test_compare_to_gsm(): - """ Compare output of python version to fortran version. + """Compare output of python version to fortran version. Compares against output of original GSM. Note that the interpolation functions used in this and in the original differ. So, the outputs @@ -55,13 +60,13 @@ def test_compare_to_gsm(): Note: Because each output """ - gsm = GlobalSkyModel(basemap='haslam', interpolation='pchip') - gsm_fortran = h5py.File(GSM2008_TEST_DATA, 'r') + gsm = GlobalSkyModel(basemap="haslam", interpolation="pchip") + gsm_fortran = h5py.File(GSM2008_TEST_DATA, "r") for freq in (10, 100, 408, 1000, 1420, 2326, 10000, 23000, 40000, 90000, 94000): print("\nComparing at %i MHz..." % freq) - a = gsm_fortran['map_%imhz.out' % freq][:] + a = gsm_fortran["map_%imhz.out" % freq][:] b = gsm.generate(freq) - frac_avg = np.average(np.abs(1 - a/b)) + frac_avg = np.average(np.abs(1 - a / b)) if frac_avg > 0.01: print("FORTRAN: ", a) print("PYTHON: ", b) @@ -70,26 +75,25 @@ def test_compare_to_gsm(): def test_set_methods(): - gsm = GlobalSkyModel() gsm.generate(40) - #gsm.view() + # gsm.view() - gsm.set_basemap('haslam') + gsm.set_basemap("haslam") gsm.view() - gsm.set_basemap('5deg') + gsm.set_basemap("5deg") gsm.view() - gsm.set_basemap('wmap') + gsm.set_basemap("wmap") gsm.view() - gsm.set_freq_unit('GHz') + gsm.set_freq_unit("GHz") gsm.view() def test_speed(): - gsm = GlobalSkyModel(basemap='haslam') + gsm = GlobalSkyModel(basemap="haslam") t1 = time.time() for ff in np.linspace(10, 10000, 100): @@ -104,46 +108,51 @@ def test_write_fits(): gsm.write_fits("test_write_fits.fits") d_fits = hp.read_map("test_write_fits.fits") - d_gsm = gsm.generated_map_data + d_gsm = gsm.generated_map_data assert d_fits.shape == d_gsm.shape assert np.allclose(d_fits, d_gsm) os.remove("test_write_fits.fits") + def test_cmb_removal(): - g = GlobalSkyModel(freq_unit='MHz', include_cmb=False) + g = GlobalSkyModel(freq_unit="MHz", include_cmb=False) sky_no_cmb = g.generate(400) - g = GlobalSkyModel(freq_unit='MHz', include_cmb=True) + g = GlobalSkyModel(freq_unit="MHz", include_cmb=True) sky_with_cmb = g.generate(400) T_cmb = (sky_with_cmb - sky_no_cmb).mean() print(T_cmb) assert np.isclose(T_cmb, 2.725) + def test_get_sky_temperature(): - gc = SkyCoord(0, 0, unit='deg', frame='galactic') + gc = SkyCoord(0, 0, unit="deg", frame="galactic") freqs = (50, 100, 150) g = GlobalSkyModel() T = g.get_sky_temperature(gc, freqs) - T_gold = np.array([258368.7463149 , 50466.45058671, 18968.12555978]) + T_gold = np.array([258368.7463149, 50466.45058671, 18968.12555978]) assert np.allclose(T, T_gold) + def test_stupid_values(): with pytest.raises(RuntimeError): - g = GlobalSkyModel(basemap='haslamalan') + g = GlobalSkyModel(basemap="haslamalan") with pytest.raises(RuntimeError): - g = GlobalSkyModel(interpolation='linear') + g = GlobalSkyModel(interpolation="linear") with pytest.raises(RuntimeError): g = GlobalSkyModel() g.generate(9999999999) + def test_set_interpolation_method(): - g = GlobalSkyModel(interpolation='pchip') - assert g.interpolation_method == 'pchip' - g.set_interpolation_method('cubic') - assert g.interpolation_method == 'cubic' + g = GlobalSkyModel(interpolation="pchip") + assert g.interpolation_method == "pchip" + g.set_interpolation_method("cubic") + assert g.interpolation_method == "cubic" + if __name__ == "__main__": test_stupid_values() diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index ca9cf70..6d7f7bd 100644 --- a/tests/test_gsm2016.py +++ b/tests/test_gsm2016.py @@ -1,31 +1,30 @@ import sys + sys.path.append("/workdata/pygdsm") -import pylab as plt -import healpy as hp -import numpy as np from datetime import datetime +import healpy as hp +import numpy as np +import pylab as plt import pytest -from pygdsm import GlobalSkyModel16, GSMObserver16 -from pygdsm import GlobalSkyModel, GSMObserver +from pygdsm import GlobalSkyModel, GlobalSkyModel16, GSMObserver, GSMObserver16 def test_compare_gsm_to_old(): - g = GlobalSkyModel16(freq_unit='GHz', resolution='hi', data_unit='TCMB') + g = GlobalSkyModel16(freq_unit="GHz", resolution="hi", data_unit="TCMB") d = g.generate(0.408) g.view() - g_old = GlobalSkyModel(freq_unit='GHz', basemap='haslam') + g_old = GlobalSkyModel(freq_unit="GHz", basemap="haslam") d_old = g_old.generate(0.408) g_old.view() def test_observer_test(): - # Setup observatory location - in this case, Parkes Australia - (latitude, longitude, elevation) = ('-32.998370', '148.263659', 100) + (latitude, longitude, elevation) = ("-32.998370", "148.263659", 100) ov = GSMObserver16() ov.lon = longitude ov.lat = latitude @@ -45,18 +44,19 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + def test_interp(): f = np.arange(40, 80, 5) - for interp in ('pchip', 'cubic'): + for interp in ("pchip", "cubic"): for SkyModel in (GlobalSkyModel, GlobalSkyModel16): - name = str(SkyModel).strip("<>").split('.')[-1].strip("' ") - gsm = SkyModel(freq_unit='MHz', interpolation=interp) + name = str(SkyModel).strip("<>").split(".")[-1].strip("' ") + gsm = SkyModel(freq_unit="MHz", interpolation=interp) d = gsm.generate(f) sky_spec = d.mean(axis=1) fit = np.poly1d(np.polyfit(f, sky_spec, 5))(f) - plt.plot(f, sky_spec - fit, label=f'{name}: {interp}') + plt.plot(f, sky_spec - fit, label=f"{name}: {interp}") plt.xlabel("Frequency [MHz]") plt.ylabel("Residual [K]") @@ -64,29 +64,33 @@ def test_interp(): plt.show() - def test_gsm_opts(): - g = GlobalSkyModel16(freq_unit='GHz', resolution='hi', data_unit='TCMB') + g = GlobalSkyModel16(freq_unit="GHz", resolution="hi", data_unit="TCMB") d = g.generate(0.408) - g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='MJysr') + g = GlobalSkyModel16(freq_unit="GHz", resolution="lo", data_unit="MJysr") d = g.generate(0.408) - g = GlobalSkyModel16(freq_unit='MHz', resolution='lo', data_unit='TRJ') + g = GlobalSkyModel16(freq_unit="MHz", resolution="lo", data_unit="TRJ") d = g.generate(408) with pytest.raises(RuntimeError): g.generate(5e12) with pytest.raises(RuntimeError): - g = GlobalSkyModel16(resolution='oh_hai') + g = GlobalSkyModel16(resolution="oh_hai") with pytest.raises(RuntimeError): - g = GlobalSkyModel16(data_unit='furlongs/fortnight') + g = GlobalSkyModel16(data_unit="furlongs/fortnight") + def test_cmb_removal(): - g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='TCMB', include_cmb=False) + g = GlobalSkyModel16( + freq_unit="GHz", resolution="lo", data_unit="TCMB", include_cmb=False + ) sky_no_cmb = g.generate(400) - g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='TCMB', include_cmb=True) - sky_with_cmb = g.generate(400) + g = GlobalSkyModel16( + freq_unit="GHz", resolution="lo", data_unit="TCMB", include_cmb=True + ) + sky_with_cmb = g.generate(400) T_cmb = (sky_with_cmb - sky_no_cmb).mean() print(T_cmb) @@ -97,4 +101,4 @@ def test_cmb_removal(): test_compare_gsm_to_old() test_observer_test() test_gsm_opts() - test_interp() \ No newline at end of file + test_interp() diff --git a/tests/test_gsm_observer.py b/tests/test_gsm_observer.py index 4382d42..5ed762b 100644 --- a/tests/test_gsm_observer.py +++ b/tests/test_gsm_observer.py @@ -1,15 +1,17 @@ -from pygdsm import GSMObserver, GSMObserver16, LFSMObserver -import pylab as plt -import healpy as hp +import os from datetime import datetime + +import healpy as hp import numpy as np -import os +import pylab as plt from astropy.time import Time +from pygdsm import GSMObserver, GSMObserver16, LFSMObserver + + def test_gsm_observer(show=False): - """ Test GSMObserver() is working - """ - (latitude, longitude, elevation) = ('37.2', '-118.2', 1222) + """Test GSMObserver() is working""" + (latitude, longitude, elevation) = ("37.2", "-118.2", 1222) ov = GSMObserver() ov.lon = longitude ov.lat = latitude @@ -21,7 +23,7 @@ def test_gsm_observer(show=False): if show: plt.show() - (latitude, longitude, elevation) = ('37.2', '-118.2', 1222) + (latitude, longitude, elevation) = ("37.2", "-118.2", 1222) ov = GSMObserver16() ov.lon = longitude ov.lat = latitude @@ -33,7 +35,7 @@ def test_gsm_observer(show=False): if show: plt.show() - (latitude, longitude, elevation) = ('37.2', '-118.2', 1222) + (latitude, longitude, elevation) = ("37.2", "-118.2", 1222) ov = LFSMObserver() ov.lon = longitude ov.lat = latitude @@ -45,48 +47,49 @@ def test_gsm_observer(show=False): if show: plt.show() + def test_observed_mollview(): - """ Generate animated maps of observing coverage over 24 hours """ + """Generate animated maps of observing coverage over 24 hours""" - (latitude, longitude, elevation) = ('37.2', '-118.2', 1222) + (latitude, longitude, elevation) = ("37.2", "-118.2", 1222) ov = GSMObserver() ov.lon = longitude ov.lat = latitude 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") for ii in range(0, 24, 4): ov.date = datetime(2000, 1, 1, ii, 0) ov.generate(50) sky = ov.view_observed_gsm(logged=True, show=False, min=9, max=20) - plt.savefig('generated_sky/galactic-%02d.png' % ii) + plt.savefig("generated_sky/galactic-%02d.png" % ii) plt.close() - hp.mollview(sky, coord=['G', 'E'], min=9, max=20) - plt.savefig('generated_sky/ecliptic-%02d.png' % ii) + hp.mollview(sky, coord=["G", "E"], min=9, max=20) + plt.savefig("generated_sky/ecliptic-%02d.png" % ii) plt.close() - hp.mollview(sky, coord=['G', 'C'], min=9, max=20) - plt.savefig('generated_sky/equatorial-%02d.png' % ii) + hp.mollview(sky, coord=["G", "C"], min=9, max=20) + plt.savefig("generated_sky/equatorial-%02d.png" % ii) plt.close() ov.view(logged=True, show=False, min=9, max=20) - plt.savefig('generated_sky/ortho-%02d.png' % ii) + plt.savefig("generated_sky/ortho-%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/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(): - """ Test generating without frequency argument """ - (latitude, longitude, elevation) = ('37.2', '-118.2', 1222) + """Test generating without frequency argument""" + (latitude, longitude, elevation) = ("37.2", "-118.2", 1222) ov = GSMObserver() ov.lon = longitude ov.lat = latitude @@ -107,6 +110,7 @@ def test_generate_with_and_without_args(): ov.generate(obstime=now, freq=53) ov.generate(obstime=now, freq=52) + if __name__ == "__main__": test_gsm_observer(show=True) test_observed_mollview() diff --git a/tests/test_haslam.py b/tests/test_haslam.py index a8aaf26..12a28d8 100644 --- a/tests/test_haslam.py +++ b/tests/test_haslam.py @@ -1,28 +1,30 @@ import sys + import numpy as np sys.path.append("/workdata/pygdsm") -import pylab as plt -import healpy as hp from datetime import datetime -from pygdsm import HaslamSkyModel, HaslamObserver, GSMObserver +import healpy as hp +import pylab as plt -def test_compare_gsm_to_old(): +from pygdsm import GSMObserver, HaslamObserver, HaslamSkyModel - gl = HaslamSkyModel(freq_unit='MHz') + +def test_compare_gsm_to_old(): + gl = HaslamSkyModel(freq_unit="MHz") dl = gl.generate(408) gl.view() import pylab as plt + plt.show() def test_observer_test(): - # Setup observatory location - in this case, Parkes Australia - (latitude, longitude, elevation) = ('-32.998370', '148.263659', 100) + (latitude, longitude, elevation) = ("-32.998370", "148.263659", 100) ov = HaslamObserver() ov.lon = longitude ov.lat = latitude @@ -42,11 +44,12 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + def test_cmb_removal(): - g = HaslamSkyModel(freq_unit='MHz', include_cmb=False) + g = HaslamSkyModel(freq_unit="MHz", include_cmb=False) sky_no_cmb = g.generate(400) - g = HaslamSkyModel(freq_unit='MHz', include_cmb=True) - sky_with_cmb = g.generate(400) + g = HaslamSkyModel(freq_unit="MHz", include_cmb=True) + sky_with_cmb = g.generate(400) T_cmb = (sky_with_cmb - sky_no_cmb).mean() print(T_cmb) diff --git a/tests/test_init.py b/tests/test_init.py index cf25ffe..e87e68a 100644 --- a/tests/test_init.py +++ b/tests/test_init.py @@ -5,24 +5,34 @@ Tests for GSM init commands """ -from pygdsm import GlobalSkyModel, GlobalSkyModel08, GlobalSkyModel16, LowFrequencySkyModel, HaslamSkyModel -from pygdsm import GSMObserver, GSMObserver16, LFSMObserver, HaslamObserver -from pygdsm import init_gsm, init_observer +from pygdsm import ( + GlobalSkyModel, + GlobalSkyModel08, + GlobalSkyModel16, + GSMObserver, + GSMObserver16, + HaslamObserver, + HaslamSkyModel, + LFSMObserver, + LowFrequencySkyModel, + init_gsm, + init_observer, +) def test_init(): - assert isinstance(init_gsm('gsm'), GlobalSkyModel) - assert isinstance(init_gsm('gsm08'), GlobalSkyModel) - assert isinstance(init_gsm('gsm16'), GlobalSkyModel16) - assert isinstance(init_gsm('lfsm'), LowFrequencySkyModel) - assert isinstance(init_gsm('haslam'), HaslamSkyModel) + assert isinstance(init_gsm("gsm"), GlobalSkyModel) + assert isinstance(init_gsm("gsm08"), GlobalSkyModel) + assert isinstance(init_gsm("gsm16"), GlobalSkyModel16) + assert isinstance(init_gsm("lfsm"), LowFrequencySkyModel) + assert isinstance(init_gsm("haslam"), HaslamSkyModel) - assert isinstance(init_observer('gsm'), GSMObserver) - assert isinstance(init_observer('gsm08'), GSMObserver) - assert isinstance(init_observer('gsm16'), GSMObserver16) - assert isinstance(init_observer('lfsm'), LFSMObserver) - assert isinstance(init_observer('haslam'), HaslamObserver) + assert isinstance(init_observer("gsm"), GSMObserver) + assert isinstance(init_observer("gsm08"), GSMObserver) + assert isinstance(init_observer("gsm16"), GSMObserver16) + assert isinstance(init_observer("lfsm"), LFSMObserver) + assert isinstance(init_observer("haslam"), HaslamObserver) if __name__ == "__main__": - test_init() \ No newline at end of file + test_init() diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index 1735c27..bda238d 100644 --- a/tests/test_lfsm.py +++ b/tests/test_lfsm.py @@ -1,30 +1,37 @@ import sys + import numpy as np sys.path.append("/workdata/pygdsm") -import pylab as plt -import healpy as hp from datetime import datetime -from pygdsm import GlobalSkyModel16, GlobalSkyModel, LowFrequencySkyModel -from pygdsm import GSMObserver16, GSMObserver, LFSMObserver +import healpy as hp +import pylab as plt +from pygdsm import ( + GlobalSkyModel, + GlobalSkyModel16, + GSMObserver, + GSMObserver16, + LFSMObserver, + LowFrequencySkyModel, +) -def test_compare_gsm_to_old(): - gl = LowFrequencySkyModel(freq_unit='MHz') +def test_compare_gsm_to_old(): + gl = LowFrequencySkyModel(freq_unit="MHz") dl = gl.generate(408) gl.view() import pylab as plt + plt.show() def test_observer_test(): - # Setup observatory location - in this case, Parkes Australia - (latitude, longitude, elevation) = ('-32.998370', '148.263659', 100) + (latitude, longitude, elevation) = ("-32.998370", "148.263659", 100) ov = LFSMObserver() ov.lon = longitude ov.lat = latitude @@ -44,16 +51,18 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + def test_cmb_removal(): - g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=False) + g = LowFrequencySkyModel(freq_unit="MHz", include_cmb=False) sky_no_cmb = g.generate(400) - g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=True) - sky_with_cmb = g.generate(400) + g = LowFrequencySkyModel(freq_unit="MHz", include_cmb=True) + sky_with_cmb = g.generate(400) T_cmb = (sky_with_cmb - sky_no_cmb).mean() print(T_cmb) assert np.isclose(T_cmb, 2.725) + if __name__ == "__main__": test_compare_gsm_to_old() # test_observer_test() diff --git a/tests/test_utils.py b/tests/test_utils.py index 8308c2c..c6806c8 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,16 +1,17 @@ -import numpy as np import healpy as hp +import numpy as np + from pygdsm.utils import hpix2sky, sky2hpix def test_pix2sky(): - """ Small test routine for converting healpix pixel_id to and from SkyCoords """ + """Small test routine for converting healpix pixel_id to and from SkyCoords""" NSIDE = 32 pix = np.arange(hp.nside2npix(NSIDE)) - sc = hpix2sky(NSIDE, pix) + sc = hpix2sky(NSIDE, pix) pix_roundtrip = sky2hpix(NSIDE, sc) assert np.allclose(pix, pix_roundtrip) if __name__ == "__main__": - test_pix2sky() \ No newline at end of file + test_pix2sky()