DM-10155: Implement in-place SourceCatalog operators in PhotoCalib

include/lsst/afw/image/PhotoCalib.h

@@ -250,8 +250,6 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
      *                            "PsfFlux_instFluxErr"
      * @param[in]  outField       The field to write the nJy and magnitude errors to.
      *                            Keys of the form "*_instFlux" and "*_instFluxErr" must exist in the schema.
-     *
-     * @warning Not implemented yet: See DM-10155.
      */
     void instFluxToNanojansky(afw::table::SourceCatalog &sourceCatalog, std::string const &instFluxField,
                               std::string const &outField) const;
@@ -330,16 +328,14 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
      * @param[in]  outField       The field to write the magnitudes and magnitude errors to.
      *                            Keys of the form "*_instFlux", "*_instFluxErr", *_mag", and "*_magErr"
      *                            must exist in the schema.
-     *
-     * @warning Not implemented yet: See DM-10155.
      */
     void instFluxToMagnitude(afw::table::SourceCatalog &sourceCatalog, std::string const &instFluxField,
                              std::string const &outField) const;
 
     /**
      * Return a flux calibrated image, with pixel values in nJy.
      *
-     * Mask pixels are propogated directly from the input image.
+     * Mask pixels are propagated directly from the input image.
      *
      * @param maskedImage The masked image to calibrate.
      * @param includeScaleUncertainty Include the uncertainty on the calibration in the resulting variance?
@@ -349,6 +345,29 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
     MaskedImage<float> calibrateImage(MaskedImage<float> const &maskedImage,
                                       bool includeScaleUncertainty = true) const;
 
+    /**
+     * Return a flux calibrated catalog, with new `_flux`, `_fluxErr`, `_mag`, and `_magErr` fields.
+     *
+     * If the input catalog already has `*_flux` and `*_mag` fields matching `instFluxFields`, they will be
+     * replaced with the new fields.
+     *
+     * @param catalog The catalog to compute calibrated fluxes of. This catalog is exactly reproduced in the
+     * output catalog.
+     * @param instFluxFields The fields to calibrate. If empty, every field ending with `_instFlux` will have
+     * its corresponding calibrated fields produced in the output catalog. If `_instFluxErr` also exists, it
+     * will be used to compute the `_fluxErr` and `_magErr` fields too.
+     *
+     * @return The calibrated catalog, with new flux and magnitude (and their errors) fields.
+     *
+     * @throws lsst::pex::exceptions::NotFoundError if any item in `instfluxFields` does not have a
+     * corresponding `*_instFlux` field in catalog.schema.
+     */
+    afw::table::SourceCatalog calibrateCatalog(afw::table::SourceCatalog const &catalog,
+                                               std::vector<std::string> const &instFluxFields) const;
+
+    /// @overload calibrateCatalog(table::SourceCatalog const &, std::vector<std::string> const &) const;
+    afw::table::SourceCatalog calibrateCatalog(afw::table::SourceCatalog const &catalog) const;
+
     /**
      * Convert AB magnitude to instFlux (ADU).
      *
@@ -362,6 +381,7 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
      * @returns    Source instFlux in ADU.
      */
     double magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const;
+
     /// @overload magnitudeToInstFlux(double, lsst::geom::Point<double, 2> const &) const;
     double magnitudeToInstFlux(double magnitude) const;
 
@@ -555,6 +575,15 @@ class PhotoCalib : public table::io::PersistableFacade<PhotoCalib>, public table
      * Helper function to manage constant vs. non-constant PhotoCalibs
      */
     double evaluate(lsst::geom::Point<double, 2> const &point) const;
+    /**
+     * Return the calibration evaluated at a sequence of points.
+     */
+    ndarray::Array<double, 1> evaluateArray(ndarray::Array<double, 1> const &xx,
+                                            ndarray::Array<double, 1> const &yy) const;
+    /**
+     * Return the calibration evaluated at the centroids of a SourceCatalog.
+     */
+    ndarray::Array<double, 1> evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const;
 
     /// Returns the spatially-constant calibration (for setting _calibrationMean)
     double computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const;
@@ -589,8 +618,8 @@ std::shared_ptr<PhotoCalib> makePhotoCalibFromMetadata(daf::base::PropertySet &m
  *
  * @param instFluxMag0 The instrumental flux at zero magnitude. If 0, the resulting `PhotoCalib` will have
  *                     infinite calibrationMean and non-finite (inf or NaN) calibrationErr.
- * @param instFluxMag0Err The instrumental flux at zero magnitude error. If 0, the resulting `PhotoCalib` will
- *                        have 0 calibrationErr.
+ * @param instFluxMag0Err The instrumental flux at zero magnitude error. If 0, the resulting `PhotoCalib`
+ * will have 0 calibrationErr.
  *
  * @returns Pointer to the constructed PhotoCalib.
  */

python/lsst/afw/image/photoCalib.cc

@@ -159,6 +159,14 @@ PYBIND11_MODULE(photoCalib, mod) {
     cls.def("calibrateImage", &PhotoCalib::calibrateImage, "maskedImage"_a,
             "includeScaleUncertainty"_a = true);
 
+    cls.def("calibrateCatalog",
+            py::overload_cast<afw::table::SourceCatalog const &, std::vector<std::string> const &>(
+                    &PhotoCalib::calibrateCatalog, py::const_),
+            "maskedImage"_a, "fluxFields"_a);
+    cls.def("calibrateCatalog",
+            py::overload_cast<afw::table::SourceCatalog const &>(&PhotoCalib::calibrateCatalog, py::const_),
+            "maskedImage"_a);
+
     /* Operators */
     cls.def("__eq__", &PhotoCalib::operator==, py::is_operator());
     cls.def("__ne__", &PhotoCalib::operator!=, py::is_operator());

src/image/PhotoCalib.cc

@@ -25,7 +25,7 @@
 #include "lsst/geom/Point.h"
 #include "lsst/afw/image/PhotoCalib.h"
 #include "lsst/afw/math/BoundedField.h"
-#include "lsst/afw/table/Source.h"
+#include "lsst/afw/table.h"
 #include "lsst/afw/table/io/CatalogVector.h"
 #include "lsst/afw/table/io/OutputArchive.h"
 #include "lsst/daf/base/PropertySet.h"
@@ -293,6 +293,86 @@ MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedIm
     return result;
 }
 
+afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
+                                                       std::vector<std::string> const &instFluxFields) const {
+    auto const &inSchema = catalog.getSchema();
+    afw::table::SchemaMapper mapper(inSchema, true);  // true: share the alias map
+    mapper.addMinimalSchema(inSchema);
+
+    using FieldD = afw::table::Field<double>;
+
+    struct Keys {
+        table::Key<double> instFlux;
+        table::Key<double> instFluxErr;
+        table::Key<double> flux;
+        table::Key<double> fluxErr;
+        table::Key<double> mag;
+        table::Key<double> magErr;
+    };
+
+    std::vector<Keys> keys;
+    keys.reserve(instFluxFields.size());
+    for (auto const &field : instFluxFields) {
+        Keys newKey;
+        newKey.instFlux = inSchema[inSchema.join(field, "instFlux")];
+        newKey.flux =
+                mapper.addOutputField(FieldD(inSchema.join(field, "flux"), "calibrated flux", "nJy"), true);
+        newKey.mag = mapper.addOutputField(
+                FieldD(inSchema.join(field, "mag"), "calibrated magnitude", "mag(AB)"), true);
+        try {
+            newKey.instFluxErr = inSchema.find<double>(inSchema.join(field, "instFluxErr")).key;
+            newKey.fluxErr = mapper.addOutputField(
+                    FieldD(inSchema.join(field, "fluxErr"), "calibrated flux uncertainty", "nJy"), true);
+            newKey.magErr = mapper.addOutputField(
+                    FieldD(inSchema.join(field, "magErr"), "calibrated magnitude uncertainty", "mag(AB)"),
+                    true);
+        } catch (pex::exceptions::NotFoundError &) {
+            ;  // Keys struct defaults to invalid keys; that marks the error as missing.
+        }
+        keys.emplace_back(newKey);
+    }
+
+    // Create the new catalog
+    afw::table::SourceCatalog output(mapper.getOutputSchema());
+    output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
+
+    auto calibration = evaluateCatalog(output);
+
+    // fill in the catalog values
+    int iRec = 0;
+    for (auto &rec : output) {
+        for (auto &key : keys) {
+            double instFlux = rec.get(key.instFlux);
+            double nanojansky = toNanojansky(instFlux, calibration[iRec]);
+            rec.set(key.flux, nanojansky);
+            rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
+            if (key.instFluxErr.isValid()) {
+                double instFluxErr = rec.get(key.instFluxErr);
+                rec.set(key.fluxErr, toNanojanskyErr(instFlux, instFluxErr, calibration[iRec],
+                                                     _calibrationErr, nanojansky));
+                rec.set(key.magErr,
+                        toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
+            }
+        }
+        ++iRec;
+    }
+
+    return output;
+}
+
+afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog) const {
+    std::vector<std::string> instFluxFields;
+    static std::string const SUFFIX = "_instFlux";
+    for (auto const &name : catalog.getSchema().getNames()) {
+        // Pick every field ending in "_instFlux", grabbing everything before that prefix.
+        if (name.size() > SUFFIX.size() + 1 &&
+            name.compare(name.size() - SUFFIX.size(), SUFFIX.size(), SUFFIX) == 0) {
+            instFluxFields.emplace_back(name.substr(0, name.size() - 9));
+        }
+    }
+    return calibrateCatalog(catalog, instFluxFields);
+}
+
 // ------------------- persistence -------------------
 
 namespace {
@@ -465,38 +545,66 @@ double PhotoCalib::evaluate(lsst::geom::Point<double, 2> const &point) const {
         return _calibration->evaluate(point);
 }
 
+ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1> const &xx,
+                                                    ndarray::Array<double, 1> const &yy) const {
+    if (_isConstant) {
+        ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
+        result.deep() = _calibrationMean;
+        return result;
+    } else {
+        return _calibration->evaluate(xx, yy);
+    }
+}
+
+ndarray::Array<double, 1> PhotoCalib::evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const {
+    ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
+    ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
+    size_t i = 0;
+    for (auto const &rec : sourceCatalog) {
+        auto point = rec.getCentroid();
+        xx[i] = point.getX();
+        yy[i] = point.getY();
+        ++i;
+    }
+    return evaluateArray(xx, yy);
+}
+
 void PhotoCalib::instFluxToNanojanskyArray(afw::table::SourceCatalog const &sourceCatalog,
                                            std::string const &instFluxField,
                                            ndarray::Array<double, 2, 2> result) const {
-    double instFlux, instFluxErr, nanojansky, calibration;
     auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
     auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
+
+    auto calibration = evaluateCatalog(sourceCatalog);
+    int i = 0;
     auto iter = result.begin();
     for (auto const &rec : sourceCatalog) {
-        instFlux = rec.get(instFluxKey);
-        instFluxErr = rec.get(instFluxErrKey);
-        calibration = evaluate(rec.getCentroid());
-        nanojansky = toNanojansky(instFlux, calibration);
+        double instFlux = rec.get(instFluxKey);
+        double instFluxErr = rec.get(instFluxErrKey);
+        double nanojansky = toNanojansky(instFlux, calibration[i]);
         (*iter)[0] = nanojansky;
-        (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
-        iter++;
+        (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration[i], _calibrationErr, nanojansky);
+        ++iter;
+        ++i;
     }
 }
 
 void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourceCatalog,
                                           std::string const &instFluxField,
                                           ndarray::Array<double, 2, 2> result) const {
-    double instFlux, instFluxErr, calibration;
     auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
     auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
+
+    auto calibration = evaluateCatalog(sourceCatalog);
     auto iter = result.begin();
+    int i = 0;
     for (auto const &rec : sourceCatalog) {
-        instFlux = rec.get(instFluxKey);
-        instFluxErr = rec.get(instFluxErrKey);
-        calibration = evaluate(rec.getCentroid());
-        (*iter)[0] = toMagnitude(instFlux, calibration);
-        (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
-        iter++;
+        double instFlux = rec.get(instFluxKey);
+        double instFluxErr = rec.get(instFluxErrKey);
+        (*iter)[0] = toMagnitude(instFlux, calibration[i]);
+        (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration[i], _calibrationErr);
+        ++iter;
+        ++i;
     }
 }