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

Passed

Push — master ( 1d1675...226afd )

by Axel

created 2021-11-19 20:14 UTC

test_make_map_no_count() A

↳ Parent: gammapy.makers.tests.test_map

Complexity

Conditions

Size

Total Lines	9
Code Lines	8

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
eloc	8
nop	1
dl	0
loc	9
rs	10
c	0
b	0
f	0

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import pytest
import numpy as np
from numpy.testing import assert_allclose
from regions import CircleSkyRegion
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.time import Time
from gammapy.data import GTI, DataStore, EventList, Observation
from gammapy.datasets import MapDataset
from gammapy.irf import EDispKernelMap, EDispMap, PSFMap
from gammapy.makers import MapDatasetMaker, SafeMaskMaker, FoVBackgroundMaker
from gammapy.maps import Map, MapAxis, WcsGeom, HpxGeom
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)

@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={},
        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(
        SkyCoord("0 deg", "0 deg")
    ).pdf_matrix
    pdfmatrix_postinterp = dataset.edisp.get_edisp_kernel(
        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)


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

gammapy / gammapy

Push — master ( 1d1675...226afd )

test_make_map_no_count() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like