%matplotlib inline
import numpy as np
import pylab as plt
from pygdsm 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("gdsm_500mhz.fits")
setting the output map dtype to [dtype('float64')]
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)
Text(0, 0.5, 'Temperature [K]')
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 pygdsm 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)