specutils icon indicating copy to clipboard operation
specutils copied to clipboard

Add IRTF SpeX loader

Open kelle opened this issue 4 months ago • 0 comments

Prism loader already written in AstrodbKit.

test file: filename = "https://s3.amazonaws.com/bdnyc/SpeX/Prism/U10013_SpeX.fits"

def _identify_spex(filename):
    """
    Check whether the given file is a SpeX data product.
    """
    try:
        with fits.open(filename, memmap=False) as hdulist:
            return "spex" in hdulist[0].header["INSTRUME"].lower() and "irtf" in hdulist[0].header["TELESCOP"].lower()
    except Exception:  # pylint: disable=broad-except,
        return False


def identify_spex_prism(origin, *args, **kwargs):
    """
    Confirm this is a SpeX Prism FITS file.
    See FITS keyword reference at http://irtfweb.ifa.hawaii.edu/~spex/observer/
    Notes: GRAT has values of: ShortXD, Prism, LXD_long, LXD_short, SO_long, SO_short
    """
    is_spex = _identify_spex(args[0])
    if is_spex:
        with fits.open(args[0], memmap=False) as hdulist:
            return (
                isinstance(args[0], str)
                and os.path.splitext(args[0].lower())[1] == ".fits"
                and is_spex
                and ("lowres" in hdulist[0].header["GRAT"].lower() or "prism" in hdulist[0].header["GRAT"].lower())
            )
    else:
        return is_spex


@data_loader("Spex Prism", identifier=identify_spex_prism, extensions=["fits"], dtype=Spectrum)
def spex_prism_loader(filename, **kwargs):
    """Open a SpeX Prism file and convert it to a Spectrum1D object"""

    with fits.open(filename, **kwargs) as hdulist:
        header = hdulist[0].header

        tab = hdulist[0].data

        # Handle missing/incorrect units
        try:
            flux_unit = header["YUNITS"].replace("ergs", "erg ").strip()
            wave_unit = header["XUNITS"].replace("Microns", "um")
        except (KeyError, ValueError):
            # For now, assume some default units
            flux_unit = "erg"
            wave_unit = "um"

        wave, data = tab[0] * Unit(wave_unit), tab[1] * Unit(flux_unit)

        if tab.shape[0] == 3:
            uncertainty = StdDevUncertainty(tab[2])
        else:
            uncertainty = None

        meta = {"header": header}

    return Spectrum(flux=data, spectral_axis=wave, uncertainty=uncertainty, meta=meta)

and relevant tests from AstrodbKit

@mock.patch("astrodbkit.spectra.fits.open")
def test_identify_spex_prism(mock_fits_open, good_spex_file):
    mock_fits_open.return_value = good_spex_file

    filename = "https://s3.amazonaws.com/bdnyc/SpeX/Prism/U10013_SpeX.fits"
    assert identify_spex_prism("read", filename)
    filename = "I am not a valid spex prism file"
    assert not identify_spex_prism("read", filename)


@mock.patch("astrodbkit.spectra.fits.open")
def test_identify_spex(mock_fits_open, good_spex_file, bad_spex_file):
    mock_fits_open.return_value = good_spex_file
    assert _identify_spex("filename")
    mock_fits_open.return_value = bad_spex_file
    assert not _identify_spex("filename")


@mock.patch("astrodbkit.spectra.fits.open")
def test_load_spex_prism(mock_fits_open, good_spex_file, bad_spex_file):
    # Test good example
    mock_fits_open.return_value = good_spex_file
    spectrum = spex_prism_loader("filename")
    assert spectrum.unit == Unit("erg / (A cm2 s)")
    # Test bad example
    mock_fits_open.return_value = bad_spex_file
    spectrum = spex_prism_loader("filename")
    assert spectrum.unit == Unit("erg")

kelle avatar Jul 31 '25 14:07 kelle