Skip to content

Improve data download from datacentral.org.au #30

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
telegraphic merged 3 commits into from
Dec 21, 2024
Merged

Improve data download from datacentral.org.au #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0ff605a
Improve data download from datacentral.org.au
telegraphic Dec 21, 2024
e14fc37
Fix get_cache_dir_path import location
telegraphic Dec 21, 2024
3db1c0e
Use get_cache_dir()
telegraphic Dec 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 26 additions & 26 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)
Expand All @@ -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).

Expand All @@ -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.

Expand All @@ -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
----------
Expand Down Expand Up @@ -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):
Expand All @@ -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).

49 changes: 21 additions & 28 deletions pygdsm/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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,
Expand All @@ -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()
case "haslam":
return HaslamObserver()
50 changes: 27 additions & 23 deletions pygdsm/base_observer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
"""
Expand All @@ -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))

Expand All @@ -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
----------
Expand Down Expand Up @@ -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
Expand All @@ -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
----------
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Loading
Loading