Convert and Ingest Zhang18 Spectra

.gitignore

@@ -13,3 +13,4 @@ build/
 simple/version.py
 *.zip
 *.fits
+*sty2054_supplemental_files/

data/reference/Instruments.json

@@ -516,5 +516,12 @@
         "telescope": "GTC",
         "description": "Espectrografo Multiobjeto Infra-Rojo",
         "reference": null
+    },
+    {
+        "instrument": "Xshooter",
+        "mode": "Echelle-Smoothed",
+        "telescope": "ESO VLT",
+        "description": null,
+        "reference": null
     }
 ]

scripts/ingests/ingest_zhang18.py

@@ -0,0 +1,134 @@
+import logging
+import os
+from pathlib import Path
+from astropy.io import fits
+from astrodb_utils import AstroDBError, load_astrodb
+from astrodb_utils.spectra import ingest_spectrum
+from astrodb_utils.instruments import ingest_instrument
+from simple import REFERENCE_TABLES
+from specutils import Spectrum
+from datetime import datetime
+
+
+# set up logging 
+astrodb_logger = logging.getLogger("astrodb_utils")
+astrodb_logger.setLevel(logging.INFO)
+
+# logger for ingest_zhang18
+logger = logging.getLogger("ingest_zhang18")
+logger.setLevel(logging.INFO)
+
+# Load DB
+SAVE_DB = False
+RECREATE_DB = True
+SCHEMA_PATH = "simple/schema.yaml"
+db = load_astrodb(
+    "SIMPLE.sqlite",
+    recreatedb=RECREATE_DB,
+    reference_tables=REFERENCE_TABLES,
+    felis_schema=SCHEMA_PATH,
+)
+
+# Path
+file_path = "scripts/spectra_convert/zhang18/processed_spectra"
+
+# Add new mode for Xshooter smoothed spectra
+def add_mode():
+    """
+        Telescope: ESO VLT (existed)
+        Instrument: XShooter (existed)
+        mode: Echelle-smoothed (new mode for smoothed spectrum)
+    """
+    try:
+        ingest_instrument(
+            db,
+            telescope="ESO VLT",
+            instrument="Xshooter",
+            mode="Echelle-Smoothed",
+            raise_error=True
+        )
+        print("mode added successfully! ")
+    except AstroDBError as e:
+        astrodb_logger.error(f"Error adding mode: {e}")
+
+def modify_date(obs_date):
+    if not obs_date:
+        return None
+
+    return datetime.fromisoformat(obs_date) # example: 2004-04-17 04:40:11.761000
+
+
+# Ingest spectra --
+def ingest_zhang18():
+    ingested = 0
+    failed = 0
+    failed_file = []
+
+    fits_file = list(Path(file_path).glob("*.fits"))
+
+    for file in fits_file:
+        filename = file.name
+        print(f"Reading {filename}")
+        spectrum = Spectrum.read(file)
+
+        with fits.open(file) as hdul:
+            header = hdul[0].header
+
+        # Get all neccessary info to ingest
+        source_name = header.get("OBJECT")
+        access_url = filename
+
+        # Regime
+        if "NIR" in filename:
+            regime = "nir"
+        else:
+            regime = "optical"
+
+        # Mode
+        if "SMOOTHED" in filename:
+            mode = "Echelle-Smoothed"
+        elif "Xshooter" in filename and not "SMOOTHED" in filename:
+            mode = "Echelle"
+        else:
+            mode = "Missing"
+
+        # publication
+        ref = header.get("VOREF")
+        if ref == "10.1093/mnras/stw2438": # I
+            reference = "Zhan17.3040"
+        elif ref == "10.1093/mnras/stx350": # II
+            reference = "Zhan17"
+        elif ref == "10.1093/mnras/sty1352": # III
+            reference  = "Zhan18.1352"
+        elif ref == "10.1093/mnras/sty2054": # IV
+            reference = "Zhan18.2054"
+
+        # Date-OBS
+        obs_date = modify_date(header.get("DATE-OBS"))
+
+        # iNgest
+        try:
+            ingest_spectrum(
+                db,
+                source=source_name,
+                spectrum=access_url,
+                regime=regime,
+                mode=mode,
+                telescope=header.get("TELESCOP"),
+                instrument=header.get("INSTRUME"),
+                obs_date=obs_date,
+                reference=reference
+            ) 
+            logger.info(f"Successfully ingested spectrum")
+            ingested += 1
+        except AstroDBError as e:
+            logger.error(f"Error ingesting spectrum: {e}")
+            failed +=1
+            failed_file.append(filename)
+
+add_mode()
+ingest_zhang18()
+
+if SAVE_DB:
+    db.save_db(directory="data/")
+

scripts/spectra_convert/zhang18/add_header.py

@@ -0,0 +1,129 @@
+import os
+import numpy as np
+import pandas as pd
+from astropy.io import fits
+from astropy.time import Time
+from datetime import datetime
+import astropy.units as u
+from specutils import Spectrum
+from astrodb_utils.fits import add_missing_keywords, add_wavelength_keywords, check_header
+from astrodb_utils.spectra import check_spectrum_plottable
+from astroquery.simbad import Simbad
+from pathlib import Path
+
+# Path config
+path = "/Users/guanying/SIMPLE db/SIMPLE-db/scripts/spectra_convert/zhang18/processed_spectra"
+
+
+def get_paper_metadata(filename):
+    if "-I_SIMPLE.fits" in filename or "-I_SMOOTHED_SIMPLE.fits" in filename:
+        title = "Primeval very low-mass stars and brown dwarfs - I. Six new L subdwarfs, classification and atmospheric properties"
+        voref = "10.1093/mnras/stw2438"
+    elif "-II_SIMPLE.fits" in filename or "-II_SMOOTHED_SIMPLE.fits" in filename:
+        title = "Primeval very low-mass stars and brown dwarfs - II. The most metal-poor substellar object"
+        voref = "10.1093/mnras/stx350"
+    elif "-III_SIMPLE.fits" in filename or "-III_SMOOTHED_SIMPLE.fits" in filename:
+        title = "Primeval very low-mass stars and brown dwarfs - III. The halo transitional brown dwarfs"
+        voref = "10.1093/mnras/sty1352"
+    elif "-IV_SIMPLE.fits" in filename or "IV_SMOOTHED_SIMPLE.fits" in filename:
+        title = "Primeval very low-mass stars and brown dwarfs - IV. New L subdwarfs, Gaia astrometry, population properties, and a blue brown dwarf binary"
+        voref = "10.1093/mnras/sty2054"
+    else:
+        title = "Unknown Paper"
+        voref = "Unknown DOI"
+
+    return title, voref
+
+
+def add_header(path):
+    missing_telescop_instrument = []
+    missing_dateobs = []
+    file_proceed = 0
+
+
+    fits_files = list(Path(path).glob("*.fits"))
+
+    for fits_file in fits_files:
+        filename = fits_file.name
+        print(f"\nProcessing {filename}...")
+        spectrum = Spectrum.read(fits_file)
+
+        title, voref = get_paper_metadata(filename)
+
+        with fits.open(fits_file, mode="update") as hdul:
+            header = hdul[0].header
+
+            # add SIMPLE FITS headers
+            header["SIMPLE"] = (True, "Conforms to FITS standard")
+            header["VOPUB"] = ("SIMPLE Archive", "Publication of the spectrum")
+
+            # object name
+            object_name = " ".join(filename.split("_")[:2])
+            header["OBJECT"] = (object_name, "Name of the object")
+
+            # Query simbad for ra/dec
+            try:
+                result = Simbad.query_object(object_name)
+                if result is not None and len(result) > 0:
+                    header["RA_TARG"] = (result[0]["ra"], "[ra] Pointing position")
+                    header["DEC_TARG"] = (result[0]["dec"], "[dec] Pointing position")
+            except Exception as e:
+                print(f"  SIMBAD lookup failed: {e}")
+
+            # Remove deprecated keyword
+            if "RADECSYS" in header:
+                header.remove("RADECSYS", ignore_missing=True)
+                print("  Removed RADECSYS")
+
+            # Check TELESCOP + INSTRUME
+            if "TELESCOP" not in header or "INSTRUME" not in header:
+                missing_telescop_instrument.append(filename)
+
+            add_wavelength_keywords(header, spectrum.spectral_axis)
+
+            # Paper metadata
+            header["TITLE"] = (title, "Title of the paper")
+            header["VOREF"] = (voref, "DOI of the paper")
+            header["AUTHOR"] = ("Zhang, Z.H. et al.", "Original authors")
+            header["DATE"] = (datetime.now().strftime("%Y-%m-%dT%H:%M:%S"))
+            header["CONTRIB1"] = ("Guan Ying Goh")
+            header["TELAPSE"] = (header.get("EXPTIME"))
+            mjd_obs = header.get("MJD-OBS")
+            exptime = header.get("EXPTIME")
+            if mjd_obs is None or exptime is None:
+                print(f"  WARNING: MJD-OBS or EXPTIME missing in header for {filename}. Cannot calculate TMID.")
+            if mjd_obs is not None and exptime is not None:
+                tmid = mjd_obs + (exptime / 2) / (60 * 60 * 24)
+                header.set("TMID", tmid, "[d] MJD of exposure mid-point")
+            print("   Added paper metadata")
+
+            # Fix flux units
+            if "BUNIT" not in header:
+                header["BUNIT"] = "erg / (cm2 s Angstrom)"
+            elif header["BUNIT"] == "erg/s/cm2/Angstrom":
+                header["BUNIT"] = "erg / (cm2 s Angstrom)"
+
+            # date-obs check
+            if "DATE-OBS" not in header or header["DATE-OBS"] in ["", "UNKNOWN", None]:
+                missing_dateobs.append(filename)
+
+            # add_missing_keywords(header)
+            check_header(hdul[0].header)
+
+            hdul.flush()
+            print(f"  Finished {filename}")
+            file_proceed += 1
+
+    # Summary
+    print(f"\File processed: {file_proceed}/{len(fits_files)}")
+    if missing_telescop_instrument:
+        print("\nFiles missing TELESCOP or INSTRUME:")
+        for f in missing_telescop_instrument:
+            print(f"  - {f}")
+
+    if missing_dateobs:
+        print("\nFiles missing DATE-OBS:")
+        for f in missing_dateobs:
+            print(f"  - {f}")
+
+add_header(path)