Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix symmetry issue in solid angle calculation for WcsGeom #2035

Merged
merged 7 commits into from Feb 18, 2019
Merged

Fix symmetry issue in solid angle calculation for WcsGeom #2035

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
43b3f86
Create test that requires symmetric solid angle calculation
watsonjj Feb 15, 2019
14893e0
Calculate separations through the centre of the pixel
watsonjj Feb 15, 2019
aad5f23
Correct fixed values in tests (change on the order of 0.01%)
watsonjj Feb 15, 2019
06c2003
Add references for solid angle computation improvements
adonath Feb 18, 2019
74b9212
Merge branch 'master' into symmetric_solid_angle
watsonjj Feb 18, 2019
6567659
Merge branch 'master' into symmetric_solid_angle
watsonjj Feb 18, 2019
f502df7
Merge remote-tracking branch 'origin/symmetric_solid_angle' into symm…
watsonjj Feb 18, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 2 additions & 2 deletions gammapy/cube/tests/test_background.py
View file
Open in desktop
Expand Up @@ -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'
Expand Down
8 changes: 4 additions & 4 deletions gammapy/cube/tests/test_make.py
View file
Open in desktop
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -67,7 +67,7 @@ def geom(ebounds):
"counts": 34366,
"exposure": 5.971096e11,
"exposure_image": 6.492968e10,
"background": 27984.78,
"background": 27981.977,
},
],
)
Expand Down
16 changes: 8 additions & 8 deletions gammapy/cube/tests/test_models.py
View file
Open in desktop
Expand Up @@ -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):
Expand All @@ -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):
Expand All @@ -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:
Expand Down Expand Up @@ -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):
Expand All @@ -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):
Expand All @@ -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)
17 changes: 17 additions & 0 deletions gammapy/maps/tests/test_wcs.py
View file
Open in desktop
Expand Up @@ -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
Expand Down
24 changes: 16 additions & 8 deletions gammapy/maps/wcs.py
View file
Open in desktop
Expand Up @@ -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)

Expand Down