oifits

Read, write, and display OIFITS interferometric observables.

cal2dict(cal, target=None, pa=0, del_pa=0, snr=4, true_flag_v2=True, true_flag_t3=False, oriented=-1, ind_hole=None)

Format calibrated data as an OIFITS-compatible dictionary.

Parameters:

Name	Type	Description	Default
cal	object	Calibrated data object.	required
target	str or None	Science-target name; uses the target metadata when omitted.	None
pa	float	Position angle and its uncertainty in degrees.	0
del_pa	float	Position angle and its uncertainty in degrees.	0
snr	float	Signal-to-noise threshold for flags.	4
true_flag_v2	bool	Whether to compute squared-visibility and closure-phase flags.	True, False
true_flag_t3	bool	Whether to compute squared-visibility and closure-phase flags.	True, False
oriented	float	East orientation; negative values put east on the left.	-1
ind_hole	int or None	Hole used to select independent closure phases.	None

Returns:

Type	Description
dict	OIFITS-compatible dictionary.

Raises:

Type	Description
KeyError	If required metadata is missing from the observables.

Source code in src/amical/oifits.py

def cal2dict(
    cal,
    target=None,
    pa=0,
    del_pa=0,
    snr=4,
    true_flag_v2=True,
    true_flag_t3=False,
    oriented=-1,
    ind_hole=None,
):
    """Format calibrated data as an OIFITS-compatible dictionary.

    Parameters
    ----------
    cal : object
        Calibrated data object.
    target : str or None, default=None
        Science-target name; uses the target metadata when omitted.
    pa, del_pa : float, default=0
        Position angle and its uncertainty in degrees.
    snr : float, default=4
        Signal-to-noise threshold for flags.
    true_flag_v2, true_flag_t3 : bool, default=True, False
        Whether to compute squared-visibility and closure-phase flags.
    oriented : float, default=-1
        East orientation; negative values put east on the left.
    ind_hole : int or None, default=None
        Hole used to select independent closure phases.

    Returns
    -------
    dict
        OIFITS-compatible dictionary.

    Raises
    ------
    KeyError
        If required metadata is missing from the observables.
    """

    res_t = cal.raw_t
    res_c = cal.raw_c
    n_baselines = res_t.mask.n_baselines
    n_bs = len(cal.cp)
    bl2h_ix = res_t.mask.bl2h_ix
    bs2bl_ix = res_t.mask.bs2bl_ix

    try:
        ins = res_t.infos.instrument
        maskname = res_t.infos.maskname
        pixscale = res_t.infos.pixscale
    except (KeyError, AttributeError) as e:
        raise KeyError(
            "'INSTRUME', 'NRMNAME' or 'PIXELSCL' are not in the header."
        ) from e

    t = _observation_time(res_t.infos)
    exp_time = _exposure_time(res_t.infos.get("hdr") or {})

    if true_flag_v2:
        flagV2 = _compute_flag(cal.vis2, cal.e_vis2, limit=snr)
    else:
        flagV2 = np.array([False] * n_baselines)

    if true_flag_t3:
        flagCP = _compute_flag(cal.cp, cal.e_cp, limit=snr)
    else:
        flagCP = np.array([False] * n_bs)

    sta_index_v2 = []
    for i in range(n_baselines):
        sta_index_v2.append(np.array(bl2h_ix[:, i]))
    sta_index_v2 = np.array(sta_index_v2)

    # Triangle (i, j, k) in the order of the closure phase: baselines ij
    # (U1/V1) and jk (U2/V2) close with ki.
    sta_index_t3 = np.stack(
        [
            bl2h_ix[0, bs2bl_ix[0, :]],
            bl2h_ix[1, bs2bl_ix[0, :]],
            bl2h_ix[1, bs2bl_ix[1, :]],
        ],
        axis=1,
    )

    if target is None:
        target = res_t.infos.target

    thepa = pa - 0.5 * del_pa
    u, v = oriented * res_t.u, res_t.v
    u1 = u * np.cos(np.deg2rad(thepa)) + v * np.sin(np.deg2rad(thepa))
    v1 = -u * np.sin(np.deg2rad(thepa)) + v * np.cos(np.deg2rad(thepa))

    if isinstance(res_c, list):
        calib_name = res_c[0].infos.target
    else:
        calib_name = res_c.infos.target

    if ind_hole is not None:
        rprint(f"[green]Select only independant CP using common hole #{ind_hole}.")
        sel_ind_cp = _peak1hole_cp(cal.raw_t, ind_hole)
    else:
        sel_ind_cp = np.arange(len(cal.cp))

    dic = {
        "OI_VIS2": {
            "VIS2DATA": cal.vis2,
            "VIS2ERR": cal.e_vis2,
            "UCOORD": u1,
            "VCOORD": v1,
            "STA_INDEX": sta_index_v2,
            "MJD": t.mjd,
            "INT_TIME": exp_time,
            "TIME": 0,
            "TARGET_ID": 1,
            "FLAG": flagV2,
            "BL": (u1**2 + v1**2) ** 0.5,
        },
        "OI_T3": {
            "MJD": t.mjd,
            "INT_TIME": exp_time,
            "T3PHI": cal.cp[sel_ind_cp],
            "T3PHIERR": cal.e_cp[sel_ind_cp],
            "TIME": 0,
            "T3AMP": np.zeros(len(cal.cp))[sel_ind_cp],
            "T3AMPERR": np.zeros(len(cal.cp))[sel_ind_cp],
            "U1COORD": u1[bs2bl_ix[0, :]][sel_ind_cp],
            "V1COORD": v1[bs2bl_ix[0, :]][sel_ind_cp],
            "U2COORD": u1[bs2bl_ix[1, :]][sel_ind_cp],
            "V2COORD": v1[bs2bl_ix[1, :]][sel_ind_cp],
            "STA_INDEX": list(sta_index_t3[sel_ind_cp]),
            "FLAG": flagCP[sel_ind_cp],
            "BL": res_t.bl_cp[sel_ind_cp],
        },
        "OI_WAVELENGTH": {
            "EFF_WAVE": np.array([res_t.wl]),
            "EFF_BAND": np.array([res_t.e_wl]),
        },
        "info": {
            "TARGET": target,
            "CALIB": calib_name,
            "FILT": res_t.infos.filtname,
            "INSTRUME": ins,
            "MASK": maskname,
            "MJD": t.mjd,
            "HDR": res_t.infos,
            "ISZ": res_t.infos.isz,
            "PSCALE": pixscale,
            "xycoord": res_t.mask.xycoord,
            "SEEING": res_t.infos.seeing,
        },
    }
    return dic

load(filename, filtname=None)

Load an OIFITS file into a dictionary.

Parameters:

Name	Type	Description	Default
filename	str or path - like	OIFITS filename.	required
filtname	str or None	Filter name used when the header does not provide one.	None

Returns:

Type	Description
dict	OIFITS data and header information.

Source code in src/amical/oifits.py

def load(filename, filtname=None):
    """Load an OIFITS file into a dictionary.

    Parameters
    ----------
    filename : str or path-like
        OIFITS filename.
    filtname : str or None, default=None
        Filter name used when the header does not provide one.

    Returns
    -------
    dict
        OIFITS data and header information.
    """
    from astropy.io import fits

    fitsHandler = fits.open(filename)
    hdr = fitsHandler[0].header

    dic = {}
    for hdu in fitsHandler[1:]:
        if hdu.header["EXTNAME"] == "OI_WAVELENGTH":
            dic["OI_WAVELENGTH"] = {
                "EFF_WAVE": hdu.data["EFF_WAVE"],
                "EFF_BAND": hdu.data["EFF_BAND"],
            }

        if hdu.header["EXTNAME"] == "OI_VIS2":
            dic["OI_VIS2"] = {
                "VIS2DATA": hdu.data["VIS2DATA"],
                "VIS2ERR": hdu.data["VIS2ERR"],
                "UCOORD": hdu.data["UCOORD"],
                "VCOORD": hdu.data["VCOORD"],
                "STA_INDEX": hdu.data["STA_INDEX"],
                "MJD": hdu.data["MJD"],
                "INT_TIME": hdu.data["INT_TIME"],
                "TIME": hdu.data["TIME"],
                "TARGET_ID": hdu.data["TARGET_ID"],
                "FLAG": np.array(hdu.data["FLAG"]),
            }
            try:
                dic["OI_VIS2"]["BL"] = hdu.data["BL"]
            except KeyError:
                dic["OI_VIS2"]["BL"] = (
                    hdu.data["UCOORD"] ** 2 + hdu.data["VCOORD"] ** 2
                ) ** 0.5

        if hdu.header["EXTNAME"] == "OI_VIS":
            dic["OI_VIS"] = {
                "TARGET_ID": hdu.data["TARGET_ID"],
                "TIME": hdu.data["TIME"],
                "MJD": hdu.data["MJD"],
                "INT_TIME": hdu.data["INT_TIME"],
                "VISAMP": hdu.data["VISAMP"],
                "VISAMPERR": hdu.data["VISAMPERR"],
                "VISPHI": hdu.data["VISPHI"],
                "VISPHIERR": hdu.data["VISPHIERR"],
                "UCOORD": hdu.data["UCOORD"],
                "VCOORD": hdu.data["VCOORD"],
                "STA_INDEX": hdu.data["STA_INDEX"],
                "FLAG": hdu.data["FLAG"],
            }
            try:
                dic["OI_VIS"]["BL"] = hdu.data["BL"]
            except KeyError:
                dic["OI_VIS"]["BL"] = (
                    hdu.data["UCOORD"] ** 2 + hdu.data["VCOORD"] ** 2
                ) ** 0.5

        if hdu.header["EXTNAME"] == "OI_T3":
            u1 = hdu.data["U1COORD"]
            u2 = hdu.data["U2COORD"]
            v1 = hdu.data["V1COORD"]
            v2 = hdu.data["V2COORD"]
            u3 = -(u1 + u2)
            v3 = -(v1 + v2)
            bl_cp = []
            for k in range(len(u1)):
                B1 = np.sqrt(u1[k] ** 2 + v1[k] ** 2)
                B2 = np.sqrt(u2[k] ** 2 + v2[k] ** 2)
                B3 = np.sqrt(u3[k] ** 2 + v3[k] ** 2)
                bl_cp.append(np.max([B1, B2, B3]))  # rad-1
            bl_cp = np.array(bl_cp)

            dic["OI_T3"] = {
                "T3PHI": hdu.data["T3PHI"],
                "T3PHIERR": hdu.data["T3PHIERR"],
                "T3AMP": hdu.data["T3AMP"],
                "T3AMPERR": hdu.data["T3AMPERR"],
                "U1COORD": hdu.data["U1COORD"],
                "V1COORD": hdu.data["V1COORD"],
                "U2COORD": hdu.data["U2COORD"],
                "V2COORD": hdu.data["V2COORD"],
                "STA_INDEX": hdu.data["STA_INDEX"],
                "MJD": hdu.data["MJD"],
                "FLAG": hdu.data["FLAG"],
                "TARGET_ID": hdu.data["TARGET_ID"],
                "TIME": hdu.data["TIME"],
                "INT_TIME": hdu.data["INT_TIME"],
            }
            try:
                dic["OI_T3"]["BL"] = hdu.data["FREQ"]
            except KeyError:
                dic["OI_T3"]["BL"] = bl_cp

    fitsHandler.close()

    dic["info"] = {h: hdr[h] for h in hdr}  # {'MJD': mjd,

    if "FILT" not in list(dic["info"].keys()):
        if filtname is None:
            print("No filter in info or as input param?")
        else:
            dic["info"]["FILT"] = filtname
    return dic

loadc(filename)

Load an OIFITS file as an attribute-accessible object.

Parameters:

Name	Type	Description	Default
filename	str or path - like	OIFITS filename.	required

Returns:

Type	Description
Munch	Loaded observables and metadata.

Source code in src/amical/oifits.py

def loadc(filename):
    """Load an OIFITS file as an attribute-accessible object.

    Parameters
    ----------
    filename : str or path-like
        OIFITS filename.

    Returns
    -------
    munch.Munch
        Loaded observables and metadata.
    """

    dic = load(filename)
    res = {}
    # Extract infos
    res["target"] = dic["info"].get("OBJECT")

    res["calib"] = dic["info"].get("CALIB")
    res["seeing"] = dic["info"].get("SEEING")
    res["mjd"] = dic["info"].get("MJD")
    if res["mjd"] is None:
        res["mjd"] = dic["info"].get("MJD-OBS")

    # Extract wavelength
    res["wl"] = dic["OI_WAVELENGTH"]["EFF_WAVE"]
    res["e_wl"] = dic["OI_WAVELENGTH"]["EFF_BAND"]

    # Extract squared visibilities
    res["vis2"] = dic["OI_VIS2"]["VIS2DATA"]
    res["e_vis2"] = dic["OI_VIS2"]["VIS2ERR"]
    res["u"] = dic["OI_VIS2"]["UCOORD"]
    res["v"] = dic["OI_VIS2"]["VCOORD"]
    res["bl"] = dic["OI_VIS2"]["BL"]
    res["flag_vis"] = dic["OI_VIS2"]["FLAG"]

    # Extract closure phases
    res["cp"] = dic["OI_T3"]["T3PHI"]
    res["e_cp"] = dic["OI_T3"]["T3PHIERR"]
    res["u1"] = dic["OI_T3"]["U1COORD"]
    res["v1"] = dic["OI_T3"]["V1COORD"]
    res["u2"] = dic["OI_T3"]["U2COORD"]
    res["v2"] = dic["OI_T3"]["V2COORD"]
    res["bl_cp"] = dic["OI_T3"]["BL"]
    res["flag_cp"] = dic["OI_T3"]["FLAG"]

    return dict2class(res)

save(observables, oifits_file=None, datadir='Saveoifits', pa=0, ind_hole=None, fake_obj=False, true_flag_v2=True, true_flag_t3=False, snr=4, verbose=False, *, origin=None, raw=False, query_simbad=None)

Save calibrated observables in OIFITS format.

Parameters:

Name	Type	Description	Default
observables	object or list of object	Calibrated observables, including IFU observables when given as a list.	required
oifits_file	str or path - like	OIFITS filename.	None
datadir	str or path - like	Output directory.	"Saveoifits"
pa	float	Position angle in degrees.	0
ind_hole	int or None	If set, save only the (N-1)(N-2)/2 independent closure phases that include this hole; otherwise save all N(N-1)(N-2)/6.	None
fake_obj	bool	Whether observables originate from simulated data. SIMBAD is then not queried and the coordinates are set to zero.	False
true_flag_v2	bool	Whether to compute squared-visibility and closure-phase flags.	True, False
true_flag_t3	bool	Whether to compute squared-visibility and closure-phase flags.	True, False
snr	float	Signal-to-noise threshold for flags.	4
verbose	bool	Whether to print status information.	False
origin	(str or None, keyword - only)	Value for the OIFITS ORIGIN key.	None
raw	(bool, keyword - only)	Whether the input is uncalibrated.	False
query_simbad	(bool or None, keyword - only)	Whether to query SIMBAD for the target coordinates, proper motion, parallax and spectral type. If False (e.g. on compute nodes without internet access), the RA/DEC keywords of the original header are used instead. If None, SIMBAD is queried unless the environment variable AMICAL_NO_SIMBAD is set to a true value (e.g. AMICAL_NO_SIMBAD=1).	None

Returns:

Name	Type	Description
dictionaries	list of dict	OIFITS-formatted dictionaries.
savedfile	str	Saved OIFITS filename.

Source code in src/amical/oifits.py

def save(
    observables,
    oifits_file=None,
    datadir="Saveoifits",
    pa=0,
    ind_hole=None,
    fake_obj=False,
    true_flag_v2=True,
    true_flag_t3=False,
    snr=4,
    verbose=False,
    *,
    origin=None,
    raw=False,
    query_simbad=None,
):
    """Save calibrated observables in OIFITS format.

    Parameters
    ----------
    observables : object or list of object
        Calibrated observables, including IFU observables when given as a list.
    oifits_file : str or path-like, optional
        OIFITS filename.
    datadir : str or path-like, default="Saveoifits"
        Output directory.
    pa : float, default=0
        Position angle in degrees.
    ind_hole : int or None, default=None
        If set, save only the (N-1)(N-2)/2 independent closure phases that
        include this hole; otherwise save all N(N-1)(N-2)/6.
    fake_obj : bool, default=False
        Whether observables originate from simulated data. SIMBAD is then not
        queried and the coordinates are set to zero.
    true_flag_v2, true_flag_t3 : bool, default=True, False
        Whether to compute squared-visibility and closure-phase flags.
    snr : float, default=4
        Signal-to-noise threshold for flags.
    verbose : bool, default=False
        Whether to print status information.
    origin : str or None, keyword-only, default=None
        Value for the OIFITS ORIGIN key.
    raw : bool, keyword-only, default=False
        Whether the input is uncalibrated.
    query_simbad : bool or None, keyword-only, default=None
        Whether to query SIMBAD for the target coordinates, proper motion,
        parallax and spectral type. If False (e.g. on compute nodes without
        internet access), the RA/DEC keywords of the original header are used
        instead. If None, SIMBAD is queried unless the environment variable
        ``AMICAL_NO_SIMBAD`` is set to a true value (e.g. ``AMICAL_NO_SIMBAD=1``).

    Returns
    -------
    dictionaries : list of dict
        OIFITS-formatted dictionaries.
    savedfile : str
        Saved OIFITS filename.
    """
    from astropy.io import fits

    if observables is None:
        rprint("[on red]\nError save : Wrong data format!", file=sys.stderr)
        return None

    if oifits_file is None:
        print("Error: oifits filename is not given, please specify oifits_file.")
        return None

    if not isinstance(observables, list):
        observables = [observables]

    if not isinstance(origin, str | type(None)):
        raise TypeError("origin should be a str or None")

    l_dic = []
    for iobs in observables:
        # If keys correspond to raw observables, make format compatible with calibrated
        if "cp" in iobs and "raw_t" not in iobs:
            if not raw and verbose:
                msg = (
                    "The input seems to contain uncalibrated observables."
                    " Saving raw observables as oifits is provided only for"
                    " convenience. To re-use the observables in amical.calibrate(),"
                    " they should be saved to a pickle file. Use raw=True to turn this"
                    " warning off."
                )
                rprint(f"[green]Warning: {msg}", file=sys.stderr)
            iobs = wrap_raw(iobs)
        idic = cal2dict(
            iobs,
            pa=pa,
            true_flag_v2=true_flag_v2,
            true_flag_t3=true_flag_t3,
            ind_hole=ind_hole,
            snr=snr,
        )
        l_dic.append(idic)
    dic = l_dic[0]

    if not os.path.exists(datadir):
        print(f"### Create {datadir} directory to save all requested Oifits ###")
        os.mkdir(datadir)

    # ------------------------------
    #       Creation OIFITS
    # ------------------------------
    if verbose:
        print(f"\n\n### Init creation of OI_FITS ({oifits_file}) :")

    # refdate = datetime.datetime(2000, 1, 1)  # Unix time reference
    hdulist = fits.HDUList()

    hdr = dic["info"].get("HDR", {})
    hdu = fits.PrimaryHDU()
    hdu.header["DATE"] = datetime.datetime.now().strftime(
        format="%F"
    )  # , 'Creation date'
    hdu.header["ORIGIN"] = origin or hdr.get("ORIGIN", "Sydney University")
    hdu.header["CONTENT"] = "OIFITS2"
    hdu.header["DATE-OBS"] = hdr.get("date-obs", "")
    hdu.header["TELESCOP"] = hdr.get("telescop", "")
    hdu.header["INSTRUME"] = hdr.get("INSTRUME", dic["info"]["INSTRUME"])
    hdu.header["OBSERVER"] = hdr.get("OBSERVER", "")
    hdu.header["OBJECT"] = hdr.get("OBJECT", dic["info"]["TARGET"])
    hdu.header["INSMODE"] = "NRM"
    hdu.header["FILT"] = dic["info"]["FILT"]
    hdu.header["MJD"] = dic["info"]["MJD"]
    hdu.header["MASK"] = dic["info"]["MASK"]
    hdu.header["SEEING"] = dic["info"].get("SEEING", "")
    hdu.header["CALIB"] = dic["info"]["CALIB"]

    hdulist.append(hdu)
    # ------------------------------
    #        OI Wavelength
    # ------------------------------

    if verbose:
        print("-> Including OI Wavelength table...")
    data = dic["OI_WAVELENGTH"]

    EFF_WAVE, EFF_BAND = [], []
    for dic in l_dic:
        d = dic["OI_WAVELENGTH"]
        EFF_WAVE.append(d["EFF_WAVE"][0])
        EFF_BAND.append(d["EFF_BAND"][0])

    iwl = len(EFF_WAVE)
    # Data
    # -> Initiation new hdu table :
    hdu = fits.BinTableHDU.from_columns(
        fits.ColDefs(
            (
                fits.Column(
                    name="EFF_WAVE", format="1E", unit="METERS", array=EFF_WAVE
                ),
                fits.Column(
                    name="EFF_BAND", format="1E", unit="METERS", array=EFF_BAND
                ),
            )
        )
    )

    # Header
    hdu.header["EXTNAME"] = "OI_WAVELENGTH"
    hdu.header["OI_REVN"] = 2  # , 'Revision number of the table definition'
    # 'Name of detector, for cross-referencing'
    hdu.header["INSNAME"] = dic["info"]["INSTRUME"]
    hdulist.append(hdu)  # Add current HDU to the final fits file.

    # ------------------------------
    #          OI Target
    # ------------------------------
    if verbose:
        print("-> Including OI Target table...")

    name_star = dic["info"]["TARGET"]
    target = _target_info(
        name_star, hdr.get("hdr"), fake_obj=fake_obj, query_simbad=query_simbad
    )
    ra, dec = [target["ra"]], [target["dec"]]
    pmra, pmdec, plx = [target["pmra"]], [target["pmdec"]], [target["plx"]]
    spectyp = [target["spectyp"]]

    if name_star == "":
        name_star = "Unknown"

    hdu = fits.BinTableHDU.from_columns(
        fits.ColDefs(
            (
                fits.Column(name="TARGET_ID", format="1I", array=[1]),
                fits.Column(name="TARGET", format="16A", array=[name_star]),
                fits.Column(name="RAEP0", format="1D", unit="DEGREES", array=ra),
                fits.Column(name="DECEP0", format="1D", unit="DEGREES", array=dec),
                fits.Column(name="EQUINOX", format="1E", unit="YEARS", array=[2000]),
                fits.Column(name="RA_ERR", format="1D", unit="DEGREES", array=[0]),
                fits.Column(name="DEC_ERR", format="1D", unit="DEGREES", array=[0]),
                fits.Column(name="SYSVEL", format="1D", unit="M/S", array=[0]),
                fits.Column(name="VELTYP", format="8A", array=["UNKNOWN"]),
                fits.Column(name="VELDEF", format="8A", array=["OPTICAL"]),
                fits.Column(name="PMRA", format="1D", unit="DEG/YR", array=pmra),
                fits.Column(name="PMDEC", format="1D", unit="DEG/YR", array=pmdec),
                fits.Column(name="PMRA_ERR", format="1D", unit="DEG/YR", array=[0]),
                fits.Column(name="PMDEC_ERR", format="1D", unit="DEG/YR", array=[0]),
                fits.Column(name="PARALLAX", format="1E", unit="DEGREES", array=plx),
                fits.Column(name="PARA_ERR", format="1E", unit="DEGREES", array=[0]),
                fits.Column(name="SPECTYP", format="16A", array=spectyp),
            )
        )
    )

    hdu.header["EXTNAME"] = "OI_TARGET"
    hdu.header["OI_REVN"] = 2, "Revision number of the table definition"
    hdulist.append(hdu)

    # ------------------------------
    #           OI Array
    # ------------------------------

    if verbose:
        print("-> Including OI Array table...")

    staxy = dic["info"]["xycoord"]
    N_ap = len(staxy)
    telName = [f"A{x}" for x in range(1, N_ap + 1)]
    staName = telName
    diameter = [0] * N_ap

    staxyz = []
    for x in staxy:
        a = list(x)
        line = [a[0], a[1], 0]
        staxyz.append(line)

    staIndex = np.arange(N_ap) + 1

    pscale = dic["info"]["PSCALE"] / 1000.0  # arcsec
    isz = dic["info"]["ISZ"]  # Size of the image to extract NRM data
    fov = [pscale * isz / 2.0] * N_ap  # radius of the extracted image
    fovtype = ["RADIUS"] * N_ap

    hdu = fits.BinTableHDU.from_columns(
        fits.ColDefs(
            (
                fits.Column(name="TEL_NAME", format="16A", array=telName),
                fits.Column(name="STA_NAME", format="16A", array=staName),
                fits.Column(name="STA_INDEX", format="1I", array=staIndex),
                fits.Column(
                    name="DIAMETER", unit="METERS", format="1E", array=diameter
                ),
                fits.Column(name="STAXYZ", unit="METERS", format="3D", array=staxyz),
                fits.Column(name="FOV", unit="ARCSEC", format="1D", array=fov),
                fits.Column(name="FOVTYPE", format="6A", array=fovtype),
            )
        )
    )

    hdu.header["EXTNAME"] = "OI_ARRAY"
    hdu.header["ARRAYX"] = float(0)
    hdu.header["ARRAYY"] = float(0)
    hdu.header["ARRAYZ"] = float(0)
    hdu.header["ARRNAME"] = dic["info"]["MASK"]
    hdu.header["FRAME"] = "SKY"
    hdu.header["OI_REVN"] = 2, "Revision number of the table definition"
    hdulist.append(hdu)

    # ------------------------------
    #           OI VIS2
    # ------------------------------
    if verbose:
        print("-> Including OI Vis2 table...")

    data = dic["OI_VIS2"]
    if not isinstance(data["TARGET_ID"], np.ndarray):
        npts = len(data["VIS2DATA"])
    else:
        npts = 1

    if isinstance(data["TARGET_ID"], int):
        npts = len(dic["OI_VIS2"]["VIS2DATA"])
        targetId = [data["TARGET_ID"]] * npts
        time = [data["TIME"]] * npts
        mjd = [data["MJD"]] * npts
        intTime = [data["INT_TIME"]] * npts
    else:
        npts = 1
        targetId = [1] * len(data["VIS2DATA"])
        time = data["TIME"]
        mjd = data["MJD"]
        intTime = data["INT_TIME"]

    staIndex = _format_staindex_v2(data["STA_INDEX"])

    VIS2DATA, VIS2ERR, FLAG = [], [], []
    for dic in l_dic:
        d = dic["OI_VIS2"]
        VIS2DATA.append(d["VIS2DATA"])
        VIS2ERR.append(d["VIS2ERR"])
        FLAG.append(d["FLAG"])
    VIS2DATA = np.array(VIS2DATA).T
    VIS2ERR = np.array(VIS2ERR).T
    FLAG = np.array(FLAG).T

    hdu = fits.BinTableHDU.from_columns(
        fits.ColDefs(
            [
                fits.Column(name="TARGET_ID", format="1I", array=targetId),
                fits.Column(name="TIME", format="1D", unit="SECONDS", array=time),
                fits.Column(name="MJD", unit="DAY", format="1D", array=mjd),
                fits.Column(
                    name="INT_TIME", format="1D", unit="SECONDS", array=intTime
                ),
                fits.Column(name="VIS2DATA", format=f"{iwl}D", array=VIS2DATA),
                fits.Column(name="VIS2ERR", format=f"{iwl}D", array=VIS2ERR),
                fits.Column(
                    name="UCOORD", format="1D", unit="METERS", array=data["UCOORD"]
                ),
                fits.Column(
                    name="VCOORD", format="1D", unit="METERS", array=data["VCOORD"]
                ),
                fits.Column(name="STA_INDEX", format="2I", array=staIndex),
                fits.Column(name="FLAG", format=f"{iwl}L", array=FLAG),
            ]
        )
    )

    hdu.header["EXTNAME"] = "OI_VIS2"
    hdu.header["INSNAME"] = dic["info"]["INSTRUME"]
    hdu.header["ARRNAME"] = dic["info"]["MASK"]
    hdu.header["OI_REVN"] = 2, "Revision number of the table definition"
    hdu.header["DATE-OBS"] = hdr.get("date-obs", "")
    hdulist.append(hdu)

    # ------------------------------
    #           OI T3
    # ------------------------------
    if verbose:
        print("-> Including OI T3 table...")

    data = dic["OI_T3"]
    try:
        check_oi = type(float(data["MJD"]))
    except TypeError:
        check_oi = int

    if check_oi is float:
        t3phi = dic["OI_T3"]["T3PHI"]
        npts = len(t3phi)
        targetId = np.ones_like(t3phi)
        time = np.zeros_like(t3phi)
        mjd = [data["MJD"]] * npts
        intTime = [data["INT_TIME"]] * npts
    else:
        npts = 1
        targetId = np.ones_like(data["T3PHI"])
        time = data["TIME"]
        mjd = data["MJD"]
        intTime = data["INT_TIME"]

    T3AMP, T3AMPERR, T3PHI, T3PHIERR, FLAG = [], [], [], [], []
    for dic in l_dic:
        d = dic["OI_T3"]
        T3AMP.append(d["T3AMP"])
        T3AMPERR.append(d["T3AMPERR"])
        T3PHI.append(d["T3PHI"])
        T3PHIERR.append(d["T3PHIERR"])
        FLAG.append(d["FLAG"])
    T3AMP = np.array(T3AMP).T
    T3AMPERR = np.array(T3AMPERR).T
    T3PHI = np.array(T3PHI).T
    T3PHIERR = np.array(T3PHIERR).T
    FLAG = np.array(FLAG).T

    staIndex = _format_staindex_t3(data["STA_INDEX"])
    hdu = fits.BinTableHDU.from_columns(
        fits.ColDefs(
            (
                fits.Column(name="TARGET_ID", format="1I", array=targetId),
                fits.Column(name="TIME", format="1D", unit="SECONDS", array=time),
                fits.Column(name="MJD", format="1D", unit="DAY", array=mjd),
                fits.Column(
                    name="INT_TIME", format="1D", unit="SECONDS", array=intTime
                ),
                fits.Column(name="T3AMP", format=f"{iwl}D", array=T3AMP),
                fits.Column(name="T3AMPERR", format=f"{iwl}D", array=T3AMPERR),
                fits.Column(
                    name="T3PHI", format=f"{iwl}D", unit="DEGREES", array=T3PHI
                ),
                fits.Column(
                    name="T3PHIERR", format=f"{iwl}D", unit="DEGREES", array=T3PHIERR
                ),
                fits.Column(
                    name="U1COORD", format="1D", unit="METERS", array=data["U1COORD"]
                ),
                fits.Column(
                    name="V1COORD", format="1D", unit="METERS", array=data["V1COORD"]
                ),
                fits.Column(
                    name="U2COORD", format="1D", unit="METERS", array=data["U2COORD"]
                ),
                fits.Column(
                    name="V2COORD", format="1D", unit="METERS", array=data["V2COORD"]
                ),
                fits.Column(name="STA_INDEX", format="3I", array=staIndex),
                fits.Column(name="FLAG", format=f"{iwl}L", array=FLAG),
            )
        )
    )

    hdu.header["EXTNAME"] = "OI_T3"
    hdu.header["INSNAME"] = dic["info"]["INSTRUME"]
    hdu.header["ARRNAME"] = dic["info"]["MASK"]
    hdu.header["OI_REVN"] = 2, "Revision number of the table definition"
    hdu.header["DATE-OBS"] = hdr.get("date-obs", "")
    hdulist.append(hdu)

    # ------------------------------
    #          Save file
    # ------------------------------
    savedfile = os.path.join(datadir, oifits_file)
    hdulist.writeto(savedfile, overwrite=True)
    if verbose:
        rprint(f"[cyan]\n\n### OIFITS CREATED ({oifits_file}).", "cyan")

    if len(l_dic) == 1:
        l_dic = l_dic[0]
    return l_dic, savedfile

show(inputList, diffWl=False, ind_hole=None, vmin=0, vmax=1.05, cmax=180, setlog=False, pa=0, unit='arcsec', unit_cp='deg', snr=4, true_flag_v2=True, true_flag_t3=False)

Display one or more OIFITS data sets.

Parameters:

Name	Type	Description	Default
inputList	object, str, dict, or list	OIFITS filenames, dictionaries, or calibrated data objects.	required
diffWl	bool	Whether to color data by wavelength.	False
ind_hole	int or None	Hole used to select independent closure phases.	None
vmin	float	Visibility plot limits.	0, 1.05
vmax	float	Visibility plot limits.	0, 1.05
cmax	float	Closure-phase plot limit in degrees.	180
setlog	bool	Whether to use a logarithmic visibility axis.	False
pa	float	Position angle in degrees.	0
unit	str	Spatial-frequency unit.	"arcsec"
unit_cp	str	Closure-phase unit.	"deg"
snr	float	Signal-to-noise threshold for flags.	4
true_flag_v2	bool	Whether to compute squared-visibility and closure-phase flags.	True, False
true_flag_t3	bool	Whether to compute squared-visibility and closure-phase flags.	True, False

Returns:

Type	Description
Figure or None	Figure containing the observables, when data are available.

Source code in src/amical/oifits.py

def show(
    inputList,
    diffWl=False,
    ind_hole=None,
    vmin=0,
    vmax=1.05,
    cmax=180,
    setlog=False,
    pa=0,
    unit="arcsec",
    unit_cp="deg",
    snr=4,
    true_flag_v2=True,
    true_flag_t3=False,
):
    """Display one or more OIFITS data sets.

    Parameters
    ----------
    inputList : object, str, dict, or list
        OIFITS filenames, dictionaries, or calibrated data objects.
    diffWl : bool, default=False
        Whether to color data by wavelength.
    ind_hole : int or None, default=None
        Hole used to select independent closure phases.
    vmin, vmax : float, default=0, 1.05
        Visibility plot limits.
    cmax : float, default=180
        Closure-phase plot limit in degrees.
    setlog : bool, default=False
        Whether to use a logarithmic visibility axis.
    pa : float, default=0
        Position angle in degrees.
    unit : str, default="arcsec"
        Spatial-frequency unit.
    unit_cp : str, default="deg"
        Closure-phase unit.
    snr : float, default=4
        Signal-to-noise threshold for flags.
    true_flag_v2, true_flag_t3 : bool, default=True, False
        Whether to compute squared-visibility and closure-phase flags.

    Returns
    -------
    matplotlib.figure.Figure or None
        Figure containing the observables, when data are available.
    """
    import matplotlib.pyplot as plt

    if not isinstance(inputList, list):
        inputList = [inputList]

    if hasattr(inputList[0], "v2"):
        l_dic = [
            cal2dict(
                x,
                pa=pa,
                true_flag_v2=true_flag_v2,
                ind_hole=ind_hole,
                true_flag_t3=true_flag_t3,
                snr=snr,
            )
            for x in inputList
        ]
        print("\n -- SHOW -- Inputs are classes from amical.calibrate:")
        print("-> (Check true_flag_v2, true_flag_t3 and snr parameters)\n")
    elif isinstance(inputList[0], str):
        l_dic = [load(x) for x in inputList]
        print("Inputs are oifits filename.")
    elif isinstance(inputList[0], dict):
        l_dic = inputList
        print("Inputs are dict from amical.load.")
    else:
        print("Wrong input format.")
        return None

    dic_color = {}
    i_c = 0
    for dic in l_dic:
        filt = dic["info"]["FILT"]
        if filt not in dic_color.keys():
            dic_color[filt] = list_color[i_c]
            i_c += 1

    props = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}

    textstr = f"PA = {pa:2.1f} deg"

    fontsize = 14
    fig = plt.figure(figsize=(16, 5.5))
    ax1 = plt.subplot2grid((2, 6), (0, 0), rowspan=2, colspan=2)
    ax2 = plt.subplot2grid((2, 6), (0, 2), colspan=4)
    ax3 = plt.subplot2grid((2, 6), (1, 2), colspan=4)

    # Plot plan UV
    # -------
    ins = l_dic[0]["info"]["INSTRUME"]
    if "IFS" in ins and len(l_dic) > 1:
        l_bmax = _plot_UV_ifu(ax1, fig, l_dic)
    else:
        l_bmax = _plot_UV(ax1, l_dic, dic_color, diffWl=False)

    Bmax = np.max(l_bmax)
    ax1.axis([Bmax, -Bmax, -Bmax, Bmax])
    ax1.spines["left"].set_visible(False)
    ax1.spines["right"].set_visible(False)
    ax1.spines["bottom"].set_visible(False)
    ax1.spines["top"].set_visible(False)
    ax1.patch.set_facecolor("#f7f9fc")
    ax1.patch.set_alpha(1)
    ax1.xaxis.set_ticks_position("none")
    ax1.yaxis.set_ticks_position("none")
    if diffWl:
        handles, labels = ax1.get_legend_handles_labels()
        labels, handles = zip(
            *sorted(zip(labels, handles, strict=False), key=lambda t: t[0]),
            strict=False,
        )
        ax1.legend(handles, labels, loc="best", fontsize=9)

    plt.text(
        0.02,
        0.98,
        textstr,
        transform=ax1.transAxes,
        fontsize=13,
        verticalalignment="top",
        bbox=props,
    )

    unitlabel = {
        "m": "m",
        "rad": "rad$^{-1}$",
        "arcsec": "arcsec$^{-1}$",
        "lambda": r"M$\lambda$",
    }

    ax1.set_xlabel(rf"U [{unitlabel[unit]}]", fontsize=fontsize)
    ax1.set_ylabel(rf"V [{unitlabel[unit]}]", fontsize=fontsize)
    ax1.grid(alpha=0.2)

    # Plot V2
    # -------
    if "IFS" in ins and len(l_dic) > 1:
        max_f_vis = _plot_V2_ifu(ax2, l_dic)
    else:
        max_f_vis = _plot_V2(ax2, l_dic, dic_color, diffWl=diffWl)

    ax2.hlines(1, 0, 1.2 * np.max(max_f_vis), lw=1, color="k", alpha=0.2, ls="--")

    ax2.set_ylim([vmin, vmax])
    ax2.set_xlim([0, 1.2 * np.max(max_f_vis)])
    ax2.set_ylabel(r"$V^2$", fontsize=fontsize)
    ax2.spines["left"].set_visible(False)
    ax2.spines["right"].set_visible(False)
    ax2.spines["bottom"].set_visible(False)
    ax2.spines["top"].set_visible(False)
    ax2.patch.set_facecolor("#f7f9fc")
    ax2.patch.set_alpha(1)
    ax2.xaxis.set_ticks_position("none")
    ax2.yaxis.set_ticks_position("none")

    if setlog:
        ax2.set_yscale("log")
    ax2.grid(which="both", alpha=0.2)

    # Plot CP
    # -------
    if unit_cp == "rad":
        conv_cp = np.pi / 180.0
        h1 = np.pi
    else:
        conv_cp = 1
        h1 = np.rad2deg(np.pi)

    cmin = -cmax

    if "IFS" in ins and len(l_dic) > 1:
        max_f_cp = _plot_CP_ifu(ax3, l_dic, conv_cp)
    else:
        max_f_cp = _plot_CP(ax3, l_dic, dic_color, conv_cp, diffWl=diffWl)

    ax3.hlines(h1, 0, 1.2 * np.max(max_f_cp), lw=1, color="k", alpha=0.2, ls="--")
    ax3.hlines(-h1, 0, 1.2 * np.max(max_f_cp), lw=1, color="k", alpha=0.2, ls="--")
    ax3.spines["left"].set_visible(False)
    ax3.spines["right"].set_visible(False)
    ax3.spines["bottom"].set_visible(False)
    ax3.spines["top"].set_visible(False)
    ax3.patch.set_facecolor("#f7f9fc")
    ax3.patch.set_alpha(1)
    ax3.xaxis.set_ticks_position("none")
    ax3.yaxis.set_ticks_position("none")
    ax3.set_xlabel("Spatial frequency [arcsec$^{-1}$]", fontsize=fontsize)
    ax3.set_ylabel(rf"Clos. $\phi$ [{unit_cp}]", fontsize=fontsize)
    ax3.axis([0, 1.2 * np.max(max_f_cp), cmin * conv_cp, cmax * conv_cp])
    ax3.grid(which="both", alpha=0.2)

    plt.subplots_adjust(
        top=0.974, bottom=0.1, left=0.05, right=0.99, hspace=0.127, wspace=0.35
    )

    return fig

simbad_enabled()

Return False if SIMBAD queries are switched off with the AMICAL_NO_SIMBAD environment variable (e.g. on HPC compute nodes without internet access).

Source code in src/amical/oifits.py

def simbad_enabled():
    """Return False if SIMBAD queries are switched off with the AMICAL_NO_SIMBAD
    environment variable (e.g. on HPC compute nodes without internet access)."""
    flag = os.environ.get("AMICAL_NO_SIMBAD", "").strip().lower()
    return flag in ("", "0", "false", "no")

wrap_raw(bs)

Wrap raw extraction observables in a calibration-compatible object.

Parameters:

Name	Type	Description	Default
bs	Munch	Raw observables returned by the extraction pipeline.	required

Returns:

Type	Description
Munch	Object compatible with calibrated-data OIFITS saving.

Source code in src/amical/oifits.py

def wrap_raw(bs):
    """Wrap raw extraction observables in a calibration-compatible object.

    Parameters
    ----------
    bs : munch.Munch
        Raw observables returned by the extraction pipeline.

    Returns
    -------
    munch.Munch
        Object compatible with calibrated-data OIFITS saving.
    """
    u1 = bs.u[bs.mask.bs2bl_ix[0, :]]
    v1 = bs.v[bs.mask.bs2bl_ix[0, :]]
    u2 = bs.u[bs.mask.bs2bl_ix[1, :]]
    v2 = bs.v[bs.mask.bs2bl_ix[1, :]]

    fake_cal = {
        "vis2": bs.vis2,
        "e_vis2": bs.e_vis2,
        "cp": bs.cp,
        "e_cp": bs.e_cp,
        "u": bs.u,
        "v": bs.v,
        "wl": bs.wl,
        "u1": u1,
        "v1": v1,
        "u2": u2,
        "v2": v2,
        "raw_t": bs,
        "raw_c": bs,
    }

    return dict2class(fake_cal)