Merge pull request #102 from lsst/tickets/DM-12370

DM-12370: Add piecewise TransmissionCurve for coadds.

include/lsst/meas/algorithms/CoaddTransmissionCurve.h

@@ -0,0 +1,57 @@
+// -*- lsst-c++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2018 LSST/AURA.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program.  If not,
+ * see <http://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#ifndef LSST_MEAS_ALGORITHMS_CoaddTransmissionCurve_h_INCLUDED
+#define LSST_MEAS_ALGORITHMS_CoaddTransmissionCurve_h_INCLUDED
+
+#include "lsst/afw/image/TransmissionCurve.h"
+#include "lsst/afw/table/Exposure.h"
+
+namespace lsst { namespace meas { namespace algorithms {
+
+/**
+ *  Create a TransmissionCurve that represents the effective throughput on a coadd.
+ *
+ *  @param[in]  coaddWcs       WCS that relates the coadd coordinate system to the sky.
+ *  @param[in]  inputSensors   A catalog containing the WCSs, bounding boxes and polygons,
+ *                             coaddition weights (in a field called 'weight'), and
+ *                             TransmissionCurves of the sensor-level images that went
+ *                             into the coadd.
+ *
+ *  @throws NotFoundError      Thrown if the 'weight' field does not exist in the schema.
+ *
+ *  @throws InvalidParameterError   Thrown if one or more inputs do not have a TransmissionCurve or
+ *                                  a Wcs (ValidPolygons may be null to
+ *                                  indicate no spatial restrictions other
+ *                                  than the bounding box).
+ *
+ *  @returns a new TransmissionCurve object.
+ */
+std::shared_ptr<afw::image::TransmissionCurve const> makeCoaddTransmissionCurve(
+    std::shared_ptr<afw::image::Wcs const> coaddWcs,
+    afw::table::ExposureCatalog const & inputSensors
+);
+
+}}} // namespace lsst::meas::algorithms
+
+#endif // !LSST_MEAS_ALGORITHMS_CoaddTransmissionCurve_h_INCLUDED

python/lsst/meas/algorithms/SConscript

@@ -4,6 +4,7 @@ scripts.BasicSConscript.pybind11(["binnedWcs",
                                   "cr",
                                   "coaddBoundedField",
                                   "coaddPsf/coaddPsf",
+                                  "coaddTransmissionCurve",
                                   "doubleGaussianPsf",
                                   "imagePsf",
                                   "interp",

python/lsst/meas/algorithms/__init__.py

@@ -39,6 +39,7 @@
 from .spatialModelPsf import *
 from .warpedPsf import *
 from .coaddPsf import *
+from .coaddTransmissionCurve import *
 from .doubleGaussianPsf import *
 
 from .defects import *

python/lsst/meas/algorithms/coaddTransmissionCurve.cc

@@ -0,0 +1,41 @@
+/*
+ * LSST Data Management System
+ * Copyright 2018 LSST/AURA.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ * See the COPYRIGHT file
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program.  If not,
+ * see <https://www.lsstcorp.org/LegalNotices/>.
+ */
+#include "pybind11/pybind11.h"
+
+#include "lsst/meas/algorithms/CoaddTransmissionCurve.h"
+#include "lsst/afw/image/Wcs.h"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+namespace lsst { namespace meas { namespace algorithms {
+
+PYBIND11_PLUGIN(coaddTransmissionCurve) {
+    py::module mod("coaddTransmissionCurve");
+    py::module::import("lsst.afw.image");
+    py::module::import("lsst.afw.table");
+    mod.def("makeCoaddTransmissionCurve", &makeCoaddTransmissionCurve, "coaddWcs"_a, "inputSensors"_a);
+    return mod.ptr();
+}
+
+}}} // lsst::meas::algorithms

python/lsst/meas/algorithms/testUtils.py

@@ -22,7 +22,7 @@
 #
 from __future__ import absolute_import, division, print_function
 
-import numpy
+import numpy as np
 
 import lsst.afw.image as afwImage
 import lsst.afw.geom as afwGeom
@@ -68,9 +68,9 @@ def plantSources(bbox, kwid, sky, coordList, addPoissonNoise=True):
 
     # add Poisson noise
     if (addPoissonNoise):
-        numpy.random.seed(seed=1)  # make results reproducible
+        np.random.seed(seed=1)  # make results reproducible
         imgArr = img.getArray()
-        imgArr[:] = numpy.random.poisson(imgArr)
+        imgArr[:] = np.random.poisson(imgArr)
 
     # bundle into a maskedimage and an exposure
     mask = afwImage.Mask(bbox)
@@ -84,3 +84,47 @@ def plantSources(bbox, kwid, sky, coordList, addPoissonNoise=True):
     exposure.setPsf(psf)
 
     return exposure
+
+
+def makeRandomTransmissionCurve(rng, minWavelength=4000.0, maxWavelength=7000.0, nWavelengths=200,
+                                maxRadius=80.0, nRadii=30, perturb=0.05):
+    """Create a random TransmissionCurve with nontrivial spatial and
+    wavelength variation.
+
+    Parameters
+    ----------
+    rng : numpy.random.RandomState
+        Random number generator.
+    minWavelength : float
+        Average minimum wavelength for generated TransmissionCurves (will be
+        randomly perturbed).
+    maxWavelength : float
+        Average maximum wavelength for generated TransmissionCurves (will be
+        randomly perturbed).
+    nWavelengths : int
+        Number of samples in the wavelength dimension.
+    maxRadius : float
+        Average maximum radius for spatial variation (will be perturbed).
+    nRadii : int
+        Number of samples in the radial dimension.
+    perturb: float
+        Fraction by which wavelength and radius bounds should be randomly
+        perturbed.
+    """
+    dWavelength = maxWavelength - minWavelength
+
+    def perturbed(x, s=perturb*dWavelength):
+        return x + 2.0*s*(rng.rand() - 0.5)
+
+    wavelengths = np.linspace(perturbed(minWavelength), perturbed(maxWavelength), nWavelengths)
+    radii = np.linspace(0.0, perturbed(maxRadius, perturb*maxRadius), nRadii)
+    throughput = np.zeros(wavelengths.shape + radii.shape, dtype=float)
+    # throughput will be a rectangle in wavelength, shifting to higher wavelengths and shrinking
+    # in height with radius, going to zero at all bounds.
+    peak0 = perturbed(0.9, 0.05)
+    start0 = perturbed(minWavelength + 0.25*dWavelength)
+    stop0 = perturbed(minWavelength + 0.75*dWavelength)
+    for i, r in enumerate(radii):
+        mask = np.logical_and(wavelengths >= start0 + r, wavelengths <= stop0 + r)
+        throughput[mask, i] = peak0*(1.0 - r/1000.0)
+    return afwImage.TransmissionCurve.makeRadial(throughput, wavelengths, radii)