Check_spectrum_plottable fix for JWST data in Jy (#151)

* add optional format

* better test for flux units

* improvements to ingest_instrument

* check_spectrum_plottable test

* log message tweak

* name the main logger

* Spectrum1D -> Spectrum

* delete dead code

* require specutils>1.20.1

* remove check_spectrum_class helper

Author: kelle

astrodb_utils/__init__.py

@@ -15,7 +15,7 @@
 
 
 
-logger = logging.getLogger(__name__)
+logger = logging.getLogger("astrodb_utils")
 
 LOGFORMAT = logging.Formatter(
     "%(levelname)-8s - %(name)-15s - %(message)s")
@@ -30,7 +30,7 @@
 handler.setFormatter(LOGFORMAT)
 logger.addHandler(handler)
 
-logger.info("astrodb_utils logger initialized")
+logger.info(f"Logger initialized: {logger.name}")
 logger.info(f"Logger level: {logging.getLevelName(logger.getEffectiveLevel()) }")
 
 warnings.filterwarnings("ignore", module="astroquery.simbad")

astrodb_utils/spectra.py

@@ -4,38 +4,84 @@
 
 import astropy.units as u
 import numpy as np
-from specutils import Spectrum1D
+from specutils import Spectrum
 
-from astrodb_utils import AstroDBError
+from astrodb_utils import AstroDBError, exit_function
 
 matplotlib_check = importlib.util.find_spec("matplotlib")
 if matplotlib_check is not None:
     import matplotlib.pyplot as plt
 
 
-__all__ = [ 
-    "check_spectrum_plottable"
-    ]
+__all__ = ["check_spectrum_plottable"]
 
 logger = logging.getLogger(__name__)
 
 
-def check_spectrum_class(spectrum, raise_error=True):
-    try:
-        Spectrum1D.read(spectrum)
-        return True
-    except Exception as error_message:
-        msg = f"Unable to load file as Spectrum1D object:{spectrum}"
-        logger.debug(f"{error_message}")
-        if raise_error:
-            logger.error(msg)
-            raise AstroDBError(msg)
-        else:
-            logger.warning(msg)
-            return False
+def check_spectrum_plottable(
+    spectrum_path: str | Spectrum, raise_error: bool = True, show_plot: bool = False, format: str = None
+):
+    """
+    Check if spectrum is readable and plottable with specutils.
+    show_plot = True requires matplotlib to be installed.
+
+    Parameters
+    ----------
+    spectrum_path : str or Spectrum
+        Path to spectrum file or Spectrum object
+
+    raise_error : bool. Default=True
+        True: Raise error if spectrum is not plottable
+        False: Do not raise error if spectrum is not plottable. Log warning instead.
+
+    show_plot : bool. Default=False
+        True: Show plot of spectrum. Matplotlib must be installed.
+
+    format : str, optional
+        Format of the spectrum file. If not provided, the format will be inferred by specutils.
+
+    Returns
+    -------
+    bool
+        True: Spectrum is plottable
+        False: Spectrum is not plottable
 
+    """
+    # check if spectrum is a Spectrum object or a file path
+    # if it's a file path, check if it can be read as a Spectrum object    
+    if isinstance(spectrum_path, Spectrum):
+        spectrum = spectrum_path
+    elif isinstance(spectrum_path, str):
+        try:
+            spectrum = Spectrum.read(spectrum_path, format=format)
+        except Exception as error_message:
+            msg = f"Unable to load file as Spectrum object:{spectrum_path}:\n{error_message}"
+            exit_function(msg, raise_error=raise_error)
+    else:
+        msg = f"Input is not a valid path or Spectrum object: {spectrum_path}"
+        exit_function(msg, raise_error=raise_error)
+
+    # checking spectrum has good units
+    wave_unit_check = _check_spectrum_wave_units(spectrum, raise_error=raise_error)
+    if not wave_unit_check:
+        return False
+
+    flux_unit_check = _check_spectrum_flux_units(spectrum, raise_error=raise_error)
+    if not flux_unit_check:
+        return False
+
+    # check for NaNs
+    nan_check = _check_spectrum_not_nans(spectrum, raise_error=raise_error)
+    if not nan_check:
+        return False
 
-def check_spectrum_not_nans(spectrum, raise_error=True):
+    if show_plot:
+        _plot_spectrum(spectrum)
+
+    return True
+    
+
+def _check_spectrum_not_nans(spectrum, raise_error=True):
     nan_check: np.ndarray = ~np.isnan(spectrum.flux) & ~np.isnan(spectrum.spectral_axis)
     wave = spectrum.spectral_axis[nan_check]
     if not len(wave):
@@ -50,7 +96,7 @@ def check_spectrum_not_nans(spectrum, raise_error=True):
         return True
 
 
-def check_spectrum_wave_units(spectrum, raise_error=True):
+def _check_spectrum_wave_units(spectrum, raise_error=True):
     try:
         spectrum.spectral_axis.to(u.micron).value
         return True
@@ -83,31 +129,18 @@ def check_spectrum_wave_units(spectrum, raise_error=True):
             return False
 
 
-def check_spectrum_flux_units(spectrum, raise_error=True):
-    try:
-        spectrum.flux.to(u.erg / u.s / u.cm**2 / u.AA).value
+def _check_spectrum_flux_units(spectrum, raise_error=True):
+    expected_units = [
+        u.get_physical_type(u.erg / u.s / u.cm**2 / u.AA),
+        u.get_physical_type(u.Jy),
+    ]
+
+    unit_type = u.get_physical_type(spectrum.flux.unit)
+
+    if unit_type in expected_units:
         return True
-    except AttributeError as e:
-        logger.debug(f"{e}")
-        msg = f"Unable to parse flux: {spectrum}"
-        if raise_error:
-            logger.error(msg)
-            raise AstroDBError(msg)
-        else:
-            logger.warning(msg)
-            return False
-    except u.UnitConversionError as e:
-        logger.debug(f"{e}")
-        msg = f"Unable to convert flux to erg/s/cm^2/Angstrom:  {spectrum}"
-        if raise_error:
-            logger.error(msg)
-            raise AstroDBError(msg)
-        else:
-            logger.warning(msg)
-            return False
-    except ValueError as e:
-        logger.debug(f"{e}")
-        msg = f"Value error: {spectrum}:"
+    else:
+        msg = f"flux units are not expected: {spectrum.flux.unit}. Expecting {expected_units}."
         if raise_error:
             logger.error(msg)
             raise AstroDBError(msg)
@@ -116,67 +149,12 @@ def check_spectrum_flux_units(spectrum, raise_error=True):
             return False
 
 
-def plot_spectrum(spectrum):
+def _plot_spectrum(spectrum):
     if "matplotlib" in sys.modules:
         plt.plot(spectrum.spectral_axis, spectrum.flux)
-        plt.xlabel("Dispersion ({spectrum.spectral_axis.unit})")
-        plt.ylabel("Flux ({spectrum.flux.unit})")
+        plt.xlabel(f"Dispersion ({spectrum.spectral_axis.unit})")
+        plt.ylabel(f"Flux ({spectrum.flux.unit})")
         plt.show()
     else:
         msg = "To display the spectrum, matplotlib most be installed."
         logger.warning(msg)
-
-
-def check_spectrum_plottable(spectrum_path, raise_error=True, show_plot=False):
-    """
-    Check if spectrum is readable and plottable with specutils.
-    show_plot = True requires matplotlib to be installed.
-
-    Parameters
-    ----------
-    spectrum_path : str
-        Path to spectrum file
-
-    raise_error : bool. Default=True
-        True: Raise error if spectrum is not plottable
-        False: Do not raise error if spectrum is not plottable. Log warning instead.
-    
-    show_plot : bool. Default=False
-        True: Show plot of spectrum. Matplotlib must be installed.
-
-    Returns
-    -------
-    bool
-        True: Spectrum is plottable
-        False: Spectrum is not plotable
-
-    """
-    # load the spectrum and make sure it's readable as a Spectrum1D object, has units, is not all NaNs.
-    if isinstance(spectrum_path, Spectrum1D):
-        spectrum = spectrum_path
-        class_check = True
-    else:
-        class_check = check_spectrum_class(spectrum_path, raise_error=raise_error)
-        if not class_check:
-            return False
-        else:
-            spectrum = Spectrum1D.read(spectrum_path)
-
-    # checking spectrum has good units
-    wave_unit_check = check_spectrum_wave_units(spectrum, raise_error=raise_error)
-    if not wave_unit_check:
-        return False
-
-    flux_unit_check = check_spectrum_flux_units(spectrum, raise_error=raise_error)
-    if not flux_unit_check:
-        return False
-
-    # check for NaNs
-    nan_check = check_spectrum_not_nans(spectrum, raise_error=raise_error)
-    if not nan_check:
-        return False
-
-    if show_plot:
-        plot_spectrum(spectrum)
-
-    return True

astrodb_utils/utils.py

@@ -142,7 +142,7 @@ def ingest_instrument(db, *, telescope=None, instrument=None, mode=None):
         raise AstroDBError(msg)
 
     msg_search = f"Searching for {telescope}, {instrument}, {mode} in database"
-    logger.info(msg_search)
+    logger.debug(msg_search)
 
     # Search for the inputs in the database
     telescope_db = (
@@ -162,7 +162,7 @@ def ingest_instrument(db, *, telescope=None, instrument=None, mode=None):
 
     if len(telescope_db) == 1 and len(mode_db) == 1:
         msg_found = (
-            f"{telescope}, {instrument}, and {mode} are already in the database."
+            f"{telescope}-{instrument}-{mode} is already in the database. Nothing added."
         )
         logger.info(msg_found)
         return
@@ -174,10 +174,10 @@ def ingest_instrument(db, *, telescope=None, instrument=None, mode=None):
             with db.engine.connect() as conn:
                 conn.execute(db.Telescopes.insert().values(telescope_add))
                 conn.commit()
-            msg_telescope = f"{telescope} was successfully ingested in the database"
+            msg_telescope = f"{telescope} was successfully added to the Telescopes table."
             logger.info(msg_telescope)
         except sqlalchemy.exc.IntegrityError as e:  # pylint: disable=invalid-name
-            msg = "Telescope could not be ingested"
+            msg = f"{telescope} could not be added to the Telescopes table."
             logger.error(msg)
             raise AstroDBError(msg) from e
 
@@ -195,10 +195,10 @@ def ingest_instrument(db, *, telescope=None, instrument=None, mode=None):
             with db.engine.connect() as conn:
                 conn.execute(db.Instruments.insert().values(instrument_add))
                 conn.commit()
-            msg_instrument = f"{instrument} was successfully ingested in the database."
+            msg_instrument = f"{telescope}-{instrument}-{mode} was successfully added to the Instruments table."
             logger.info(msg_instrument)
         except sqlalchemy.exc.IntegrityError as e:  # pylint: disable=invalid-name
-            msg = "Instrument/Mode could not be ingested"
+            msg = f"{telescope}-{instrument}-{mode} could not be added to the Instruments table."
             logger.error(msg)
             raise AstroDBError(msg) from e
 

pyproject.toml

@@ -28,6 +28,7 @@ dependencies = [
     "tqdm",
     "ads",
     "dateparser",
+    "specutils>=1.20.1",
 ]
 dynamic = ["version"]
 

tests/conftest.py

@@ -52,7 +52,7 @@ def db():
         ignore_ads=True,
     )
 
-    ingest_publication(db, doi="10.1086/161442", reference="Prob83")
+    ingest_publication(db, doi="10.1086/161442", reference="Prob83", ignore_ads=True)
 
     return db