Enable MEF extraction for raw data ingest

python/lsst/obs/cfht/__init__.py

@@ -22,3 +22,4 @@
 from .version import *
 from .megacamMapper import *
 from ._instrument import MegaPrime
+from .ingest import MegaPrimeRawIngestTask

python/lsst/obs/cfht/ingest.py

@@ -20,12 +20,17 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
-__all__ = ["MegacamParseTask"]
+__all__ = ["MegacamParseTask", "MegaPrimeRawIngestTask"]
 
 import re
 
+import lsst.obs.base
+from lsst.obs.base.ingest import RawFileData
+from astro_metadata_translator import fix_header
+
 from lsst.pipe.tasks.ingest import ParseTask
 import lsst.pex.exceptions
+from ._instrument import MegaPrime
 
 filters = {'u.MP9301': 'u',
            'u.MP9302': 'u2',
@@ -41,6 +46,34 @@
            }
 
 
+class MegaPrimeRawIngestTask(lsst.obs.base.RawIngestTask):
+    """Task for ingesting raw MegaPrime multi-extension FITS data into Gen3.
+    """
+    def extractMetadata(self, filename: str) -> RawFileData:
+        datasets = []
+        fitsData = lsst.afw.fits.Fits(filename, "r")
+
+        # NOTE: The primary header (HDU=0) does not contain detector data.
+        for i in range(1, fitsData.countHdus()):
+            fitsData.setHdu(i)
+            header = fitsData.readMetadata()
+            if not header["EXTNAME"].startswith("ccd"):
+                continue
+            fix_header(header)
+            datasets.append(self._calculate_dataset_info(header, filename))
+
+        # The data model currently assumes that whilst multiple datasets
+        # can be associated with a single file, they must all share the
+        # same formatter.
+        instrument = MegaPrime()
+        FormatterClass = instrument.getRawFormatter(datasets[0].dataId)
+
+        self.log.info(f"Found images for {len(datasets)} detectors in {filename}")
+        return RawFileData(datasets=datasets, filename=filename,
+                           FormatterClass=FormatterClass,
+                           instrumentClass=type(instrument))
+
+
 class MegacamParseTask(ParseTask):
 
     def translate_ccd(self, md):

python/lsst/obs/cfht/rawFormatter.py

@@ -19,21 +19,121 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
+import lsst.afw.image
+import lsst.afw.fits
 from lsst.obs.base import FitsRawFormatterBase
-from astro_metadata_translator import MegaPrimeTranslator
+from astro_metadata_translator import MegaPrimeTranslator, fix_header
+import logging
 
 from .cfhtFilters import MEGAPRIME_FILTER_DEFINITIONS
 from ._instrument import MegaPrime
 
 __all__ = ("MegaPrimeRawFormatter",)
 
+log = logging.getLogger(__name__)
+
 
 class MegaPrimeRawFormatter(FitsRawFormatterBase):
     """Gen3 Butler Formatters for HSC raw data.
+
+    MegaPrime uses multi-extension FITS files.
     """
 
     translatorClass = MegaPrimeTranslator
     filterDefinitions = MEGAPRIME_FILTER_DEFINITIONS
 
     def getDetector(self, id):
         return MegaPrime().getCamera()[id]
+
+    def _toExtName(self, detectorId):
+        """Return the extension name associated with the given detector.
+
+        Parameters
+        ----------
+        detectorId : `int`
+            The detector ID to search for.
+
+        Returns
+        -------
+        name : `str`
+            The name of the extension associated with this detector.
+        """
+        return f"ccd{detectorId:02d}"
+
+    def _scanHdus(self, filename, detectorId):
+        """Scan through a file for the HDU containing data from one detector.
+
+        Parameters
+        ----------
+        filename : `str`
+            The file to search through.
+        detectorId : `int`
+            The detector id to search for.
+
+        Returns
+        -------
+        index : `int`
+            The index of the HDU with the requested data.
+        metadata: `lsst.daf.base.PropertyList`
+            The metadata read from the header for that detector id.
+
+        Raises
+        ------
+        ValueError
+            Raised if detectorId is not found in any of the file HDUs
+        """
+        fitsData = lsst.afw.fits.Fits(filename, 'r')
+        # NOTE: The primary header (HDU=0) does not contain detector data.
+        extname = self._toExtName(detectorId)
+        for i in range(1, fitsData.countHdus()):
+            fitsData.setHdu(i)
+            metadata = fitsData.readMetadata()
+            if metadata.get("EXTNAME") == extname:
+                log.debug("Found detector %s in extension %s", detectorId, i)
+                return i, metadata
+        else:
+            raise ValueError(f"Did not find detectorId={detectorId} in any HDU of {filename}.")
+
+    def _determineHDU(self, detectorId):
+        """Determine the correct HDU number for a given detector id.
+
+        Parameters
+        ----------
+        detectorId : `int`
+            The detector id to search for.
+
+        Returns
+        -------
+        index : `int`
+            The index of the HDU with the requested data.
+        metadata : `lsst.daf.base.PropertyList`
+            The metadata read from the header for that detector id.
+
+        Raises
+        ------
+        ValueError
+            Raised if detectorId is not found in any of the file HDUs
+        """
+        filename = self.fileDescriptor.location.path
+        # We start by assuming that ccdN is HDU N+1
+        index = detectorId + 1
+        metadata = lsst.afw.fits.readMetadata(filename, index)
+
+        # There may be two EXTNAME headers but the second one is the one
+        # we want (the first indicates compression).
+        if metadata.get("EXTNAME") == self._toExtName(detectorId):
+            return index, metadata
+
+        log.info("Did not find detector=%s at expected HDU=%s in %s: scanning through all HDUs.",
+                 detectorId, index, filename)
+        return self._scanHdus(filename, detectorId)
+
+    def readMetadata(self):
+        # Headers are duplicated so no need to merge with primary
+        index, metadata = self._determineHDU(self.dataId["detector"])
+        fix_header(metadata)
+        return metadata
+
+    def readImage(self):
+        index, metadata = self._determineHDU(self.dataId["detector"])
+        return lsst.afw.image.ImageI(self.fileDescriptor.location.path, index)

tests/test_ingest.py

@@ -42,12 +42,13 @@ class MegaPrimeIngestTestCase(IngestTestBase, lsst.utils.tests.TestCase):
     curatedCalibrationDatasetTypes = ()
     ingestDir = os.path.dirname(__file__)
     instrumentClassName = "lsst.obs.cfht.MegaPrime"
+    rawIngestTask = "lsst.obs.cfht.MegaPrimeRawIngestTask"
 
     @property
     def file(self):
         return os.path.join(testDataDirectory, "DATA/raw/08BL05/w2.+2+2/2008-11-01/i2/1038843o.fits.fz")
 
-    dataIds = [dict(instrument="MegaPrime", exposure=1038843, detector=0)]
+    dataIds = [dict(instrument="MegaPrime", exposure=1038843, detector=i) for i in range(36)]
 
     @property
     def visits(self):