Code for making maps of the sky, projecting HEALPix arrays
@ipynb_doc
def fits_image_demo():
"""
### FITSimage demo
This loads the FITS image file from astropy, {filename},
loads it with this class, then displays it:
{fig}
An object of the class is a functor: the value at the GC is {gcval:.0f}.
"""
from astropy.utils.data import get_pkg_data_filename
fullfilename = get_pkg_data_filename('galactic_center/gc_2mass_k.fits')
filename = Path(fullfilename).name
fi2 = FITSimage(fullfilename)
fi2.plot();
fig = plt.gcf()
gcval = fi2(SkyCoord(0,0,unit='deg',frame='galactic'))
return locals()
fits_image_demo()
hpm = None
@ipynb_doc
def all_sky_count_map(bmin=4, nside=128, pct=(50, 99.98)):
"""
## Count plots with the full data set
Time to extract data: {elapsed}
Data set: from {start} to {stop}.
Plotted with nside={nside}, band index >= {bmin}, {photon_count:,} total photons
Note that value limits, {vlim}, correspond to specified percentiles {pct}.
### All-sky: AIT projection
{fig1}
### Individual sources: ZEA projections
{fig2}
"""
global hpm # save this for subsquent demos
# from wtlike.exposure import WeightedAeff
dvall = DataView(nside=nside, bmin=bmin)
a, b = dvall.time_range
start, stop = UTC(a)[:10], UTC(b)[:10]
with Timer() as elapsed:
cntmap = dvall.count_map() #beam_window=WeightedAeff(Config(), spectrum=None).beam_window() )
photon_count = sum(cntmap)
vlim = np.percentile(cntmap, pct).round()
hpm = HPmap(cntmap, f'MJD {dvall.time_range}\nbmin {bmin}',
unit=f'counts per nside={nside} pixel')
fig1 = hpm.ait_plot(log=True, tick_labels=False, pixelsize=0.2, figsize=(20,8.),
vlim=vlim )
fig2 = plt.figure(figsize=(10,10))
plt.subplots_adjust(wspace=0.2, hspace=0.3)
kw = dict( fig=fig2, size=20, colorbar=False, log=True, vmin=1e2, vmax=1e4, cmap='inferno')
for i,name in enumerate(['Geminga','3C 273', 'BL Lac', '3C 454.3']):
hpm.zea_plot(name, title=name, axpos=221+i, **kw)
return locals()
if valid: all_sky_count_map()
Count plots with the full data set
Time to extract data: elapsed time: 104.3s (1.7 min)
Data set: from 2008-08-05 to 2022-10-27.
Plotted with nside=128, band index >= 4, 122,350,773 total photons
Note that value limits, [ 279. 25031.], correspond to specified percentiles (50, 99.98).
All-sky: AIT projection
Individual sources: ZEA projections
@ipynb_doc
def squaremap_demo(name, hpmap, size, vmin=10, vmax=1000):
r"""
### Demonstrate SquareWCS
This looks at a region of the all-sky count map generated above.
Start with a ${size}^\circ$ plot in Galactic coordinates, centered on the source "{name}".
The `SquareWCS.plot` function is called with a callback specified by a kwarg `annotator`.
Here the `annotate` function identifies Geminga and the Crab.
{fig1}
Repeat, using the equatorial "fk5" frame for this one.
{fig2}
"""
sc = findsource(name, gal=True)
swcs = SquareWCS(sc, size=90)
def annotate(ax, frame):
def source_lonlat(name, gal=True):
sky = findsource(name, gal=gal)
if gal: return sky.galactic.l.deg, sky.galactic.b.deg
return sky.ra.deg, sky.dec.deg
tf =ax.get_transform(frame) # need the transformation from pixel to world
for name in ['Geminga', 'Crab pulsar']:
glb = source_lonlat(name, gal=frame=='galactic')
ax.scatter(*glb , transform=tf, s=500, edgecolor='black', facecolor='none')
ax.text(*glb, ' '+name, color='black', transform=tf);
cb_kw = dict(ticks = [0.1, 1, 10, 100, 1000],
ticklabels='0.1 1 10 100 1000'.split(), shrink=0.8)
fig1 = swcs.plot(hpmap, log=True, vmin=vmin, vmax=vmax, unit='counts/pixel',
cb_kw=cb_kw, annotator = annotate)
swcsc = SquareWCS(sc.fk5, size=90)
fig2 = swcsc.plot(hpmap, log=True, vmin=vmin, vmax=vmax, unit='counts/pixel',
cb_kw=cb_kw, annotator=annotate)
return locals()
if valid and hpm is not None: squaremap_demo('Crab pulsar', hpmap=hpm.map, size=90, fignum=1)
Demonstrate SquareWCS
This looks at a region of the all-sky count map generated above.
Start with a $90^\circ$ plot in Galactic coordinates, centered on the source "Crab pulsar".
The SquareWCS.plot function is called with a callback specified by a kwarg annotator.
Here the annotate function identifies Geminga and the Crab.
Repeat, using the equatorial "fk5" frame for this one.
@ipynb_doc
def livetime_demo():
"""
## Livetime demo
Demonstrate the code to generate a 30-day map of livetime.
To generate an animation:
```
from wtlike.config import FermiMonth
p = Path('./livetime')
for k,t in enumerate(FermiMonth()):
filename = p/f'month_{k:03d}.png'
fig = make_ltmap(t, vmin=0, vmax=200, colorbar=False);
fig.text(0.22,0.14, f'month {k+1}', fontsize=12)
fig.savefig(filename, bbox_inches='tight')
fig.clear()
```
The following command will make an animated gif
```
ffmpeg -i livetime/month_%03d.png -loop -framerate 1 monthly_livetime.gif
```
The framerate (frames/sec) does not seem to work
__First and last months__
The orange band is the path of the Sun during the 30 days.
{fig1}
{fig2}
"""
from wtlike.config import FermiMonth
fm = FermiInterval(30)
def get_livetime(time_range):
return DataView(time_range).livetime_map(sigma=2)
def make_month(k):
fig = draw_map(get_livetime, fm[k], fignum=k,
vmin=0, vmax=200, colorbar=False);
month = k+1 if k>=0 else len(fm)+k+1
fig.text(0.22,0.14, f'month {month:3}', fontsize=12)
return fig
fig1 = make_month(0)
fig2 = make_month(-1)
return locals()
if valid: livetime_demo()
Livetime demo
Demonstrate the code to generate a 30-day map of livetime. To generate an animation:
from wtlike.config import FermiMonth
p = Path('./livetime')
for k,t in enumerate(FermiMonth()):
filename = p/f'month_{k:03d}.png'
fig = make_ltmap(t, vmin=0, vmax=200, colorbar=False);
fig.text(0.22,0.14, f'month {k+1}', fontsize=12)
fig.savefig(filename, bbox_inches='tight')
fig.clear()The following command will make an animated gif
ffmpeg -i livetime/month_%03d.png -loop -framerate 1 monthly_livetime.gifThe framerate (frames/sec) does not seem to work
First and last months The orange band is the path of the Sun during the 30 days.
