Skip to content

Update data loader to the latest PFS datamodel (as of Nov 2024) #1202

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Update data loader to the latest PFS datamodel (as of Nov 2024) #1202

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 58 additions & 35 deletions specutils/io/default_loaders/subaru_pfs_spec.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,26 @@

https://github.com/Subaru-PFS/datamodel/blob/master/datamodel.txt
"""

import os
import re

from astropy.units import Unit
from astropy.nddata import StdDevUncertainty

import numpy as np
from astropy.nddata import StdDevUncertainty
from astropy.units import Unit

from ...spectra import Spectrum1D
from ..registers import data_loader
from ..parsing_utils import _fits_identify_by_name, read_fileobj_or_hdulist
from ..registers import data_loader

__all__ = ['identify_pfs_spec', 'pfs_spec_loader']
__all__ = ["identify_pfs_spec", "pfs_spec_loader"]

# This RE matches the file name pattern defined in Subaru-PFS' datamodel.txt :
# "pfsObject-%05d-%s-%3d-%08x-%02d-0x%08x.fits" % (tract, patch, catId, objId,
# nVisit % 100, pfsVisitHash)
_spec_pattern = re.compile(r'pfsObject-(?P<tract>\d{5})-(?P<patch>.{3})-'
r'(?P<catId>\d{3})-(?P<objId>\d{8})-'
r'(?P<nVisit>\d{2})-(?P<pfsVisitHash>0x\w{8})'
r'\.fits')
# "pfsObject-%05d-%05d-%s-%016x-%03d-0x%016x.fits" % (catId, tract, patch, objId,
# nVisit % 1000, pfsVisitHash)
_spec_pattern = re.compile(
r"pfsObject-(?P<catId>\d{5})-(?P<tract>\d{5})-(?P<patch>.{3})-(?P<objId>\d{16})-(?P<nVisit>\d{3})-(?P<pfsVisitHash>0x\w{16})\.fits"
)


def identify_pfs_spec(origin, *args, **kwargs):
Expand All @@ -35,7 +34,9 @@ def identify_pfs_spec(origin, *args, **kwargs):


@data_loader(
label="Subaru-pfsObject", identifier=identify_pfs_spec, extensions=['fits'],
label="Subaru-pfsObject",
identifier=identify_pfs_spec,
extensions=["fits"],
priority=10,
)
def pfs_spec_loader(file_obj, **kwargs):
Expand Down Expand Up @@ -64,28 +65,50 @@ def pfs_spec_loader(file_obj, **kwargs):

with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist:
header = hdulist[0].header
meta = {'header': header,
'tract': m['tract'],
'patch': m['patch'],
'catId': m['catId'],
'objId': m['objId'],
'nVisit': m['nVisit'],
'pfsVisitHash': m['pfsVisitHash']}

# spectrum is in HDU 2
data = hdulist[2].data['flux']
unit = Unit('nJy')

error = hdulist[2].data['fluxVariance']
meta = {
"header": header,
"tract": m["tract"],
"patch": m["patch"],
"catId": m["catId"],
"objId": m["objId"],
"nVisit": m["nVisit"],
"pfsVisitHash": m["pfsVisitHash"],
}

# PFS datamodel: https://github.com/Subaru-PFS/datamodel/blob/master/datamodel.txt
#
# HDU #8 FLUXTABLE Binary table [FITS BINARY TABLE] NOBS*NROW
# Columns for:
# wavelength in units of nm (vacuum) [64-bit FLOAT]
# intensity in units of nJy [FLOAT]
# intensity error in same units as intensity [FLOAT]
# mask [32-bit INT]
#
# In the datamodel, the FLUXTABLE is the 8th HDU, but in fact it's in the 10th HDU with EXTNAME of FLUX_TABLE.
#
# NOTE: No backward compatibility is guaranteed for the FLUXTABLE format.
data = hdulist[10].data["flux"]
unit = Unit("nJy")

error = hdulist[10].data["error"]
uncertainty = StdDevUncertainty(np.sqrt(error))

wave = hdulist[2].data['lambda']
wave_unit = Unit('nm')

mask = hdulist[2].data['mask'] != 0

return Spectrum1D(flux=data * unit,
spectral_axis=wave * wave_unit,
uncertainty=uncertainty,
meta=meta,
mask=mask)
wave = hdulist[10].data["wavelength"]
wave_unit = Unit("nm")

# NOTE: REFLINE mask is the 15th bit in the mask bits, and can be ignored.
# TODO: Mask condition needs to be refined once more information becomes available.
mask = ~np.logical_or(
# mask==0: good data
hdulist[10].data["mask"] == 0,
# mask==2**REFLINE (32768) can be ignored
hdulist[10].data["mask"] == 2 ** (hdulist[10].header["MP_REFLINE"]),
)

return Spectrum1D(
flux=data * unit,
spectral_axis=wave * wave_unit,
uncertainty=uncertainty,
meta=meta,
mask=mask,
)
Loading