Merge pull request #142 from lsst/tickets/DM-19802

DM-19802: Fix jointcal ra/dec bounding box calculations

include/lsst/jointcal/Associations.h

@@ -35,6 +35,7 @@
 #include "lsst/afw/image/VisitInfo.h"
 #include "lsst/daf/base/PropertySet.h"
 #include "lsst/afw/geom/Box.h"
+#include "lsst/sphgeom/Circle.h"
 
 #include "lsst/jointcal/RefStar.h"
 #include "lsst/jointcal/FittedStar.h"
@@ -124,6 +125,11 @@ class Associations {
                         std::shared_ptr<afw::cameraGeom::Detector> detector, int visit, int ccd,
                         lsst::jointcal::JointcalControl const &control);
 
+    /**
+     * Add a pre-constructed ccdImage to the ccdImageList.
+     */
+    void addCcdImage(std::shared_ptr<CcdImage> const ccdImage) { ccdImageList.push_back(ccdImage); }
+
     //! incrementaly builds a merged catalog of all image catalogs
     void associateCatalogs(const double matchCutInArcsec = 0, const bool useFittedList = false,
                            const bool enlargeFittedList = true);
@@ -167,8 +173,13 @@ class Associations {
     //! Number of different bands in the input image list. Not implemented so far
     unsigned getNFilters() const { return 1; }
 
-    // Return the bounding box in (ra, dec) coordinates containing the whole catalog
-    const lsst::afw::geom::Box2D getRaDecBBox();
+    /**
+     * Return the bounding circle in on-sky (RA, Dec) coordinates containing all CcdImages.
+     *
+     * Requires that computeCommonTangentPoint() be called first, so that sensor bounding boxes can be
+     * transformed into the common tangent plane.
+     */
+    lsst::sphgeom::Circle computeBoundingCircle() const;
 
     /**
      * @brief      return the number of CcdImages with non-empty catalogs to-be-fit.

python/lsst/jointcal/associations.cc

@@ -27,6 +27,7 @@
 
 #include "lsst/jointcal/Associations.h"
 #include "lsst/jointcal/CcdImage.h"
+#include "lsst/sphgeom/Circle.h"
 
 namespace py = pybind11;
 using namespace pybind11::literals;
@@ -52,14 +53,14 @@ void declareAssociations(py::module &mod) {
     cls.def("nFittedStarsWithAssociatedRefStar", &Associations::nFittedStarsWithAssociatedRefStar);
 
     cls.def("createCcdImage", &Associations::createCcdImage);
+    cls.def("addCcdImage", &Associations::addCcdImage);
     cls.def("prepareFittedStars", &Associations::prepareFittedStars);
 
     cls.def("getCcdImageList", &Associations::getCcdImageList, py::return_value_policy::reference_internal);
     cls.def_property_readonly("ccdImageList", &Associations::getCcdImageList,
                               py::return_value_policy::reference_internal);
 
-    cls.def("getRaDecBBox", &Associations::getRaDecBBox);
-    cls.def_property_readonly("raDecBBox", &Associations::getRaDecBBox);
+    cls.def("computeBoundingCircle", &Associations::computeBoundingCircle);
 
     cls.def("getCommonTangentPoint", &Associations::getCommonTangentPoint);
     cls.def("setCommonTangentPoint", &Associations::setCommonTangentPoint);
@@ -68,6 +69,7 @@ void declareAssociations(py::module &mod) {
 
 PYBIND11_MODULE(associations, mod) {
     py::module::import("lsst.jointcal.ccdImage");
+    py::module::import("lsst.sphgeom");
     declareAssociations(mod);
 }
 }  // namespace

python/lsst/jointcal/jointcal.py

@@ -23,11 +23,11 @@
 import numpy as np
 import astropy.units as u
 
+import lsst.geom
 import lsst.utils
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 from lsst.afw.image import fluxErrFromABMagErr
-import lsst.afw.geom as afwGeom
 import lsst.pex.exceptions as pexExceptions
 import lsst.afw.table
 import lsst.meas.algorithms
@@ -518,20 +518,9 @@ def runDataRef(self, dataRefs, profile_jointcal=False):
 
         associations.computeCommonTangentPoint()
 
-        # Use external reference catalogs handled by LSST stack mechanism
-        # Get the bounding box overlapping all associated images
-        # ==> This is probably a bad idea to do it this way <== To be improved
-        bbox = associations.getRaDecBBox()
-        bboxCenter = bbox.getCenter()
-        center = afwGeom.SpherePoint(bboxCenter[0], bboxCenter[1], afwGeom.degrees)
-        bboxMax = bbox.getMax()
-        corner = afwGeom.SpherePoint(bboxMax[0], bboxMax[1], afwGeom.degrees)
-        radius = center.separation(corner).asRadians()
-
-        # Get astrometry_net_data path
-        anDir = lsst.utils.getPackageDir('astrometry_net_data')
-        if anDir is None:
-            raise RuntimeError("astrometry_net_data is not setup")
+        boundingCircle = associations.computeBoundingCircle()
+        center = lsst.geom.SpherePoint(boundingCircle.getCenter())
+        radius = lsst.geom.Angle(boundingCircle.getOpeningAngle().asRadians(), lsst.geom.radians)
 
         # Determine a default filter associated with the catalog. See DM-9093
         defaultFilter = filters.most_common(1)[0][0]
@@ -636,7 +625,7 @@ def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
             refCoordErr = self.config.astrometryReferenceErr
 
         associations.collectRefStars(refCat,
-                                     self.config.matchCut*afwGeom.arcseconds,
+                                     self.config.matchCut*lsst.geom.arcseconds,
                                      fluxField,
                                      refCoordinateErr=refCoordErr,
                                      rejectBadFluxes=reject_bad_fluxes)
@@ -693,7 +682,7 @@ def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radiu
             The name of the reference catalog flux field appropriate for ``filterName``.
         """
         skyCircle = refObjLoader.loadSkyCircle(center,
-                                               afwGeom.Angle(radius, afwGeom.radians),
+                                               radius,
                                                filterName)
 
         selected = referenceSelector.run(skyCircle.refCat)

python/lsst/jointcal/testUtils.py

@@ -50,7 +50,9 @@ def canRunTests():
 
 
 def createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeCcdId=12,
-                           photoCalibMean1=1e-2, photoCalibMean2=1.2e-2):
+                           photoCalibMean1=1e-2, photoCalibMean2=1.2e-2,
+                           fakeWcses=(None, None),
+                           fakeVisitInfos=(None, None)):
     """Return two fake ccdImages built on CFHT Megacam metadata.
 
     If ``num1 == num2``, the catalogs will align on-sky so each source will
@@ -73,6 +75,10 @@ def createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeCcdId=12,
     photoCalibMean1, photoCalibMean2: `float`, optional
         The mean photometric calibration to pass to each ccdImage construction.
         Note: this value is 1/instFluxMag0, so it should be less than 1.
+    fakeWcses : `list` [`lsst.afw.geom.SkyWcs`], optional
+        The SkyWcses to use instead of the ones read from disk.
+    fakeWcses : `list` [`lsst.afw.image.VisitInfo`], optional
+        The VisitInfos to use instead of the ones read from disk.
 
     Returns
     -------
@@ -107,19 +113,23 @@ def createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeCcdId=12,
     camera = butler.get('camera', visit=visit1)
 
     struct1 = createFakeCcdImage(butler, visit1, num1, fluxFieldName,
-                                 photoCalibMean=photoCalibMean1, photoCalibErr=1.0, fakeCcdId=fakeCcdId)
+                                 photoCalibMean=photoCalibMean1, photoCalibErr=1.0, fakeCcdId=fakeCcdId,
+                                 fakeWcs=fakeWcses[0], fakeVisitInfo=fakeVisitInfos[0])
     struct2 = createFakeCcdImage(butler, visit2, num2, fluxFieldName,
-                                 photoCalibMean=photoCalibMean2, photoCalibErr=5.0, fakeCcdId=fakeCcdId)
+                                 photoCalibMean=photoCalibMean2, photoCalibErr=5.0, fakeCcdId=fakeCcdId,
+                                 fakeWcs=fakeWcses[1], fakeVisitInfo=fakeVisitInfos[1])
 
     return lsst.pipe.base.Struct(camera=camera,
                                  catalogs=[struct1.catalog, struct2.catalog],
                                  ccdImageList=[struct1.ccdImage, struct2.ccdImage],
                                  bbox=struct1.bbox,
+                                 skyWcs=[struct1.skyWcs, struct2.skyWcs],
                                  fluxFieldName=fluxFieldName)
 
 
 def createFakeCcdImage(butler, visit, num, fluxFieldName,
-                       photoCalibMean=1e-2, photoCalibErr=1.0, fakeCcdId=12):
+                       photoCalibMean=1e-2, photoCalibErr=1.0, fakeCcdId=12,
+                       fakeWcs=None, fakeVisitInfo=None):
     """Create a fake CcdImage by making a fake catalog.
 
     Parameters
@@ -141,6 +151,10 @@ def createFakeCcdImage(butler, visit, num, fluxFieldName,
         Value to set for calibrationErr in the created PhotoCalib.
     fakeCcdId : `int`, optional
         Use this as the ccdId in the returned CcdImage.
+    fakeWcs : `lsst.afw.geom.SkyWcs`, optional
+        A SkyWcs to use instead of one read from disk.
+    fakeVisitInfo : `lsst.afw.image.VisitInfo`, optional
+        A VisitInfo to use instead of one read from disk.
 
     Returns
     -------
@@ -152,12 +166,13 @@ def createFakeCcdImage(butler, visit, num, fluxFieldName,
        - `ccdImage` : CcdImage containing the metadata and fake sources
            (`lsst.jointcal.CcdImage`).
        - `bbox` : Bounding Box of the image (`lsst.afw.geom.Box2I`).
+       - `skyWcs` : SkyWcs of the image (`lsst.afw.geom.SkyWcs`).
     """
     ccdId = 12  # we only have data for ccd=12
 
     dataId = dict(visit=visit, ccd=ccdId)
-    skyWcs = butler.get('calexp_wcs', dataId=dataId)
-    visitInfo = butler.get('calexp_visitInfo', dataId=dataId)
+    skyWcs = fakeWcs if fakeWcs is not None else butler.get('calexp_wcs', dataId=dataId)
+    visitInfo = fakeVisitInfo if fakeVisitInfo is not None else butler.get('calexp_visitInfo', dataId=dataId)
     bbox = butler.get('calexp_bbox', dataId=dataId)
     detector = butler.get('calexp_detector', dataId=dataId)
     filt = butler.get("calexp_filter", dataId=dataId).getName()
@@ -167,7 +182,7 @@ def createFakeCcdImage(butler, visit, num, fluxFieldName,
     ccdImage = lsst.jointcal.ccdImage.CcdImage(catalog, skyWcs, visitInfo, bbox, filt, photoCalib,
                                                detector, visit, fakeCcdId, fluxFieldName)
 
-    return lsst.pipe.base.Struct(catalog=catalog, ccdImage=ccdImage, bbox=bbox)
+    return lsst.pipe.base.Struct(catalog=catalog, ccdImage=ccdImage, bbox=bbox, skyWcs=skyWcs)
 
 
 def createFakeCatalog(num, bbox, fluxFieldName, skyWcs=None, refCat=False):

src/Associations.cc

@@ -46,6 +46,10 @@
 #include "lsst/afw/geom/Point.h"
 #include "lsst/afw/image/Calib.h"
 
+#include "lsst/sphgeom/LonLat.h"
+#include "lsst/sphgeom/Circle.h"
+#include "lsst/sphgeom/ConvexPolygon.h"
+
 namespace jointcal = lsst::jointcal;
 
 namespace {
@@ -85,6 +89,35 @@ void Associations::setCommonTangentPoint(lsst::afw::geom::Point2D const &commonT
     for (auto &ccdImage : ccdImageList) ccdImage->setCommonTangentPoint(_commonTangentPoint);
 }
 
+lsst::sphgeom::Circle Associations::computeBoundingCircle() const {
+    // Compute the frame on the common tangent plane that contains all input images.
+    Frame tangentPlaneFrame;
+
+    for (auto const &ccdImage : ccdImageList) {
+        Frame CTPFrame = ccdImage->getPixelToCommonTangentPlane()->apply(ccdImage->getImageFrame(), false);
+        if (tangentPlaneFrame.getArea() == 0)
+            tangentPlaneFrame = CTPFrame;
+        else
+            tangentPlaneFrame += CTPFrame;
+    }
+
+    // Convert tangent plane coordinates to RaDec.
+    AstrometryTransformLinear identity;
+    TanPixelToRaDec commonTangentPlaneToRaDec(identity, _commonTangentPoint);
+    Frame raDecFrame = commonTangentPlaneToRaDec.apply(tangentPlaneFrame, false);
+
+    std::vector<sphgeom::UnitVector3d> points;
+    points.reserve(4);
+    // raDecFrame is in on-sky (RA,Dec) degrees stored as an x/y box:
+    // the on-sky bounding box it represents is given by the corners of that box.
+    points.emplace_back(sphgeom::LonLat::fromDegrees(raDecFrame.xMin, raDecFrame.yMin));
+    points.emplace_back(sphgeom::LonLat::fromDegrees(raDecFrame.xMax, raDecFrame.yMin));
+    points.emplace_back(sphgeom::LonLat::fromDegrees(raDecFrame.xMin, raDecFrame.yMax));
+    points.emplace_back(sphgeom::LonLat::fromDegrees(raDecFrame.xMax, raDecFrame.yMax));
+
+    return sphgeom::ConvexPolygon::convexHull(points).getBoundingCircle();
+}
+
 void Associations::associateCatalogs(const double matchCutInArcSec, const bool useFittedList,
                                      const bool enlargeFittedList) {
     // clear reference stars
@@ -236,30 +269,6 @@ void Associations::collectRefStars(afw::table::SimpleCatalog &refCat, afw::geom:
     associateRefStars(matchCut.asArcseconds(), &raDecToCommonTangentPlane);
 }
 
-const lsst::afw::geom::Box2D Associations::getRaDecBBox() {
-    // compute the frame on the CTP that contains all input images
-    Frame tangentPlaneFrame;
-
-    for (auto const &ccdImage : ccdImageList) {
-        Frame CTPFrame = ccdImage->getPixelToCommonTangentPlane()->apply(ccdImage->getImageFrame(), false);
-        if (tangentPlaneFrame.getArea() == 0)
-            tangentPlaneFrame = CTPFrame;
-        else
-            tangentPlaneFrame += CTPFrame;
-    }
-
-    // convert tangent plane coordinates to RaDec:
-    AstrometryTransformLinear identity;
-    TanPixelToRaDec commonTangentPlaneToRaDec(identity, _commonTangentPoint);
-    Frame raDecFrame = commonTangentPlaneToRaDec.apply(tangentPlaneFrame, false);
-
-    lsst::afw::geom::Point<double> min(raDecFrame.xMin, raDecFrame.yMin);
-    lsst::afw::geom::Point<double> max(raDecFrame.xMax, raDecFrame.yMax);
-    lsst::afw::geom::Box2D box(min, max);
-
-    return box;
-}
-
 void Associations::associateRefStars(double matchCutInArcSec, const AstrometryTransform *transform) {
     // associate with FittedStars
     // 3600 because coordinates are in degrees (in CTP).

tests/test_jointcal.py

@@ -29,6 +29,7 @@
 
 import lsst.afw.table
 import lsst.daf.persistence
+from lsst.daf.base import DateTime
 import lsst.geom
 from lsst.meas.algorithms import getRefFluxField, LoadIndexedReferenceObjectsTask, DatasetConfig
 import lsst.pipe.base
@@ -269,6 +270,81 @@ def test_fit_astrometry_writeChi2(self):
             fit.return_value.saveChi2Contributions.assert_has_calls(expected)
 
 
+class TestComputeBoundingCircle(lsst.utils.tests.TestCase):
+    """Tests of Associations.computeBoundingCircle()"""
+    def _checkPointsInCircle(self, points, center, radius):
+        """Check that all points are within the (center, radius) circle.
+
+        The test is whether the max(points - center) separation is equal to
+        (or slightly less than) radius.
+        """
+        maxSeparation = 0*lsst.geom.degrees
+        for point in points:
+            maxSeparation = max(maxSeparation, center.separation(point))
+        self.assertAnglesAlmostEqual(maxSeparation, radius, maxDiff=3*lsst.geom.arcseconds)
+        self.assertLess(maxSeparation, radius)
+
+    def _testPoints(self, ccdImage1, ccdImage2, skyWcs1, skyWcs2, bbox):
+        """Fill an Associations object and test that it computes the correct
+        bounding circle for the input data.
+
+        Parameters
+        ----------
+        ccdImage1, ccdImage2 : `lsst.jointcal.CcdImage`
+            The CcdImages to add to the Associations object.
+        skyWcs1, skyWcs2 : `lsst.afw.geom.SkyWcs`
+            The WCS of each of the above images.
+        bbox : `lsst.geom.Box2D`
+            The ccd bounding box of both images.
+        """
+        lsst.log.setLevel('jointcal', lsst.log.DEBUG)
+        associations = lsst.jointcal.Associations()
+        associations.addCcdImage(ccdImage1)
+        associations.addCcdImage(ccdImage2)
+        associations.computeCommonTangentPoint()
+
+        circle = associations.computeBoundingCircle()
+        center = lsst.geom.SpherePoint(circle.getCenter())
+        radius = lsst.geom.Angle(circle.getOpeningAngle().asRadians(), lsst.geom.radians)
+        points = [lsst.geom.SpherePoint(skyWcs1.pixelToSky(lsst.geom.Point2D(x)))
+                  for x in bbox.getCorners()]
+        points.extend([lsst.geom.SpherePoint(skyWcs2.pixelToSky(lsst.geom.Point2D(x)))
+                      for x in bbox.getCorners()])
+        self._checkPointsInCircle(points, center, radius)
+
+    def testPoints(self):
+        """Test for points in an "easy" area, far from RA=0 or the poles."""
+        struct = lsst.jointcal.testUtils.createTwoFakeCcdImages()
+        self._testPoints(struct.ccdImageList[0], struct.ccdImageList[1],
+                         struct.skyWcs[0], struct.skyWcs[1], struct.bbox)
+
+    def testPointsRA0(self):
+        """Test for CcdImages crossing RA=0; this demonstrates a fix for
+        the bug described in DM-19802.
+        """
+        # Use the same pixel origins as the cfht_minimal data, but put the sky origin at RA=0
+        crpix = lsst.geom.Point2D(931.517869, 2438.572109)
+        cd = np.array([[5.19513851e-05, -2.81124812e-07],
+                      [-3.25186974e-07, -5.19112119e-05]])
+        crval1 = lsst.geom.SpherePoint(0.01, -0.01, lsst.geom.degrees)
+        crval2 = lsst.geom.SpherePoint(-0.01, 0.01, lsst.geom.degrees)
+        wcs1 = lsst.afw.geom.makeSkyWcs(crpix, crval1, cd)
+        wcs2 = lsst.afw.geom.makeSkyWcs(crpix, crval2, cd)
+
+        # Put the visit boresights at the WCS origin, for consistency
+        visitInfo1 = lsst.afw.image.VisitInfo(exposureId=30577512,
+                                              date=DateTime(65321.1),
+                                              boresightRaDec=wcs1.getSkyOrigin())
+        visitInfo2 = lsst.afw.image.VisitInfo(exposureId=30621144,
+                                              date=DateTime(65322.1),
+                                              boresightRaDec=wcs1.getSkyOrigin())
+
+        struct = lsst.jointcal.testUtils.createTwoFakeCcdImages(fakeWcses=[wcs1, wcs2],
+                                                                fakeVisitInfos=[visitInfo1, visitInfo2])
+        self._testPoints(struct.ccdImageList[0], struct.ccdImageList[1],
+                         struct.skyWcs[0], struct.skyWcs[1], struct.bbox)
+
+
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass