Merge pull request #14 from fact-project/exposure

Exposure

Author: maxnoe

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

@@ -0,0 +1,88 @@
+import astropy.units as u
+from astropy.coordinates import Angle, SkyCoord
+from astropy import wcs
+from regions import CircleSkyRegion
+import numpy as np
+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.
+    '''
+    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()
+
+    exposure_time = len(delta_s) * scale
+    loc = min(delta_s)
+
+    # this percentile is somewhat abritrary but seems to work well.
+    live_time_fraction = 1 - expon.ppf(0.1, loc=loc, scale=scale)
+
+    return (exposure_time * live_time_fraction) * u.s
+
+
+@u.quantity_input(ra=u.hourangle, dec=u.deg, fov=u.deg)
+def build_exposure_regions(pointing_coords, fov=4.5 * u.deg):
+    '''
+    Takes a list of pointing positions and a field of view and returns
+    the unique pointing positions and the astropy.regions.
+
+    For an observation with N wobble positions this will return N unique
+    pointing positions and N circular regions.
+    '''
+    unique_pointing_positions = SkyCoord(
+        ra=np.unique(pointing_coords.ra),
+        dec=np.unique(pointing_coords.dec)
+    )
+
+    regions = [CircleSkyRegion(
+        center=pointing,
+        radius=Angle(fov) / 2
+    ) for pointing in unique_pointing_positions]
+
+    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, wcs=None, shape=(1000, 1000)):
+    '''
+    Takes pointing coordinates for each event and the corresponding timestamp.
+    Returns a masked array containing the estimated exposure time in hours
+    and a WCS object so the mask can be plotted. 
+    '''
+    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)
+
+    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])
+        times.append(exposure_time)
+
+    masks = [r.to_pixel(wcs).to_mask().to_image(shape) for r in regions]
+    cutout = sum(masks).astype(bool)
+    mask = sum([w.to('h').value * m for w, m in zip(times, masks)])
+
+    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

@@ -0,0 +1,54 @@
+import os
+import pytest
+import fact.io as fio
+import pandas as pd
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+from irf import estimate_exposure_time, build_exposure_map
+
+
+FIXTURE_DIR = os.path.join(
+    os.path.dirname(os.path.realpath(__file__)),
+    'test_files',
+)
+
+
+@pytest.fixture
+def events():
+    return fio.read_data(os.path.join(FIXTURE_DIR, 'crab_dl3_sample.hdf5'), key='events')
+
+
+@pytest.fixture
+def runs():
+    return fio.read_data(os.path.join(FIXTURE_DIR, 'crab_dl3_sample.hdf5'), key='runs')
+
+
+
+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)
+
+    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)
+
+    live_time_crab = runs.ontime.sum() / 3600  # to hours
+
+    expected_live_time_sample = n_test_sample / n_crab_sample * live_time_crab
+
+    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)