Merge pull request #139 from vandalt/niriss-pa

Add PA correction for NIRISS

amical/mf_pipeline/bispect.py

@@ -1185,6 +1185,10 @@ def extract_bs(
 
     with fits.open(filename) as hdu:
         hdr = hdu[0].header
+        try:
+            sci_hdr = hdu["SCI"].header
+        except KeyError:
+            sci_hdr = None
 
     infos = _check_input_infos(
         hdr, targetname=targetname, filtname=filtname, instrum=instrum, verbose=False
@@ -1383,7 +1387,7 @@ def extract_bs(
 
     # 13. Compute the absolute oriention (North-up, East-left)
     # ------------------------------------------------------------------------
-    pa = compute_pa(hdr, n_ps, display=display, verbose=verbose)
+    pa = compute_pa(hdr, n_ps, display=display, verbose=verbose, sci_hdr=sci_hdr)
 
     # Compile informations in the storage infos class
     infos = _add_infos_header(infos, hdr, mf, pa, filename, maskname, npix)

amical/tests/test_cli.py

@@ -12,6 +12,11 @@
 
 valid_commands = ["clean", "extract", "calibrate"]
 
+# The test data has no SCI header and does not need one: ignore warning
+pytestmark = pytest.mark.filterwarnings(
+    "ignore: No SCI header for NIRISS. No PA correction will be applied."
+)
+
 
 @pytest.fixture()
 def close_figures():

amical/tests/test_io.py

@@ -14,6 +14,11 @@
 # from amical.mf_pipeline.ami_function import find_bad_BL_BS
 # from amical.mf_pipeline.ami_function import find_bad_holes
 
+# The test data has no SCI header and does not need one: ignore warning
+pytestmark = pytest.mark.filterwarnings(
+    "ignore: No SCI header for NIRISS. No PA correction will be applied."
+)
+
 
 @pytest.fixture()
 def close_figures():

amical/tests/test_tools.py

@@ -51,6 +51,56 @@ def test_SPHERE_parang(global_datadir):
     with fits.open(global_datadir / "hdr_sphere.fits") as hdu:
         hdr = hdu[0].header
     n_ps = 1
-    pa = tools.compute_pa(hdr, n_ps=n_ps)
+    pa = tools.sphere_parang(hdr, n_dit_ifs=n_ps)
     true_pa = 109  # Human value
     assert pa == pytest.approx(true_pa, 0.01)
+
+
+def test_NIRISS_parang(global_datadir):
+    with fits.open(global_datadir / "hdr_niriss_mirage.fits") as hdu:
+        hdr = hdu["SCI"].header
+    pa = tools.niriss_parang(hdr)
+    true_pa = 157.9079  # Human value
+    assert pa == pytest.approx(true_pa, 0.01)
+
+
+def test_compute_pa_sphere(global_datadir):
+    with fits.open(global_datadir / "hdr_sphere.fits") as hdul:
+        hdr = hdul[0].header
+    n_ps = 1
+    pa = tools.compute_pa(hdr, n_ps)
+    true_pa = 109  # Human value
+    assert pa == pytest.approx(true_pa, 0.01)
+
+
+def test_compute_pa_niriss(global_datadir):
+    with fits.open(global_datadir / "hdr_niriss_mirage.fits") as hdul:
+        hdr = hdul[0].header
+        sci_hdr = hdul["SCI"].header
+        n_ps = hdul["SCI"].data.shape[-1]
+    pa = tools.compute_pa(hdr, n_ps, sci_hdr=sci_hdr)
+    true_pa = 157.9079  # Human value
+    assert pa == pytest.approx(true_pa, 0.01)
+
+
+def test_NIRISS_parang_amisim():
+    # ami_sim file has no SCI
+    hdr = None
+    with pytest.warns(RuntimeWarning) as record:
+        pa = tools.niriss_parang(hdr)
+    assert len(record) == 1
+    assert (
+        record[0].message.args[0]
+        == "No SCI header for NIRISS. No PA correction will be applied."
+    )
+    assert pa == 0.0
+
+
+def test_compute_pa_niriss_amisim(global_datadir):
+    with fits.open(global_datadir / "hdr_niriss_amisim.fits") as hdul:
+        hdr = hdul[0].header
+        n_ps = hdul[0].data.shape[-1]
+    with pytest.warns(RuntimeWarning) as record:
+        pa = tools.compute_pa(hdr, n_ps)
+    assert len(record) == 1
+    assert pa == 0.0

amical/tools.py

@@ -10,6 +10,7 @@
 --------------------------------------------------------------------
 """
 import math as m
+import warnings
 
 import numpy as np
 from termcolor import cprint
@@ -366,10 +367,11 @@ def jd2lst(lng, jd):
     return lst
 
 
-def compute_pa(hdr, n_ps, verbose=False, display=False):
+def compute_pa(hdr, n_ps, verbose=False, display=False, *, sci_hdr=None):
 
     list_fct_pa = {
-        "SPHERE": sphere_parang,
+        "SPHERE": (sphere_parang, {"hdr": hdr, "n_dit_ifs": n_ps}),
+        "NIRISS": (niriss_parang, {"hdr": sci_hdr}),
     }
 
     instrument = hdr["INSTRUME"]
@@ -387,7 +389,8 @@ def compute_pa(hdr, n_ps, verbose=False, display=False):
         pa_exist = False
         l_pa = np.zeros(nframe)
     else:
-        l_pa = list_fct_pa[instrument](hdr, n_ps)
+        fct_pa, kwargs_pa = list_fct_pa[instrument]
+        l_pa = fct_pa(**kwargs_pa)
         pa_exist = True
 
     pa = np.mean(l_pa)
@@ -407,6 +410,19 @@ def compute_pa(hdr, n_ps, verbose=False, display=False):
     return pa
 
 
+def niriss_parang(hdr):
+    if hdr is None:
+        warnings.warn(
+            "No SCI header for NIRISS. No PA correction will be applied.",
+            RuntimeWarning,
+        )
+        return 0.0
+    v3i_yang = hdr["V3I_YANG"]  # Angle from V3 axis to ideal y axis (deg)
+    roll_ref_pa = hdr["ROLL_REF"]  # Offset between V3 and N in local aperture coord
+
+    return roll_ref_pa - v3i_yang
+
+
 def sphere_parang(hdr, n_dit_ifs=None):
     """
     Reads the header and creates an array giving the paralactic angle for each frame,