Skip to content

Commit

Permalink
Merge pull request #63 from cta-observatory/ogadf
Browse files Browse the repository at this point in the history
Verify hdus using ogadf-schema
  • Loading branch information
maxnoe committed Sep 29, 2020
2 parents a0e2a24 + 24912f8 commit bc826fa
Show file tree
Hide file tree
Showing 2 changed files with 141 additions and 56 deletions.
190 changes: 135 additions & 55 deletions pyirf/io/tests/test_gadf.py
Expand Up @@ -8,109 +8,189 @@
import tempfile


def test_effective_area2d():
e_bins = np.geomspace(0.1, 100, 31) * u.TeV
migra_bins = np.linspace(0.2, 5, 101)
fov_bins = [0, 1, 2, 3] * u.deg
source_bins = np.linspace(0, 1, 101) * u.deg


@pytest.fixture
def aeff2d_hdus():
from pyirf.io import create_aeff2d_hdu

area = np.full((len(e_bins) - 1, len(fov_bins) - 1), 1e6) * u.m**2

hdus = [
create_aeff2d_hdu(area, e_bins, fov_bins, point_like=point_like)
for point_like in [True, False]
]

return area, hdus


@pytest.fixture
def edisp_hdus():
from pyirf.io import create_energy_dispersion_hdu

edisp = np.zeros((len(e_bins) - 1, len(migra_bins) - 1, len(fov_bins) - 1))
edisp[:, 50, :] = 1.0

hdus = [
create_energy_dispersion_hdu(
edisp, e_bins, migra_bins, fov_bins, point_like=point_like
)
for point_like in [True, False]
]

return edisp, hdus


@pytest.fixture
def psf_hdu():
from pyirf.io import create_psf_table_hdu
from pyirf.utils import cone_solid_angle

psf = np.zeros((len(e_bins) - 1, len(source_bins) - 1, len(fov_bins) - 1))
psf[:, 0, :] = 1
psf = psf / cone_solid_angle(source_bins[1])

hdu = create_psf_table_hdu(
psf, e_bins, source_bins, fov_bins, point_like=False
)
return psf, hdu


@pytest.fixture
def bg_hdu():
from pyirf.io import create_background_2d_hdu

background = np.column_stack([
np.geomspace(1e9, 1e3, 3),
np.geomspace(0.5e9, 0.5e3, 3),
np.geomspace(1e8, 1e2, 3),
]) * u.Unit('TeV-1 s-1 sr-1')

hdu = create_background_2d_hdu(background, e_bins, fov_bins)

return background, hdu


@pytest.fixture
def rad_max_hdu():
from pyirf.io import create_rad_max_hdu

rad_max = np.full((len(e_bins) - 1, len(fov_bins) - 1), 0.1) * u.deg
print(rad_max.shape)
hdu = create_rad_max_hdu(e_bins, fov_bins, rad_max)

return rad_max, hdu


def test_effective_area2d_gammapy(aeff2d_hdus):
'''Test our effective area is readable by gammapy'''
pytest.importorskip('gammapy')
from pyirf.io import create_aeff2d_hdu
from gammapy.irf import EffectiveAreaTable2D

e_bins = np.geomspace(0.1, 100, 31) * u.TeV
fov_bins = [0, 1, 2, 3] * u.deg
area = np.full((30, 3), 1e6) * u.m**2
area, hdus = aeff2d_hdus

for point_like in [True, False]:
for hdu in hdus:
with tempfile.NamedTemporaryFile(suffix='.fits') as f:
hdu = create_aeff2d_hdu(area, e_bins, fov_bins, point_like=point_like)

fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)

# test reading with gammapy works
aeff2d = EffectiveAreaTable2D.read(f.name)
assert u.allclose(area, aeff2d.data.data, atol=1e-16 * u.m**2)


def test_energy_dispersion():
def test_effective_area2d_schema(aeff2d_hdus):
'''Test our effective area is readable by gammapy'''
from ogadf_schema.irfs import AEFF_2D

_, hdus = aeff2d_hdus

for hdu in hdus:
AEFF_2D.validate_hdu(hdu)


def test_energy_dispersion_gammapy(edisp_hdus):
'''Test our energy dispersion is readable by gammapy'''
pytest.importorskip('gammapy')
from pyirf.io import create_energy_dispersion_hdu
from gammapy.irf import EnergyDispersion2D

e_bins = np.geomspace(0.1, 100, 31) * u.TeV
migra_bins = np.linspace(0.2, 5, 101)
fov_bins = [0, 1, 2, 3] * u.deg
edisp = np.zeros((30, 100, 3))
edisp[:, 50, :] = 1.0
edisp, hdus = edisp_hdus


for point_like in [True, False]:
for hdu in hdus:
with tempfile.NamedTemporaryFile(suffix='.fits') as f:
hdu = create_energy_dispersion_hdu(
edisp, e_bins, migra_bins, fov_bins, point_like=point_like
)

fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)

# test reading with gammapy works
edisp2d = EnergyDispersion2D.read(f.name, 'EDISP')
assert u.allclose(edisp, edisp2d.data.data, atol=1e-16)


def test_psf_table():
def test_energy_dispersion_schema(edisp_hdus):
from ogadf_schema.irfs import EDISP_2D

_, hdus = edisp_hdus

for hdu in hdus:
EDISP_2D.validate_hdu(hdu)


def test_psf_table_gammapy(psf_hdu):
'''Test our psf is readable by gammapy'''
pytest.importorskip('gammapy')
from pyirf.io import create_psf_table_hdu
from pyirf.utils import cone_solid_angle
from gammapy.irf import PSF3D

e_bins = np.geomspace(0.1, 100, 31) * u.TeV
source_bins = np.linspace(0, 1, 101) * u.deg
fov_bins = [0, 1, 2, 3] * u.deg
psf = np.zeros((30, 100, 3))
psf[:, 0, :] = 1
psf = psf / cone_solid_angle(source_bins[1])
psf, hdu = psf_hdu

for point_like in [True, False]:
with tempfile.NamedTemporaryFile(suffix='.fits') as f:
hdu = create_psf_table_hdu(
psf, e_bins, source_bins, fov_bins, point_like=point_like
)
with tempfile.NamedTemporaryFile(suffix='.fits') as f:
fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)

fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)
# test reading with gammapy works
psf3d = PSF3D.read(f.name, 'PSF')

# gammapy does not transpose psf when reading from fits,
# unlike how it handles effective area and edisp
# see https://github.com/gammapy/gammapy/issues/3025
assert u.allclose(psf, psf3d.psf_value.T, atol=1e-16 / u.sr)

# test reading with gammapy works
psf3d = PSF3D.read(f.name, 'PSF')

# gammapy does not transpose psf when reading from fits,
# unlike how it handles effective area and edisp
# see https://github.com/gammapy/gammapy/issues/3025
assert u.allclose(psf, psf3d.psf_value.T, atol=1e-16 / u.sr)
def test_psf_schema(psf_hdu):
from ogadf_schema.irfs import PSF_TABLE

_, hdu = psf_hdu
PSF_TABLE.validate_hdu(hdu)


# gammapy uses inconsistent axis order, should be fixed before gammapy 1.0
# see https://github.com/gammapy/gammapy/issues/2067
# TODO: remove xfail when this is fixed in gammapy and bump the required gammapy version
@pytest.mark.xfail
def test_background_2d():
def test_background_2d_gammapy(bg_hdu):
'''Test our background hdu is readable by gammapy'''

pytest.importorskip('gammapy')
from pyirf.io import create_background_2d_hdu
from gammapy.irf import Background2D

e_bins = np.geomspace(0.1, 100, 31) * u.TeV
fov_bins = [0, 1, 2, 3] * u.deg
background = np.column_stack([
np.geomspace(1e9, 1e3, 3),
np.geomspace(0.5e9, 0.5e3, 3),
np.geomspace(1e8, 1e2, 3),
]) * u.Unit('TeV-1 s-1 sr-1')
background, hdu = bg_hdu

with tempfile.NamedTemporaryFile(suffix='.fits') as f:
hdu = create_background_2d_hdu(background, e_bins, fov_bins)

fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)

# test reading with gammapy works
bg2d = Background2D.read(f.name, 'BACKGROUND')

assert u.allclose(background, bg2d.data.data, atol=1e-16 * u.Unit('TeV-1 s-1 sr-1'))


def test_background_2d_schema(bg_hdu):
from ogadf_schema.irfs import BKG_2D

_, hdu = bg_hdu
BKG_2D.validate_hdu(hdu)


def test_rad_max_schema(rad_max_hdu):
from ogadf_schema.irfs import RAD_MAX

_, hdu = rad_max_hdu
RAD_MAX.validate_hdu(hdu)
7 changes: 6 additions & 1 deletion setup.py
Expand Up @@ -13,7 +13,12 @@
"nbsphinx",
"uproot~=3.0",
],
"tests": ["pytest", "pytest-cov", "gammapy~=0.17"],
"tests": [
"pytest",
"pytest-cov",
"gammapy~=0.17",
"ogadf-schema~=0.2.3",
],
}

extras_require["all"] = extras_require["tests"] + extras_require["docs"]
Expand Down

0 comments on commit bc826fa

Please sign in to comment.