add script and unit tests

Author: kbruegge

irf/__init__.py

@@ -1,4 +1,4 @@
 from .collection_area import collection_area
+from .exposure_map import estimate_exposure_time, build_exposure_map
 
-
-__all__ = ['collection_area']
+__all__ = ['collection_area', 'estimate_exposure_time', 'build_exposure_map']

irf/exposure_map.py

@@ -1,20 +1,17 @@
 import astropy.units as u
 from astropy.coordinates import Angle, SkyCoord
+from astropy import wcs
 from regions import CircleSkyRegion
-from astroquery.skyview import SkyView
-from astropy.wcs import WCS
-import matplotlib.pyplot as plt
-import fact.io as fio
 import numpy as np
-import pandas as pd
 from scipy.stats import expon
 
 
 def estimate_exposure_time(timestamps):
     '''
-    Takes numpy datetime64[ns] timestamps and returns an estimates of the exposure time in seconds
+    Takes numpy datetime64[ns] timestamps and returns an estimates of the exposure time in seconds.
     '''
     delta_s = np.diff(timestamps).astype(int) * float(1e-9)
+
     # take only events that came within 30 seconds or so
     delta_s = delta_s[delta_s < 30]
     scale = delta_s.mean()
@@ -43,71 +40,38 @@ def build_exposure_regions(pointing_coords, fov=4.5 * u.deg):
     return unique_pointing_positions, regions
 
 
+def _build_standard_wcs(image_center, shape, naxis=2, fov=9 * u.deg):
+    width, height = shape
+
+    w = wcs.WCS(naxis=2)
+    w.wcs.crpix = [width / 2 + 0.5, height / 2 + 0.5]
+    w.wcs.cdelt = np.array([-fov.value / width, fov.value / height])
+    w.wcs.crval = [image_center.ra.deg, image_center.dec.deg]
+    w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
+    w.wcs.radesys = 'FK5'
+    w.wcs.equinox = 2000.0
+    w.wcs.cunit = ['deg', 'deg']
+    w._naxis = [width, height]
+    return w
+
+
 @u.quantity_input(event_ra=u.hourangle, event_dec=u.deg, fov=u.deg)
-def build_exposure_map(pointing_coords, event_time, fov=4.5 * u.deg):
-    unique_pointing_positions, regions = build_exposure_regions(pointing_coords)
+def build_exposure_map(pointing_coords, event_time, fov=4.5 * u.deg, wcs=None, shape=(1000, 1000)):
 
-    img_center = SkyCoord(ra=pointing_coords.ra.mean(), dec=pointing_coords.dec.mean())
-    print('getting image')
-    hdu = SkyView.get_images(position=img_center, pixels=1000, survey=['DSS'], width=2 * fov, height=2 * fov)[0][0]
-    img = hdu.data
-    wcs = WCS(hdu.header)
+    if not wcs:
+        image_center = SkyCoord(ra=pointing_coords.ra.mean(), dec=pointing_coords.dec.mean())
+        wcs = _build_standard_wcs(image_center, shape, fov=2 * fov)
+
+    unique_pointing_positions, regions = build_exposure_regions(pointing_coords)
 
-    print('starting loop')
     times = []
     for p in unique_pointing_positions:
         m = (pointing_coords.ra == p.ra) & (pointing_coords.dec == p.dec)
         exposure_time = estimate_exposure_time(event_time[m])
-        print(exposure_time.to('h'))
         times.append(exposure_time)
 
-    print(sum(times).to('h'))
-
-    # import IPython; IPython.embed()
-    masks = [r.to_pixel(wcs).to_mask().to_image(img.shape) for r in regions]
-    #
+    masks = [r.to_pixel(wcs).to_mask().to_image(shape) for r in regions]
     cutout = sum(masks).astype(bool)
-    mask = sum([w * m for w, m in zip(times, masks)])
-
-    return np.ma.masked_array(mask.to('h'), mask=~cutout), wcs, img
-
-
-def plot_exposure(img, mask, wcs):
-    ax = plt.subplot(projection=wcs)
-
-    ax.imshow(img, cmap='gray')
-    d = ax.imshow(mask, alpha=0.7)
-    cb = plt.colorbar(d)
-    cb.set_label('Live Time / Hours')
-    ax.set_xlabel('Galactic Longitude')
-    ax.set_ylabel('Galactic Latitude')
-
-    crab = SkyCoord.from_name('Crab Nebula')
-    ax.scatter(crab.ra.deg, crab.dec.deg, transform=ax.get_transform('icrs'), label='Crab Nebula')
-    ax.legend()
-
-
-if __name__ == '__main__':
-    runs = fio.read_data('crab_dl3.hdf5', key='runs')
-    dl3 = fio.read_data('crab_dl3.hdf5', key='events')
-
-    data = pd.merge(runs, dl3, on=['run_id', 'night'] )
-
-    timestamps = pd.to_datetime(data.timestamp).values
-    total_ontime = estimate_exposure_time(timestamps)
-    print('total exposure time: {}'.format(total_ontime.to(u.h)))
-
-    ra = data.ra_prediction.values * u.hourangle
-    dec = data.dec_prediction.values * u.deg
-
-    ra_pointing = data.right_ascension.values * u.hourangle
-    dec_pointing = data.declination.values * u.deg
-
-    event_coords = SkyCoord(ra=ra, dec=dec)
-    pointing_coords = SkyCoord(ra=ra_pointing, dec=dec_pointing)
-
-    mask, wcs, img = build_exposure_map(pointing_coords, timestamps)
+    mask = sum([w.to('h').value * m for w, m in zip(times, masks)])
 
-    ax = plot_exposure(img, mask, wcs)
-    plt.savefig('exposure.pdf')
-    # plt.show()
+    return np.ma.masked_array(mask, mask=~cutout), wcs

irf/scripts/create_fact_exposure_map.py

@@ -0,0 +1,122 @@
+from irf import estimate_exposure_time, build_exposure_map
+import click
+import fact.io as fio
+import pandas as pd
+from astroquery.skyview import SkyView
+from astropy.wcs import WCS
+from astropy.visualization import ImageNormalize, ZScaleInterval, AsinhStretch
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+import matplotlib.pyplot as plt
+
+
+@u.quantity_input(fov=u.deg)
+def get_sdss_sky_image(img_center, fov=9 * u.deg, n_pix=1000):
+    '''
+    A small helper method which uses astroquery to get an image from the sdss.
+    This requires internet access and fails pretty often due to http timeouts.
+    '''
+    hdu = SkyView.get_images(
+        position=img_center,
+        pixels=n_pix,
+        survey=['DSS'],
+        width=fov,
+        height=fov)[0][0]
+
+    img = hdu.data
+    wcs = WCS(hdu.header)
+    return img, wcs
+
+
+@click.command()
+@click.argument(
+    'dl3_path',
+    type=click.Path(file_okay=True, dir_okay=False),
+)
+@click.argument(
+    'output_path',
+    type=click.Path(file_okay=True, dir_okay=False, exists=False),
+)
+@click.option(
+    '-n',
+    '--n_pix',
+    default=1000,
+    help='number of pixels along the edge of the produced image',
+)
+@click.option(
+    '-s',
+    '--source_name',
+    default=None,
+    help=
+    'If supplied, e.g. "Crab Nebula" will draw the position of that source into the image',
+)
+@click.option(
+    '--background/--no-background',
+    default=True,
+    help='If true, downloads SDSS image for backgrund image in the plot')
+def main(dl3_path, output_path, n_pix, source_name, background):
+    '''
+    Takes FACT dl3 output and plots a skymap which is being saved to the output_path.
+    '''
+    runs = fio.read_data(dl3_path, key='runs')
+    dl3 = fio.read_data(dl3_path, key='events')
+
+    data = pd.merge(runs, dl3, on=['run_id', 'night'])
+
+    timestamps = pd.to_datetime(data.timestamp).values
+    total_ontime = estimate_exposure_time(timestamps)
+    print('Total estimated exposure time: {}'.format(total_ontime.to(u.h)))
+
+    ra_pointing = data.right_ascension.values * u.hourangle
+    dec_pointing = data.declination.values * u.deg
+    pointing = SkyCoord(ra=ra_pointing, dec=dec_pointing)
+
+    img = None
+    wcs = None
+    if background:
+        img_center = SkyCoord(ra=pointing.ra.mean(), dec=pointing.dec.mean())
+        img, wcs = get_sdss_sky_image(
+            img_center=img_center, n_pix=n_pix, fov=9 * u.deg)
+
+    mask, wcs = build_exposure_map(pointing, timestamps, shape=(n_pix, n_pix))
+
+    ax = plot_exposure(mask, wcs, image=img)
+    if source_name:
+        source = SkyCoord.from_name('Crab Nebula')
+        ax.scatter(
+            source.ra.deg,
+            source.dec.deg,
+            transform=ax.get_transform('icrs'),
+            label=source_name,
+            s=10**2,
+            facecolors='none',
+            edgecolors='r',
+        )
+        ax.legend()
+
+    plt.savefig(output_path, dpi=200)
+
+
+def plot_exposure(mask, wcs, image=None):
+    plt.figure(figsize=(10, 10))
+    ax = plt.subplot(projection=wcs)
+
+    if image is not None:
+        norm = ImageNormalize(
+            image,
+            interval=ZScaleInterval(contrast=0.05),
+            stretch=AsinhStretch(a=0.2))
+        ax.imshow(image, cmap='gray', norm=norm, interpolation='nearest')
+
+    d = ax.imshow(mask, alpha=0.7)
+    cb = plt.colorbar(d)
+    cb.set_label('Live Time / Hours')
+    ax.set_xlabel('Galactic Longitude')
+    ax.set_ylabel('Galactic Latitude')
+
+    return ax
+
+
+if __name__ == '__main__':
+    main()

irf/tests/test_exposure_map.py

@@ -2,7 +2,11 @@
 import pytest
 import fact.io as fio
 import pandas as pd
-from irf import estimate_exposure_time
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+from irf import estimate_exposure_time, build_exposure_map
 
 
 FIXTURE_DIR = os.path.join(
@@ -22,15 +26,29 @@ def runs():
 
 
 
-def test_exposure_time(runs, events):
-        n_crab_sample = 687226
-        n_test_sample = len(events)
+def test_exposure_map(runs, events):
+    data = pd.merge(runs, events, on=['run_id', 'night'])
+
+    timestamps = pd.to_datetime(data.timestamp).values
+
+    ra_pointing = data.right_ascension.values * u.hourangle
+    dec_pointing = data.declination.values * u.deg
+    pointing = SkyCoord(ra=ra_pointing, dec=dec_pointing)
 
-        live_time_crab = runs.ontime.sum() / 3600  # to hours
+    mask, wcs = build_exposure_map(pointing, timestamps, fov=4.5 * u.deg)
+
+    assert mask.mask.any()  # assert whether some region is masked
+
+
+
+def test_exposure_time(runs, events):
+    n_crab_sample = 687226
+    n_test_sample = len(events)
 
-        expected_live_time_sample = n_test_sample / n_crab_sample * live_time_crab
+    live_time_crab = runs.ontime.sum() / 3600  # to hours
 
-        timestamps = pd.to_datetime(events.timestamp).values
-        exposure = estimate_exposure_time(timestamps).to('h').value
+    expected_live_time_sample = n_test_sample / n_crab_sample * live_time_crab
 
-        assert expected_live_time_sample == pytest.approx(exposure, 0.05)
+    timestamps = pd.to_datetime(events.timestamp).values
+    exposure = estimate_exposure_time(timestamps).to('h').value
+    assert expected_live_time_sample == pytest.approx(exposure, 0.05)