%matplotlib inline
Firstly, we import the pygsm
module:
from pygsm import GlobalSkyModel
We can very quickly get started with generating a sky map and plotting it:
gsm = GlobalSkyModel()
gsm.generate(500)
gsm.view(logged=True)
We might then want to write this to a FITS file (in healpix format):
gsm.write_fits("gsm_500mhz.fits")
WARNING: AstropyDeprecationWarning: The new_table function is deprecated and may be removed in a future version. Use :meth:`BinTableHDU.from_columns` for new BINARY tables or :meth:`TableHDU.from_columns` for new ASCII tables instead. [astropy.utils.misc] WARNING:astropy:AstropyDeprecationWarning: The new_table function is deprecated and may be removed in a future version. Use :meth:`BinTableHDU.from_columns` for new BINARY tables or :meth:`TableHDU.from_columns` for new ASCII tables instead. WARNING: AstropyDeprecationWarning: The use of header.update() to add new keywords to a header is deprecated. Instead, use either header.set() or simply `header[keyword] = value` or `header[keyword] = (value, comment)`. header.set() is only necessary to use if you also want to use the before/after keyword arguments. [astropy.io.fits.header] WARNING:astropy:AstropyDeprecationWarning: The use of header.update() to add new keywords to a header is deprecated. Instead, use either header.set() or simply `header[keyword] = value` or `header[keyword] = (value, comment)`. header.set() is only necessary to use if you also want to use the before/after keyword arguments.
There are also a few options we can play around with. For example:
gsm = GlobalSkyModel(freq_unit='GHz', interpolation='pchip', basemap='wmap')
gsm.generate(30) # Generate at 30 GHz
gsm.view(logged=True)
You can set things on the fly, too, once you've got a gsm
object:
gsm.set_basemap('5deg')
gsm.view(logged=True)
Just in case you're interested in spectral cubes, you can even pass a frequency range:
gsm = GlobalSkyModel()
freqs = np.linspace(10, 1000, 100)
map_cube = gsm.generate(freqs)
print map_cube.shape
plt.loglog(map_cube[:,0])
plt.loglog(map_cube[:,12345])
plt.xlabel("Frequency [MHz]")
plt.ylabel("Temperature [K]")
(100, 3145728)
<matplotlib.text.Text at 0x7fb2966d4950>
Be careful with this one though, as it can eat up a lot of memory!
A common task is to find out what the sky looks like at a given lat, long and date. PyGSM
gives a quick method to do this:
from pygsm import GSMObserver
from datetime import datetime
# Setup observatory location - in this case, Parkes Australia
(latitude, longitude, elevation) = ('-32.998370', '148.263659', 100)
ov = GSMObserver()
ov.lon = longitude
ov.lat = latitude
ov.elev = elevation
ov.date = datetime(2000, 1, 1, 23, 0)
Now generate a sky model and view an all-sky orthographic plot
ov.generate(1400)
d = ov.view(logged=True)
And we can view galactic Mollweide projection too, with below-horizon data blanked out:
d = ov.view_observed_gsm(logged=True)