diff --git a/gammapy/cube/tests/test_background.py b/gammapy/cube/tests/test_background.py index 982c0f5de8..30d8cf05fc 100644 --- a/gammapy/cube/tests/test_background.py +++ b/gammapy/cube/tests/test_background.py @@ -35,12 +35,12 @@ def geom(map_type, ebounds): { "geom": geom(map_type="wcs", ebounds=[0.1, 1, 10]), "shape": (2, 3, 4), - "sum": 744.281309, + "sum": 744.268356, }, { "geom": geom(map_type="wcs", ebounds=[0.1, 10]), "shape": (1, 3, 4), - "sum": 799.760344, + "sum": 799.743734, }, # TODO: make this work for HPX # 'HpxGeom' object has no attribute 'separation' diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index e8171371a3..92b0ebf059 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -39,7 +39,7 @@ def geom(ebounds): "counts": 34366, "exposure": 3.99815e11, "exposure_image": 7.921993e10, - "background": 27984.78, + "background": 27981.977, }, { # Test single energy bin @@ -48,7 +48,7 @@ def geom(ebounds): "counts": 34366, "exposure": 1.16866e11, "exposure_image": 1.16866e11, - "background": 30418.951, + "background": 30415.916, }, { # Test single energy bin with exclusion mask @@ -58,7 +58,7 @@ def geom(ebounds): "counts": 34366, "exposure": 1.16866e11, "exposure_image": 1.16866e11, - "background": 30418.951, + "background": 30415.916, }, { # Test for different e_true and e_reco bins @@ -67,7 +67,7 @@ def geom(ebounds): "counts": 34366, "exposure": 5.971096e11, "exposure_image": 6.492968e10, - "background": 27984.78, + "background": 27981.977, }, ], ) diff --git a/gammapy/cube/tests/test_models.py b/gammapy/cube/tests/test_models.py index 5e9db5b908..2ace0a9a8e 100644 --- a/gammapy/cube/tests/test_models.py +++ b/gammapy/cube/tests/test_models.py @@ -268,7 +268,7 @@ def test_compute_flux(diffuse_evaluator): assert out.shape == (3, 4, 5) out = out.to("cm-2 s-1") assert_allclose(out.value.sum(), 633263.444803, rtol=1e-5) - assert_allclose(out.value[0, 0, 0], 1164.578565, rtol=1e-5) + assert_allclose(out.value[0, 0, 0], 1164.656176, rtol=1e-5) @staticmethod def test_apply_psf(diffuse_evaluator): @@ -277,7 +277,7 @@ def test_apply_psf(diffuse_evaluator): out = diffuse_evaluator.apply_psf(npred) assert out.data.shape == (3, 4, 5) assert_allclose(out.data.sum(), 1.106404e12, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 5.586252e08, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 5.586508e08, rtol=1e-5) @staticmethod def test_apply_edisp(diffuse_evaluator): @@ -286,14 +286,14 @@ def test_apply_edisp(diffuse_evaluator): out = diffuse_evaluator.apply_edisp(npred) assert out.data.shape == (2, 4, 5) assert_allclose(out.data.sum(), 1.606345e12, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 1.164579e09, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 1.164656e09, rtol=1e-5) @staticmethod def test_compute_npred(diffuse_evaluator): out = diffuse_evaluator.compute_npred() assert out.data.shape == (2, 4, 5) assert_allclose(out.data.sum(), 1.106403e12, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 5.586252e08, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 5.586508e08, rtol=1e-5) class TestSkyModelMapEvaluator: @@ -350,7 +350,7 @@ def test_compute_flux(evaluator): out = evaluator.compute_flux().to_value("cm-2 s-1") assert out.shape == (3, 4, 5) assert_allclose(out.sum(), 1.291427605881036e-12, rtol=1e-5) - assert_allclose(out[0, 0, 0], 4.630536473119458e-14, rtol=1e-5) + assert_allclose(out[0, 0, 0], 4.630845e-14, rtol=1e-5) @staticmethod def test_apply_psf(evaluator): @@ -359,7 +359,7 @@ def test_apply_psf(evaluator): out = evaluator.apply_psf(npred) assert out.data.shape == (3, 4, 5) assert_allclose(out.data.sum(), 2.2530737e-06, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 2.4071529e-08, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 2.407252e-08, rtol=1e-5) @staticmethod def test_apply_edisp(evaluator): @@ -368,11 +368,11 @@ def test_apply_edisp(evaluator): out = evaluator.apply_edisp(npred) assert out.data.shape == (2, 4, 5) assert_allclose(out.data.sum(), 3.2758538225058153e-06, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 4.630536473119458e-08, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 4.630845e-08, rtol=1e-5) @staticmethod def test_compute_npred(evaluator): out = evaluator.compute_npred() assert out.data.shape == (2, 4, 5) assert_allclose(out.data.sum(), 2.253073467739508e-06, rtol=1e-5) - assert_allclose(out.data[0, 0, 0], 2.4071528770264194e-08, rtol=1e-5) + assert_allclose(out.data[0, 0, 0], 2.407252e-08, rtol=1e-5) diff --git a/gammapy/maps/tests/test_wcs.py b/gammapy/maps/tests/test_wcs.py index a557643f10..b16f3e1340 100644 --- a/gammapy/maps/tests/test_wcs.py +++ b/gammapy/maps/tests/test_wcs.py @@ -167,6 +167,23 @@ def test_wcsgeom_solid_angle(): assert_allclose(solid_angle.value[0, 9, 5], 0.0003038, rtol=1e-3) +def test_wcsgeom_solid_angle_symmetry(): + geom = WcsGeom.create( + skydir=(0, 0), + coordsys="GAL", + npix=(3, 3), + binsz=20.0 * u.deg, + ) + + sa = geom.solid_angle() + + assert_allclose(sa[1, :], sa[1, 0]) # Constant along lon + assert_allclose(sa[0, 1], sa[2, 1]) # Symmetric along lat + with pytest.raises(AssertionError): + # Not constant along lat due to changes in solid angle (great circle) + assert_allclose(sa[:, 1], sa[0, 1]) + + def test_wcsgeom_solid_angle_ait(): # Pixels that don't correspond to locations on ths sky # should have solid angles set to NaN diff --git a/gammapy/maps/wcs.py b/gammapy/maps/wcs.py index d60aa71a6b..277b686911 100644 --- a/gammapy/maps/wcs.py +++ b/gammapy/maps/wcs.py @@ -773,15 +773,23 @@ def solid_angle(self): lon = coord.lon * np.pi / 180.0 lat = coord.lat * np.pi / 180.0 - # Compute solid angle using the approximation that it's - # the product between angular separation of pixel corners. + # Compute solid angle across centres of the pixels, approximating it + # as a rectangle # First index is "y", second index is "x" - ylo_xlo = lon[..., :-1, :-1], lat[..., :-1, :-1] - ylo_xhi = lon[..., :-1, 1:], lat[..., :-1, 1:] - yhi_xlo = lon[..., 1:, :-1], lat[..., 1:, :-1] - - dx = angular_separation(*(ylo_xlo + ylo_xhi)) - dy = angular_separation(*(ylo_xlo + yhi_xlo)) + # TODO: Calculate actual solid angle between two great circles? Here are two references + # suggesting more precise methods: + # https://mail.python.org/pipermail/astropy/2013-December/002632.html + # https://cta-redmine.irap.omp.eu/issues/1017 + lon_centres = (lon[..., :-1, :-1] + lon[..., 1:, 1:]) / 2 + lat_centres = (lat[..., :-1, :-1] + lat[..., 1:, 1:]) / 2 + + ymid_xlo = lon[..., :-1, :-1], lat_centres + ymid_xhi = lon[..., :-1, 1:], lat_centres + ylo_xmid = lon_centres, lat[..., :-1, 1:] + yhi_xmid = lon_centres, lat[..., 1:, :-1] + + dx = angular_separation(*(ymid_xlo + ymid_xhi)) + dy = angular_separation(*(ylo_xmid + yhi_xmid)) return u.Quantity(dx * dy, "sr", copy=False)