gammapy.makers.tests.test_map.test_dataset_hawc() - Code Metrics - Inspection of "Merge pull request #4111 from adonath/codestyle" - gammapy/gammapy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — master ( 2176c3...4ab053 )

by Axel

created 2022-09-26 21:49 UTC

gammapy.makers.tests.test_map.test_dataset_hawc() B

↳ Parent: gammapy.makers.tests.test_map

Complexity

Conditions

Size

Total Lines	66
Code Lines	42

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	2
eloc	42
nop	0
dl	0
loc	66
rs	8.872
c	0
b	0
f	0

How to fix Long Method

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import pytest
import numpy as np
from numpy.testing import assert_allclose
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.time import Time
from regions import CircleSkyRegion
from gammapy.data import (
    GTI,
    DataStore,
    EventList,
    HDUIndexTable,
    Observation,
    ObservationTable,
)
from gammapy.datasets import MapDataset
from gammapy.datasets.map import RAD_AXIS_DEFAULT
from gammapy.irf import EDispKernelMap, EDispMap, PSFMap
from gammapy.makers import FoVBackgroundMaker, MapDatasetMaker, SafeMaskMaker
from gammapy.maps import HpxGeom, Map, MapAxis, WcsGeom
from gammapy.utils.testing import requires_data, requires_dependency


@pytest.fixture(scope="session")
def observations():
    data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
    obs_id = [110380, 111140]
    return data_store.get_observations(obs_id)


def geom(ebounds, binsz=0.5):
    skydir = SkyCoord(0, -1, unit="deg", frame="galactic")
    energy_axis = MapAxis.from_edges(ebounds, name="energy", unit="TeV", interp="log")
    return WcsGeom.create(
        skydir=skydir, binsz=binsz, width=(10, 5), frame="galactic", axes=[energy_axis]
    )


@pytest.fixture(scope="session")
def geom_config_hpx():
    energy_axis = MapAxis.from_energy_bounds("0.5 TeV", "30 TeV", nbin=3)

    energy_axis_true = MapAxis.from_energy_bounds(
        "0.3 TeV", "30 TeV", nbin=3, per_decade=True, name="energy_true"
    )

    geom_hpx = HpxGeom.create(
        binsz=0.1, frame="galactic", axes=[energy_axis], region="DISK(0, 0, 5.)"
    )
    return {"geom": geom_hpx, "energy_axis_true": energy_axis_true}


@requires_data()
@pytest.mark.parametrize(
    "pars",
    [
        {
            # Default, same e_true and reco
            "geom": geom(ebounds=[0.1, 1, 10]),
            "e_true": None,
            "counts": 34366,
            "exposure": 9.995376e08,
            "exposure_image": 3.99815e11,
            "background": 27989.05,
            "binsz_irf": 0.5,
            "migra": None,
        },
        {
            # Test single energy bin
            "geom": geom(ebounds=[0.1, 10]),
            "e_true": None,
            "counts": 34366,
            "exposure": 5.843302e08,
            "exposure_image": 1.16866e11,
            "background": 30424.451,
            "binsz_irf": 0.5,
            "migra": None,
        },
        {
            # Test single energy bin with exclusion mask
            "geom": geom(ebounds=[0.1, 10]),
            "e_true": None,
            "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])),
            "counts": 34366,
            "exposure": 5.843302e08,
            "exposure_image": 1.16866e11,
            "background": 30424.451,
            "binsz_irf": 0.5,
            "migra": None,
        },
        {
            # Test for different e_true and e_reco bins
            "geom": geom(ebounds=[0.1, 1, 10]),
            "e_true": MapAxis.from_edges(
                [0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
            ),
            "counts": 34366,
            "exposure": 9.951827e08,
            "exposure_image": 5.971096e11,
            "background": 28760.283,
            "background_oversampling": 2,
            "binsz_irf": 0.5,
            "migra": None,
        },
        {
            # Test for different e_true and e_reco and spatial bins
            "geom": geom(ebounds=[0.1, 1, 10]),
            "e_true": MapAxis.from_edges(
                [0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
            ),
            "counts": 34366,
            "exposure": 9.951827e08,
            "exposure_image": 5.971096e11,
            "background": 28760.283,
            "background_oversampling": 2,
            "binsz_irf": 1.0,
            "migra": None,
        },
        {
            # Test for different e_true and e_reco and use edispmap
            "geom": geom(ebounds=[0.1, 1, 10]),
            "e_true": MapAxis.from_edges(
                [0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
            ),
            "counts": 34366,
            "exposure": 9.951827e08,
            "exposure_image": 5.971096e11,
            "background": 28760.283,
            "background_oversampling": 2,
            "binsz_irf": 0.5,
            "migra": MapAxis.from_edges(
                np.linspace(0.0, 3.0, 100), name="migra", unit=""
            ),
        },
    ],
)
def test_map_maker(pars, observations):
    stacked = MapDataset.create(
        geom=pars["geom"],
        energy_axis_true=pars["e_true"],
        binsz_irf=pars["binsz_irf"],
        migra_axis=pars["migra"],
    )

    maker = MapDatasetMaker(background_oversampling=pars.get("background_oversampling"))
    safe_mask_maker = SafeMaskMaker(methods=["offset-max"], offset_max="2 deg")

    for obs in observations:
        cutout = stacked.cutout(position=obs.pointing_radec, width="4 deg")
        dataset = maker.run(cutout, obs)
        dataset = safe_mask_maker.run(dataset, obs)
        stacked.stack(dataset)

    counts = stacked.counts
    assert counts.unit == ""
    assert_allclose(counts.data.sum(), pars["counts"], rtol=1e-5)

    exposure = stacked.exposure
    assert exposure.unit == "m2 s"
    assert_allclose(exposure.data.mean(), pars["exposure"], rtol=3e-3)

    background = stacked.npred_background()
    assert background.unit == ""
    assert_allclose(background.data.sum(), pars["background"], rtol=1e-4)

    image_dataset = stacked.to_image()

    counts = image_dataset.counts
    assert counts.unit == ""
    assert_allclose(counts.data.sum(), pars["counts"], rtol=1e-4)

    exposure = image_dataset.exposure
    assert exposure.unit == "m2 s"
    assert_allclose(exposure.data.sum(), pars["exposure_image"], rtol=1e-3)

    background = image_dataset.npred_background()
    assert background.unit == ""
    assert_allclose(background.data.sum(), pars["background"], rtol=1e-4)


@requires_data()
def test_map_maker_obs(observations):
    # Test for different spatial geoms and etrue, ereco bins

    geom_reco = geom(ebounds=[0.1, 1, 10])
    e_true = MapAxis.from_edges(
        [0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
    )

    reference = MapDataset.create(
        geom=geom_reco, energy_axis_true=e_true, binsz_irf=1.0
    )

    maker_obs = MapDatasetMaker()

    map_dataset = maker_obs.run(reference, observations[0])
    assert map_dataset.counts.geom == geom_reco
    assert map_dataset.npred_background().geom == geom_reco
    assert isinstance(map_dataset.edisp, EDispKernelMap)
    assert map_dataset.edisp.edisp_map.data.shape == (3, 2, 5, 10)
    assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10)
    assert map_dataset.psf.psf_map.data.shape == (3, 66, 5, 10)
    assert map_dataset.psf.exposure_map.data.shape == (3, 1, 5, 10)
    assert_allclose(map_dataset.gti.time_delta, 1800.0 * u.s)


@requires_data()
def test_map_maker_obs_with_migra(observations):
    # Test for different spatial geoms and etrue, ereco bins
    migra = MapAxis.from_edges(np.linspace(0, 2.0, 50), unit="", name="migra")
    geom_reco = geom(ebounds=[0.1, 1, 10])
    e_true = MapAxis.from_edges(
        [0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
    )

    reference = MapDataset.create(
        geom=geom_reco, energy_axis_true=e_true, migra_axis=migra, binsz_irf=1.0
    )

    maker_obs = MapDatasetMaker()

    map_dataset = maker_obs.run(reference, observations[0])
    assert map_dataset.counts.geom == geom_reco
    assert isinstance(map_dataset.edisp, EDispMap)
    assert map_dataset.edisp.edisp_map.data.shape == (3, 49, 5, 10)
    assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10)


@requires_data()
def test_make_meta_table(observations):
    maker_obs = MapDatasetMaker()
    map_dataset_meta_table = maker_obs.make_meta_table(observation=observations[0])

    assert_allclose(map_dataset_meta_table["RA_PNT"], 267.68121338)
    assert_allclose(map_dataset_meta_table["DEC_PNT"], -29.6075)
    assert_allclose(map_dataset_meta_table["OBS_ID"], 110380)
    assert map_dataset_meta_table["OBS_MODE"] == "POINTING"


@requires_data()
def test_make_map_no_count(observations):
    dataset = MapDataset.create(geom((0.1, 1, 10)))
    maker_obs = MapDatasetMaker(selection=["exposure"])
    map_dataset = maker_obs.run(dataset, observation=observations[0])

    assert map_dataset.counts is not None
    assert_allclose(map_dataset.counts.data, 0)
    assert map_dataset.counts.geom == dataset.counts.geom


@requires_data()
@requires_dependency("healpy")
def test_map_dataset_maker_hpx(geom_config_hpx, observations):
    reference = MapDataset.create(**geom_config_hpx, binsz_irf=5 * u.deg)

    maker = MapDatasetMaker()
    safe_mask_maker = SafeMaskMaker(
        offset_max="2.5 deg", methods=["aeff-default", "offset-max"]
    )

    dataset = maker.run(reference, observation=observations[0])
    dataset = safe_mask_maker.run(dataset, observation=observations[0]).to_masked()

    assert_allclose(dataset.counts.data.sum(), 4264)
    assert_allclose(dataset.background.data.sum(), 2964.5369, rtol=1e-5)
    assert_allclose(dataset.exposure.data[4, 1000], 5.987e09, rtol=1e-4)

    coords = SkyCoord([0, 3], [0, 0], frame="galactic", unit="deg")
    coords = {"skycoord": coords, "energy": 1 * u.TeV}
    assert_allclose(dataset.mask_safe.get_by_coord(coords), [True, False])

    kernel = dataset.edisp.get_edisp_kernel()

    assert_allclose(kernel.data.sum(axis=1)[3], 1, rtol=0.01)


def test_interpolate_map_dataset():
    energy = MapAxis.from_energy_bounds("1 TeV", "300 TeV", nbin=5, name="energy")
    energy_true = MapAxis.from_nodes(
        np.logspace(-1, 3, 20), name="energy_true", interp="log", unit="TeV"
    )

    # make dummy map IRFs
    geom_allsky = WcsGeom.create(
        npix=(5, 3), proj="CAR", binsz=60, axes=[energy], skydir=(0, 0)
    )
    geom_allsky_true = geom_allsky.drop("energy").to_cube([energy_true])

    # background
    geom_background = WcsGeom.create(
        skydir=(0, 0), width=(5, 5), binsz=0.2 * u.deg, axes=[energy]
    )
    value = 30
    bkg_map = Map.from_geom(geom_background, unit="")
    bkg_map.data = value * np.ones(bkg_map.data.shape)

    # effective area - with a gradient that also depends on energy
    aeff_map = Map.from_geom(geom_allsky_true, unit="cm2 s")
    ra_arr = np.arange(aeff_map.data.shape[1])
    dec_arr = np.arange(aeff_map.data.shape[2])
    for i in np.arange(aeff_map.data.shape[0]):
        aeff_map.data[i, :, :] = (
            (i + 1) * 10 * np.meshgrid(dec_arr, ra_arr)[0]
            + 10 * np.meshgrid(dec_arr, ra_arr)[1]
            + 10
        )
    aeff_map.meta["TELESCOP"] = "HAWC"

    # psf map
    width = 0.2 * u.deg
    rad_axis = MapAxis.from_nodes(np.linspace(0, 2, 50), name="rad", unit="deg")
    psfMap = PSFMap.from_gauss(energy_true, rad_axis, width)

    # edispmap
    edispmap = EDispKernelMap.from_gauss(
        energy, energy_true, sigma=0.1, bias=0.0, geom=geom_allsky
    )

    # events and gti
    nr_ev = 10
    ev_t = Table()
    gti_t = Table()

    ev_t["EVENT_ID"] = np.arange(nr_ev)
    ev_t["TIME"] = nr_ev * [Time("2011-01-01 00:00:00", scale="utc", format="iso")]
    ev_t["RA"] = np.linspace(-1, 1, nr_ev) * u.deg
    ev_t["DEC"] = np.linspace(-1, 1, nr_ev) * u.deg
    ev_t["ENERGY"] = np.logspace(0, 2, nr_ev) * u.TeV

    gti_t["START"] = [Time("2010-12-31 00:00:00", scale="utc", format="iso")]
    gti_t["STOP"] = [Time("2011-01-02 00:00:00", scale="utc", format="iso")]

    events = EventList(ev_t)
    gti = GTI(gti_t)

    # define observation
    obs = Observation(
        obs_id=0,
        obs_info={"RA_PNT": 0.0, "DEC_PNT": 0.0},
        gti=gti,
        aeff=aeff_map,
        edisp=edispmap,
        psf=psfMap,
        bkg=bkg_map,
        events=events,
        obs_filter=None,
    )

    # define analysis geometry
    geom_target = WcsGeom.create(
        skydir=(0, 0), width=(5, 5), binsz=0.1 * u.deg, axes=[energy]
    )

    maker = MapDatasetMaker(
        selection=["exposure", "counts", "background", "edisp", "psf"]
    )
    dataset = MapDataset.create(
        geom=geom_target, energy_axis_true=energy_true, rad_axis=rad_axis, name="test"
    )
    dataset = maker.run(dataset, obs)

    # test counts
    assert dataset.counts.data.sum() == nr_ev

    # test background
    assert np.floor(np.sum(dataset.npred_background().data)) == np.sum(bkg_map.data)
    coords_bg = {"skycoord": SkyCoord("0 deg", "0 deg"), "energy": energy.center[0]}
    assert_allclose(
        dataset.npred_background().get_by_coord(coords_bg)[0], 7.5, atol=1e-4
    )

    # test effective area
    coords_aeff = {
        "skycoord": SkyCoord("0 deg", "0 deg"),
        "energy_true": energy_true.center[0],
    }
    assert_allclose(
        aeff_map.get_by_coord(coords_aeff)[0],
        dataset.exposure.interp_by_coord(coords_aeff)[0],
        atol=1e-3,
    )

    # test edispmap
    pdfmatrix_preinterp = edispmap.get_edisp_kernel(
        position=SkyCoord("0 deg", "0 deg")
    ).pdf_matrix
    pdfmatrix_postinterp = dataset.edisp.get_edisp_kernel(
        position=SkyCoord("0 deg", "0 deg")
    ).pdf_matrix
    assert_allclose(pdfmatrix_preinterp, pdfmatrix_postinterp, atol=1e-7)

    # test psfmap
    geom_psf = geom_target.drop("energy").to_cube([energy_true])
    psfkernel_preinterp = psfMap.get_psf_kernel(
        position=SkyCoord("0 deg", "0 deg"), geom=geom_psf, max_radius=2 * u.deg
    ).data
    psfkernel_postinterp = dataset.psf.get_psf_kernel(
        position=SkyCoord("0 deg", "0 deg"), geom=geom_psf, max_radius=2 * u.deg
    ).data
    assert_allclose(psfkernel_preinterp, psfkernel_postinterp, atol=1e-4)


@requires_data()
@pytest.mark.xfail
def test_minimal_datastore():
    """ "Check that a standard analysis runs on a minimal datastore"""

    energy_axis = MapAxis.from_energy_bounds(
        1, 10, nbin=3, per_decade=False, unit="TeV", name="energy"
    )
    geom = WcsGeom.create(
        skydir=(83.633, 22.014),
        binsz=0.5,
        width=(2, 2),
        frame="icrs",
        proj="CAR",
        axes=[energy_axis],
    )

    data_store = DataStore.from_dir("$GAMMAPY_DATA/tests/minimal_datastore")

    observations = data_store.get_observations()
    maker = MapDatasetMaker()
    offset_max = 2.3 * u.deg
    maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=offset_max)
    circle = CircleSkyRegion(
        center=SkyCoord("83.63 deg", "22.14 deg"), radius=0.2 * u.deg
    )
    exclusion_mask = ~geom.region_mask(regions=[circle])
    maker_fov = FoVBackgroundMaker(method="fit", exclusion_mask=exclusion_mask)

    stacked = MapDataset.create(geom=geom, name="crab-stacked")
    for obs in observations:
        dataset = maker.run(stacked, obs)
        dataset = maker_safe_mask.run(dataset, obs)
        dataset = maker_fov.run(dataset)
        stacked.stack(dataset)

    assert_allclose(stacked.exposure.data.sum(), 6.01909e10)
    assert_allclose(stacked.counts.data.sum(), 1446)
    assert_allclose(stacked.background.data.sum(), 1445.9841)


@requires_data()
def test_dataset_hawc():
    # create the energy reco axis
    energy_axis = MapAxis.from_edges(
        [1.00, 1.78, 3.16, 5.62, 10.0, 17.8, 31.6, 56.2, 100, 177, 316] * u.TeV,
        name="energy",
        interp="log",
    )

    # and energy true axis
    energy_axis_true = MapAxis.from_energy_bounds(
        1e-3, 1e4, nbin=140, unit="TeV", name="energy_true"
    )

    # create a geometry around the Crab location
    geom = WcsGeom.create(
        skydir=SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs"),
        width=3 * u.deg,
        axes=[energy_axis],
        binsz=0.1,
    )

    maker = MapDatasetMaker(
        selection=["counts", "background", "exposure", "edisp", "psf"]
    )
    safemask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)

    results = {}
    results["GP"] = [6.53623241669e16, 58, 0.72202391]
    results["NN"] = [6.57154247837e16, 62, 0.76743538]

    for which in ["GP", "NN"]:

        # paths and file names
        data_path = "$GAMMAPY_DATA/hawc/crab_events_pass4/"
        hdu_filename = "hdu-index-table-" + which + "-Crab.fits.gz"
        obs_filename = "obs-index-table-" + which + "-Crab.fits.gz"

        # We want the last event lass for speed
        obs_table = ObservationTable.read(data_path + obs_filename)
        hdu_table = HDUIndexTable.read(data_path + hdu_filename, hdu=9)
        data_store = DataStore(hdu_table=hdu_table, obs_table=obs_table)

        observations = data_store.get_observations()

        # create empty dataset that will contain the data
        geom_exposure = geom.to_image().to_cube([energy_axis_true])
        geom_psf = geom.to_image().to_cube([RAD_AXIS_DEFAULT, energy_axis])
        geom_edisp = geom.to_cube([energy_axis_true])

        dataset_empty = MapDataset.from_geoms(
            geom=geom,
            name="nHit-9",
            geom_exposure=geom_exposure,
            geom_psf=geom_psf,
            geom_edisp=geom_edisp,
        )

        # run the maker
        dataset = maker.run(dataset_empty, observations[0])
        dataset.exposure.meta["livetime"] = "1 s"
        dataset = safemask_maker.run(dataset)

        assert_allclose(dataset.exposure.data.sum(), results[which][0])
        assert_allclose(dataset.counts.data.sum(), results[which][1])
        assert_allclose(dataset.background.data.sum(), results[which][2])


1			# Licensed under a 3-clause BSD style license - see LICENSE.rst
2			import pytest
3			import numpy as np
4			from numpy.testing import assert_allclose
5			import astropy.units as u
6			from astropy.coordinates import SkyCoord
7			from astropy.table import Table
8			from astropy.time import Time
9			from regions import CircleSkyRegion
10			from gammapy.data import (
11			GTI,
12			DataStore,
13			EventList,
14			HDUIndexTable,
15			Observation,
16			ObservationTable,
17			)
18			from gammapy.datasets import MapDataset
19			from gammapy.datasets.map import RAD_AXIS_DEFAULT
20			from gammapy.irf import EDispKernelMap, EDispMap, PSFMap
21			from gammapy.makers import FoVBackgroundMaker, MapDatasetMaker, SafeMaskMaker
22			from gammapy.maps import HpxGeom, Map, MapAxis, WcsGeom
23			from gammapy.utils.testing import requires_data, requires_dependency
24
25
26			@pytest.fixture(scope="session")
27			def observations():
28			data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
29			obs_id = [110380, 111140]
30			return data_store.get_observations(obs_id)
31
32
33			def geom(ebounds, binsz=0.5):
34			skydir = SkyCoord(0, -1, unit="deg", frame="galactic")
35			energy_axis = MapAxis.from_edges(ebounds, name="energy", unit="TeV", interp="log")
36			return WcsGeom.create(
37			skydir=skydir, binsz=binsz, width=(10, 5), frame="galactic", axes=[energy_axis]
38			)
39
40
41			@pytest.fixture(scope="session")
42			def geom_config_hpx():
43			energy_axis = MapAxis.from_energy_bounds("0.5 TeV", "30 TeV", nbin=3)
44
45			energy_axis_true = MapAxis.from_energy_bounds(
46			"0.3 TeV", "30 TeV", nbin=3, per_decade=True, name="energy_true"
47			)
48
49			geom_hpx = HpxGeom.create(
50			binsz=0.1, frame="galactic", axes=[energy_axis], region="DISK(0, 0, 5.)"
51			)
52			return {"geom": geom_hpx, "energy_axis_true": energy_axis_true}
53
54
55			@requires_data()
56			@pytest.mark.parametrize(
57			"pars",
58			[
59			{
60			# Default, same e_true and reco
61			"geom": geom(ebounds=[0.1, 1, 10]),
62			"e_true": None,
63			"counts": 34366,
64			"exposure": 9.995376e08,
65			"exposure_image": 3.99815e11,
66			"background": 27989.05,
67			"binsz_irf": 0.5,
68			"migra": None,
69			},
70			{
71			# Test single energy bin
72			"geom": geom(ebounds=[0.1, 10]),
73			"e_true": None,
74			"counts": 34366,
75			"exposure": 5.843302e08,
76			"exposure_image": 1.16866e11,
77			"background": 30424.451,
78			"binsz_irf": 0.5,
79			"migra": None,
80			},
81			{
82			# Test single energy bin with exclusion mask
83			"geom": geom(ebounds=[0.1, 10]),
84			"e_true": None,
85			"exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])),
86			"counts": 34366,
87			"exposure": 5.843302e08,
88			"exposure_image": 1.16866e11,
89			"background": 30424.451,
90			"binsz_irf": 0.5,
91			"migra": None,
92			},
93			{
94			# Test for different e_true and e_reco bins
95			"geom": geom(ebounds=[0.1, 1, 10]),
96			"e_true": MapAxis.from_edges(
97			[0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
98			),
99			"counts": 34366,
100			"exposure": 9.951827e08,
101			"exposure_image": 5.971096e11,
102			"background": 28760.283,
103			"background_oversampling": 2,
104			"binsz_irf": 0.5,
105			"migra": None,
106			},
107			{
108			# Test for different e_true and e_reco and spatial bins
109			"geom": geom(ebounds=[0.1, 1, 10]),
110			"e_true": MapAxis.from_edges(
111			[0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
112			),
113			"counts": 34366,
114			"exposure": 9.951827e08,
115			"exposure_image": 5.971096e11,
116			"background": 28760.283,
117			"background_oversampling": 2,
118			"binsz_irf": 1.0,
119			"migra": None,
120			},
121			{
122			# Test for different e_true and e_reco and use edispmap
123			"geom": geom(ebounds=[0.1, 1, 10]),
124			"e_true": MapAxis.from_edges(
125			[0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
126			),
127			"counts": 34366,
128			"exposure": 9.951827e08,
129			"exposure_image": 5.971096e11,
130			"background": 28760.283,
131			"background_oversampling": 2,
132			"binsz_irf": 0.5,
133			"migra": MapAxis.from_edges(
134			np.linspace(0.0, 3.0, 100), name="migra", unit=""
135			),
136			},
137			],
138			)
139			def test_map_maker(pars, observations):
140			stacked = MapDataset.create(
141			geom=pars["geom"],
142			energy_axis_true=pars["e_true"],
143			binsz_irf=pars["binsz_irf"],
144			migra_axis=pars["migra"],
145			)
146
147			maker = MapDatasetMaker(background_oversampling=pars.get("background_oversampling"))
148			safe_mask_maker = SafeMaskMaker(methods=["offset-max"], offset_max="2 deg")
149
150			for obs in observations:
151			cutout = stacked.cutout(position=obs.pointing_radec, width="4 deg")
152			dataset = maker.run(cutout, obs)
153			dataset = safe_mask_maker.run(dataset, obs)
154			stacked.stack(dataset)
155
156			counts = stacked.counts
157			assert counts.unit == ""
158			assert_allclose(counts.data.sum(), pars["counts"], rtol=1e-5)
159
160			exposure = stacked.exposure
161			assert exposure.unit == "m2 s"
162			assert_allclose(exposure.data.mean(), pars["exposure"], rtol=3e-3)
163
164			background = stacked.npred_background()
165			assert background.unit == ""
166			assert_allclose(background.data.sum(), pars["background"], rtol=1e-4)
167
168			image_dataset = stacked.to_image()
169
170			counts = image_dataset.counts
171			assert counts.unit == ""
172			assert_allclose(counts.data.sum(), pars["counts"], rtol=1e-4)
173
174			exposure = image_dataset.exposure
175			assert exposure.unit == "m2 s"
176			assert_allclose(exposure.data.sum(), pars["exposure_image"], rtol=1e-3)
177
178			background = image_dataset.npred_background()
179			assert background.unit == ""
180			assert_allclose(background.data.sum(), pars["background"], rtol=1e-4)
181
182
183			@requires_data()
184			def test_map_maker_obs(observations):
185			# Test for different spatial geoms and etrue, ereco bins
186
187			geom_reco = geom(ebounds=[0.1, 1, 10])
188			e_true = MapAxis.from_edges(
189			[0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
190			)
191
192			reference = MapDataset.create(
193			geom=geom_reco, energy_axis_true=e_true, binsz_irf=1.0
194			)
195
196			maker_obs = MapDatasetMaker()
197
198			map_dataset = maker_obs.run(reference, observations[0])
199			assert map_dataset.counts.geom == geom_reco
200			assert map_dataset.npred_background().geom == geom_reco
201			assert isinstance(map_dataset.edisp, EDispKernelMap)
202			assert map_dataset.edisp.edisp_map.data.shape == (3, 2, 5, 10)
203			assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10)
204			assert map_dataset.psf.psf_map.data.shape == (3, 66, 5, 10)
205			assert map_dataset.psf.exposure_map.data.shape == (3, 1, 5, 10)
206			assert_allclose(map_dataset.gti.time_delta, 1800.0 * u.s)
207
208
209			@requires_data()
210			def test_map_maker_obs_with_migra(observations):
211			# Test for different spatial geoms and etrue, ereco bins
212			migra = MapAxis.from_edges(np.linspace(0, 2.0, 50), unit="", name="migra")
213			geom_reco = geom(ebounds=[0.1, 1, 10])
214			e_true = MapAxis.from_edges(
215			[0.1, 0.5, 2.5, 10.0], name="energy_true", unit="TeV", interp="log"
216			)
217
218			reference = MapDataset.create(
219			geom=geom_reco, energy_axis_true=e_true, migra_axis=migra, binsz_irf=1.0
220			)
221
222			maker_obs = MapDatasetMaker()
223
224			map_dataset = maker_obs.run(reference, observations[0])
225			assert map_dataset.counts.geom == geom_reco
226			assert isinstance(map_dataset.edisp, EDispMap)
227			assert map_dataset.edisp.edisp_map.data.shape == (3, 49, 5, 10)
228			assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10)
229
230
231			@requires_data()
232			def test_make_meta_table(observations):
233			maker_obs = MapDatasetMaker()
234			map_dataset_meta_table = maker_obs.make_meta_table(observation=observations[0])
235
236			assert_allclose(map_dataset_meta_table["RA_PNT"], 267.68121338)
237			assert_allclose(map_dataset_meta_table["DEC_PNT"], -29.6075)
238			assert_allclose(map_dataset_meta_table["OBS_ID"], 110380)
239			assert map_dataset_meta_table["OBS_MODE"] == "POINTING"
240
241
242			@requires_data()
243			def test_make_map_no_count(observations):
244			dataset = MapDataset.create(geom((0.1, 1, 10)))
245			maker_obs = MapDatasetMaker(selection=["exposure"])
246			map_dataset = maker_obs.run(dataset, observation=observations[0])
247
248			assert map_dataset.counts is not None
249			assert_allclose(map_dataset.counts.data, 0)
250			assert map_dataset.counts.geom == dataset.counts.geom
251
252
253			@requires_data()
254			@requires_dependency("healpy")
255			def test_map_dataset_maker_hpx(geom_config_hpx, observations):
256			reference = MapDataset.create(*geom_config_hpx, binsz_irf=5 u.deg)
257
258			maker = MapDatasetMaker()
259			safe_mask_maker = SafeMaskMaker(
260			offset_max="2.5 deg", methods=["aeff-default", "offset-max"]
261			)
262
263			dataset = maker.run(reference, observation=observations[0])
264			dataset = safe_mask_maker.run(dataset, observation=observations[0]).to_masked()
265
266			assert_allclose(dataset.counts.data.sum(), 4264)
267			assert_allclose(dataset.background.data.sum(), 2964.5369, rtol=1e-5)
268			assert_allclose(dataset.exposure.data[4, 1000], 5.987e09, rtol=1e-4)
269
270			coords = SkyCoord([0, 3], [0, 0], frame="galactic", unit="deg")
271			coords = {"skycoord": coords, "energy": 1 * u.TeV}
272			assert_allclose(dataset.mask_safe.get_by_coord(coords), [True, False])
273
274			kernel = dataset.edisp.get_edisp_kernel()
275
276			assert_allclose(kernel.data.sum(axis=1)[3], 1, rtol=0.01)
277
278
279			def test_interpolate_map_dataset():
280			energy = MapAxis.from_energy_bounds("1 TeV", "300 TeV", nbin=5, name="energy")
281			energy_true = MapAxis.from_nodes(
282			np.logspace(-1, 3, 20), name="energy_true", interp="log", unit="TeV"
283			)
284
285			# make dummy map IRFs
286			geom_allsky = WcsGeom.create(
287			npix=(5, 3), proj="CAR", binsz=60, axes=[energy], skydir=(0, 0)
288			)
289			geom_allsky_true = geom_allsky.drop("energy").to_cube([energy_true])
290
291			# background
292			geom_background = WcsGeom.create(
293			skydir=(0, 0), width=(5, 5), binsz=0.2 * u.deg, axes=[energy]
294			)
295			value = 30
296			bkg_map = Map.from_geom(geom_background, unit="")
297			bkg_map.data = value * np.ones(bkg_map.data.shape)
298
299			# effective area - with a gradient that also depends on energy
300			aeff_map = Map.from_geom(geom_allsky_true, unit="cm2 s")
301			ra_arr = np.arange(aeff_map.data.shape[1])
302			dec_arr = np.arange(aeff_map.data.shape[2])
303			for i in np.arange(aeff_map.data.shape[0]):
304			aeff_map.data[i, :, :] = (
305			(i + 1) * 10 * np.meshgrid(dec_arr, ra_arr)[0]
306			+ 10 * np.meshgrid(dec_arr, ra_arr)[1]
307			+ 10
308			)
309			aeff_map.meta["TELESCOP"] = "HAWC"
310
311			# psf map
312			width = 0.2 * u.deg
313			rad_axis = MapAxis.from_nodes(np.linspace(0, 2, 50), name="rad", unit="deg")
314			psfMap = PSFMap.from_gauss(energy_true, rad_axis, width)
315
316			# edispmap
317			edispmap = EDispKernelMap.from_gauss(
318			energy, energy_true, sigma=0.1, bias=0.0, geom=geom_allsky
319			)
320
321			# events and gti
322			nr_ev = 10
323			ev_t = Table()
324			gti_t = Table()
325
326			ev_t["EVENT_ID"] = np.arange(nr_ev)
327			ev_t["TIME"] = nr_ev * [Time("2011-01-01 00:00:00", scale="utc", format="iso")]
328			ev_t["RA"] = np.linspace(-1, 1, nr_ev) * u.deg
329			ev_t["DEC"] = np.linspace(-1, 1, nr_ev) * u.deg
330			ev_t["ENERGY"] = np.logspace(0, 2, nr_ev) * u.TeV
331
332			gti_t["START"] = [Time("2010-12-31 00:00:00", scale="utc", format="iso")]
333			gti_t["STOP"] = [Time("2011-01-02 00:00:00", scale="utc", format="iso")]
334
335			events = EventList(ev_t)
336			gti = GTI(gti_t)
337
338			# define observation
339			obs = Observation(
340			obs_id=0,
341			obs_info={"RA_PNT": 0.0, "DEC_PNT": 0.0},
342			gti=gti,
343			aeff=aeff_map,
344			edisp=edispmap,
345			psf=psfMap,
346			bkg=bkg_map,
347			events=events,
348			obs_filter=None,
349			)
350
351			# define analysis geometry
352			geom_target = WcsGeom.create(
353			skydir=(0, 0), width=(5, 5), binsz=0.1 * u.deg, axes=[energy]
354			)
355
356			maker = MapDatasetMaker(
357			selection=["exposure", "counts", "background", "edisp", "psf"]
358			)
359			dataset = MapDataset.create(
360			geom=geom_target, energy_axis_true=energy_true, rad_axis=rad_axis, name="test"
361			)
362			dataset = maker.run(dataset, obs)
363
364			# test counts
365			assert dataset.counts.data.sum() == nr_ev
366
367			# test background
368			assert np.floor(np.sum(dataset.npred_background().data)) == np.sum(bkg_map.data)
369			coords_bg = {"skycoord": SkyCoord("0 deg", "0 deg"), "energy": energy.center[0]}
370			assert_allclose(
371			dataset.npred_background().get_by_coord(coords_bg)[0], 7.5, atol=1e-4
372			)
373
374			# test effective area
375			coords_aeff = {
376			"skycoord": SkyCoord("0 deg", "0 deg"),
377			"energy_true": energy_true.center[0],
378			}
379			assert_allclose(
380			aeff_map.get_by_coord(coords_aeff)[0],
381			dataset.exposure.interp_by_coord(coords_aeff)[0],
382			atol=1e-3,
383			)
384
385			# test edispmap
386			pdfmatrix_preinterp = edispmap.get_edisp_kernel(
387			position=SkyCoord("0 deg", "0 deg")
388			).pdf_matrix
389			pdfmatrix_postinterp = dataset.edisp.get_edisp_kernel(
390			position=SkyCoord("0 deg", "0 deg")
391			).pdf_matrix
392			assert_allclose(pdfmatrix_preinterp, pdfmatrix_postinterp, atol=1e-7)
393
394			# test psfmap
395			geom_psf = geom_target.drop("energy").to_cube([energy_true])
396			psfkernel_preinterp = psfMap.get_psf_kernel(
397			position=SkyCoord("0 deg", "0 deg"), geom=geom_psf, max_radius=2 * u.deg
398			).data
399			psfkernel_postinterp = dataset.psf.get_psf_kernel(
400			position=SkyCoord("0 deg", "0 deg"), geom=geom_psf, max_radius=2 * u.deg
401			).data
402			assert_allclose(psfkernel_preinterp, psfkernel_postinterp, atol=1e-4)
403
404
405			@requires_data()
406			@pytest.mark.xfail
407			def test_minimal_datastore():
408			""" "Check that a standard analysis runs on a minimal datastore"""
409
410			energy_axis = MapAxis.from_energy_bounds(
411			1, 10, nbin=3, per_decade=False, unit="TeV", name="energy"
412			)
413			geom = WcsGeom.create(
414			skydir=(83.633, 22.014),
415			binsz=0.5,
416			width=(2, 2),
417			frame="icrs",
418			proj="CAR",
419			axes=[energy_axis],
420			)
421
422			data_store = DataStore.from_dir("$GAMMAPY_DATA/tests/minimal_datastore")
423
424			observations = data_store.get_observations()
425			maker = MapDatasetMaker()
426			offset_max = 2.3 * u.deg
427			maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=offset_max)
428			circle = CircleSkyRegion(
429			center=SkyCoord("83.63 deg", "22.14 deg"), radius=0.2 * u.deg
430			)
431			exclusion_mask = ~geom.region_mask(regions=[circle])
432			maker_fov = FoVBackgroundMaker(method="fit", exclusion_mask=exclusion_mask)
433
434			stacked = MapDataset.create(geom=geom, name="crab-stacked")
435			for obs in observations:
436			dataset = maker.run(stacked, obs)
437			dataset = maker_safe_mask.run(dataset, obs)
438			dataset = maker_fov.run(dataset)
439			stacked.stack(dataset)
440
441			assert_allclose(stacked.exposure.data.sum(), 6.01909e10)
442			assert_allclose(stacked.counts.data.sum(), 1446)
443			assert_allclose(stacked.background.data.sum(), 1445.9841)
444
445
446			@requires_data()
447			def test_dataset_hawc():
448			# create the energy reco axis
449			energy_axis = MapAxis.from_edges(
450			[1.00, 1.78, 3.16, 5.62, 10.0, 17.8, 31.6, 56.2, 100, 177, 316] * u.TeV,
451			name="energy",
452			interp="log",
453			)
454
455			# and energy true axis
456			energy_axis_true = MapAxis.from_energy_bounds(
457			1e-3, 1e4, nbin=140, unit="TeV", name="energy_true"
458			)
459
460			# create a geometry around the Crab location
461			geom = WcsGeom.create(
462			skydir=SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs"),
463			width=3 * u.deg,
464			axes=[energy_axis],
465			binsz=0.1,
466			)
467
468			maker = MapDatasetMaker(
469			selection=["counts", "background", "exposure", "edisp", "psf"]
470			)
471			safemask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
472
473			results = {}
474			results["GP"] = [6.53623241669e16, 58, 0.72202391]
475			results["NN"] = [6.57154247837e16, 62, 0.76743538]
476
477			for which in ["GP", "NN"]:
478
479			# paths and file names
480			data_path = "$GAMMAPY_DATA/hawc/crab_events_pass4/"
481			hdu_filename = "hdu-index-table-" + which + "-Crab.fits.gz"
482			obs_filename = "obs-index-table-" + which + "-Crab.fits.gz"
483
484			# We want the last event lass for speed
485			obs_table = ObservationTable.read(data_path + obs_filename)
486			hdu_table = HDUIndexTable.read(data_path + hdu_filename, hdu=9)
487			data_store = DataStore(hdu_table=hdu_table, obs_table=obs_table)
488
489			observations = data_store.get_observations()
490
491			# create empty dataset that will contain the data
492			geom_exposure = geom.to_image().to_cube([energy_axis_true])
493			geom_psf = geom.to_image().to_cube([RAD_AXIS_DEFAULT, energy_axis])
494			geom_edisp = geom.to_cube([energy_axis_true])
495
496			dataset_empty = MapDataset.from_geoms(
497			geom=geom,
498			name="nHit-9",
499			geom_exposure=geom_exposure,
500			geom_psf=geom_psf,
501			geom_edisp=geom_edisp,
502			)
503
504			# run the maker
505			dataset = maker.run(dataset_empty, observations[0])
506			dataset.exposure.meta["livetime"] = "1 s"
507			dataset = safemask_maker.run(dataset)
508
509			assert_allclose(dataset.exposure.data.sum(), results[which][0])
510			assert_allclose(dataset.counts.data.sum(), results[which][1])
511			assert_allclose(dataset.background.data.sum(), results[which][2])
512

gammapy / gammapy

Push — master ( 2176c3...4ab053 )

gammapy.makers.tests.test_map.test_dataset_hawc() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like