Review for DM-2674: Make meas_mosaic work on LSST.

.gitignore

@@ -4,6 +4,7 @@ config.log
 *.o
 *.os
 *.so
+*.dylib
 *.cfgc
 *.pyc
 *_wrap.cc
@@ -12,5 +13,10 @@ doc/html
 doc/*.tag
 doc/*.inc
 doc/doxygen.conf
+tests/.cache
 tests/.tests
 version.py
+bin/calibrateCatalog.py
+bin/calibrateExposure.py
+bin/checkMosaic.py
+bin/mosaic.py

bin.src/SConscript

@@ -0,0 +1,3 @@
+# -*- python -*-
+from lsst.sconsUtils import scripts
+scripts.BasicSConscript.shebang()

include/lsst/meas/mosaic/fluxfit.h

@@ -11,7 +11,7 @@ namespace lsst {
 
 	    class FluxFitParams {
 	    public:
-		typedef boost::shared_ptr<FluxFitParams> Ptr;
+		typedef std::shared_ptr<FluxFitParams> Ptr;
 
 		int order;
 		bool chebyshev;
@@ -58,23 +58,23 @@ namespace lsst {
 
 	    FluxFitParams::Ptr
 	      convertFluxFitParams(FluxFitParams::Ptr& ffp,
-				   lsst::afw::cameraGeom::Ccd::Ptr& ccd,
+				   PTR(lsst::afw::cameraGeom::Detector)& ccd,
 				   double x0=0.0, double y0=0.0);
 
 	    lsst::daf::base::PropertySet::Ptr
 	      metadataFromFluxFitParams(FluxFitParams::Ptr& ffp);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getFCorImg(FluxFitParams::Ptr& p,
-			 lsst::afw::cameraGeom::Ccd::Ptr& ccd,
+			 PTR(lsst::afw::cameraGeom::Detector)& ccd,
 			 Coeff::Ptr& coeff);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getFCorImg(FluxFitParams::Ptr& p, int width, int height);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getFCorImg(FluxFitParams::Ptr& p,
-			 lsst::afw::cameraGeom::Ccd::Ptr& ccd);
+			 PTR(lsst::afw::cameraGeom::Detector)& ccd);
 
 #include "chebyshev.h"
     }

include/lsst/meas/mosaic/mosaicfit.h

@@ -2,14 +2,14 @@
 #if !defined(HSC_MEAS_MOSAIC_H)
 #define HSC_MEAS_MOSAIC_H
 
+#include <cmath>
+#include <memory>
 #include <vector>
 #include "lsst/pex/exceptions.h"
 #include "lsst/afw/image.h"
 #include "lsst/afw/geom.h"
 #include "lsst/afw/cameraGeom.h"
 #include "lsst/afw/table.h"
-#include "lsst/utils/ieee.h"
-#include "boost/enable_shared_from_this.hpp"
 
 namespace lsst {
     namespace meas {
@@ -33,11 +33,11 @@ namespace lsst {
                     _err(std::numeric_limits<double>::quiet_NaN()),
                     _xerr(std::numeric_limits<double>::quiet_NaN()),
                     _yerr(std::numeric_limits<double>::quiet_NaN()),
-                    _astromBad(!lsst::utils::isfinite(_flux))
+                    _astromBad(!std::isfinite(_flux))
                     {
                         try {
-                            _err = record.get(record.getSchema().find<double>("flux.err").key);
-                        } catch (pex::exceptions::NotFoundException const&) {
+                            _err = record.get(record.getSchema().find<double>("fluxSigma").key);
+                        } catch (pex::exceptions::NotFoundError const&) {
 			    // flux.err is not availabe. sqrt of flux will be used
 			    // 1.E-08 is a scaling factor to make mag_err ~ 0.1mag for 15 mag
 			    // This is a good approximation for 2MASS J-band
@@ -111,7 +111,7 @@ namespace lsst {
 
 	    class Poly {
 	    public:
-	        typedef boost::shared_ptr<Poly> Ptr;
+	        typedef std::shared_ptr<Poly> Ptr;
 
 		int order;
 		int ncoeff;
@@ -128,7 +128,7 @@ namespace lsst {
 
 	    class Coeff {
 	    public:
-		typedef boost::shared_ptr<Coeff> Ptr;
+		typedef std::shared_ptr<Coeff> Ptr;
 
 		Poly::Ptr p;
 		int iexp;
@@ -181,7 +181,7 @@ namespace lsst {
 
 	    class Obs {
 	    public:
-		typedef boost::shared_ptr<Obs> Ptr;
+		typedef std::shared_ptr<Obs> Ptr;
 
 		double ra, dec;
 		double xi, eta;
@@ -217,20 +217,20 @@ namespace lsst {
 
 		Obs(int id, double ra, double dec, double x, double y, int ichip, int iexp);
 		Obs(int id, double ra, double dec, int ichip, int iexp);
-		void setUV(lsst::afw::cameraGeom::Ccd::Ptr const &ccd, double x0=0, double y0=0);
+		void setUV(PTR(lsst::afw::cameraGeom::Detector) &ccd, double x0=0, double y0=0);
 		void setXiEta(double ra_c, double dec_c);
 		void setFitVal(Coeff::Ptr& c, Poly::Ptr p);
 		void setFitVal2(Coeff::Ptr& c, Poly::Ptr p);
 	    };
 
 	    class KDTree
 #if !defined(SWIG)
-	      : public boost::enable_shared_from_this<KDTree>
+	      : public std::enable_shared_from_this<KDTree>
 #endif
 	    {
 	    public:
-		typedef boost::shared_ptr<KDTree> Ptr;
-		typedef boost::shared_ptr<const KDTree> ConstPtr;
+		typedef std::shared_ptr<KDTree> Ptr;
+		typedef std::shared_ptr<const KDTree> ConstPtr;
 
 		int depth;
 		int axis;
@@ -269,7 +269,7 @@ namespace lsst {
                 void _initializeMatches(SourceMatchSet& m, int depth);
             };
 
-	    typedef std::map<int, lsst::afw::cameraGeom::Ccd::Ptr> CcdSet;
+	    typedef std::map<int, PTR(lsst::afw::cameraGeom::Detector)> CcdSet;
 	    typedef std::map<int, Coeff::Ptr> CoeffSet;
 	    typedef std::vector<Obs::Ptr> ObsVec;
 
@@ -314,23 +314,23 @@ namespace lsst {
                                      std::string const & snapshotDir = ".");
 
 	    Coeff::Ptr convertCoeff(Coeff::Ptr& coeff,
-				    lsst::afw::cameraGeom::Ccd::Ptr& ccd);
+				    PTR(lsst::afw::cameraGeom::Detector)& ccd);
 
 	    lsst::afw::image::TanWcs::Ptr wcsFromCoeff(Coeff::Ptr& coeff);
 
             Coeff::Ptr coeffFromTanWcs(lsst::afw::image::Wcs::Ptr& wcs);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getJImg(Coeff::Ptr& coeff,
-		      lsst::afw::cameraGeom::Ccd::Ptr& ccd);
+		      PTR(lsst::afw::cameraGeom::Detector)& ccd);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getJImg(lsst::afw::image::Wcs::Ptr& wcs,
 		      int width, int height);
 
 	    lsst::afw::image::Image<float>::Ptr
 	      getJImg(lsst::afw::image::Wcs::Ptr& wcs,
-		      lsst::afw::cameraGeom::Ccd::Ptr& ccd);
+		      PTR(lsst::afw::cameraGeom::Detector)& ccd);
 
             //{
             // Calculate Jacobian correction
@@ -353,7 +353,6 @@ namespace lsst {
                 ndarray::Array<double const, 1> const& y  ///< y positions for correction
                 );
             //}
-
     }
   }
 }

include/lsst/meas/mosaic/shimCameraGeom.h

@@ -0,0 +1,50 @@
+// -*- lsst-c++ -*-
+#ifndef LSST_MEAS_MOSAIC_SHIM_CAMERAGEOM_H
+#define LSST_MEAS_MOSAIC_SHIM_CAMERAGEOM_H
+
+#include "lsst/afw/geom.h"
+#include "lsst/afw/cameraGeom.h"
+
+namespace lsst {
+namespace meas {
+namespace mosaic {
+
+// Get the number of quarter rotations of the detector.
+int getNQuarter(CONST_PTR(afw::cameraGeom::Detector));
+
+// Get the detector yaw.
+afw::geom::Angle getYaw(CONST_PTR(afw::cameraGeom::Detector));
+
+// Return a linear transform which scales from dimensions in mm to dimensions
+// in pixels.
+afw::geom::LinearTransform makeScalingMmToPx(afw::geom::Extent2D const pSize);
+
+// Return the position of the center of the detector in pixels on the focal plane.
+// Mimics HSC's camGeom: ccd.getCenter().getPixels(ccd.getPixelSize())
+afw::geom::Point2D getCenterInFpPixels(CONST_PTR(afw::cameraGeom::Detector));
+
+// Return the position of the center of the detector in pixels on the
+// detector.
+// Mimics HSC's camGeom: ccd.getCenterPixel()
+afw::geom::Point2D getCenterInDetectorPixels(CONST_PTR(afw::cameraGeom::Detector));
+
+// Return the width of the detector in pixels.
+int getWidth(CONST_PTR(afw::cameraGeom::Detector));
+
+// Return the height of the detector in pixels.
+int getHeight(CONST_PTR(afw::cameraGeom::Detector));
+
+// Convert a pixel position on a given detector to a pixel position on the focal plane.
+afw::geom::Point2D detPxToFpPx(CONST_PTR(afw::cameraGeom::Detector), afw::geom::Point2D const);
+
+// Convert a pixel position on a given detector to a pixel position on the focal plane
+// accounting for yaw rotation.
+// Mimics HSC's camGeom: ccd.getPositionFromPixel(point).getPixels(ccd.getPixelSize())
+afw::geom::Point2D detPxToFpPxRot(CONST_PTR(afw::cameraGeom::Detector), afw::geom::Point2D const);
+
+// Compute new position of lower left corner in Focal Plane pixels: X0, Y0
+afw::geom::Point2D computeX0Y0(CONST_PTR(afw::cameraGeom::Detector), double, double);
+
+}}} // lsst::meas::mosaic
+
+#endif

include/lsst/meas/mosaic/spatialCellSource.h

@@ -9,8 +9,8 @@ namespace lsst {
 
       class SpatialCellSource : public lsst::afw::math::SpatialCellCandidate {
       public:
-	typedef boost::shared_ptr<SpatialCellSource> Ptr;
-	typedef boost::shared_ptr<const SpatialCellSource> ConstPtr;
+	typedef std::shared_ptr<SpatialCellSource> Ptr;
+	typedef std::shared_ptr<const SpatialCellSource> ConstPtr;
 
 	SpatialCellSource(float const xCenter,
 			  float const yCenter,

python/lsst/meas/mosaic/calibrate.py

@@ -1,10 +1,17 @@
 from lsst.pipe.base import CmdLineTask, ArgumentParser
-from lsst.pex.config import Config, Field
+from lsst.pex.config import Config, Field, DictField
 from .updateExposure import applyMosaicResultsExposure, applyMosaicResultsCatalog, applyCalib
-from lsst.pipe.tasks.dataIds import PerTractCcdDataIdContainer
+from lsst.meas.base.forcedPhotCcd import PerTractCcdDataIdContainer
+from . import utils as mosaicUtils
 
 class CalibrateCatalogConfig(Config):
     doApplyCalib = Field(dtype=bool, default=True, doc="Calibrate fluxes to magnitudes?")
+    srcSchemaMap = DictField(
+        doc="Mapping between different stack (e.g. HSC vs. LSST) schema names",
+        keytype=str,
+        itemtype=str,
+        default=None,
+        optional=True)
 
 class CalibrateCatalogTask(CmdLineTask):
     ConfigClass = CalibrateCatalogConfig
@@ -19,10 +26,18 @@ def _makeArgumentParser(cls):
 
     def run(self, dataRef):
         catalog = dataRef.get("src", immediate=True)
+        calexp_md = dataRef.get('calexp_md', immediate=True)
+        # Check if we are looking at HSC stack outputs
+        hscRun = mosaicUtils.checkHscStack(calexp_md)
+        # Set the aliap map for the source catalog
+        if self.config.srcSchemaMap is not None and hscRun is not None:
+            aliasMap = catalog.schema.getAliasMap()
+            for lsstName, otherName in self.config.srcSchemaMap.iteritems():
+                aliasMap.set(lsstName, otherName)
         results = applyMosaicResultsCatalog(dataRef, catalog)
         catalog = results.catalog
         if self.config.doApplyCalib:
-            catalog = applyCalib(catalog, results.ffp.calib)
+            catalog = applyCalib(catalog, results.ffp.calib, hscRun=hscRun)
         dataRef.put(catalog, "calibrated_src")
 
     def writeConfig(self, *args, **kwargs):

python/lsst/meas/mosaic/checkMosaicTask.py

@@ -13,12 +13,13 @@
 from lsst.meas.mosaic.mosaicTask import MosaicConfig
 from lsst.meas.mosaic.mosaicTask import MosaicRunner
 import lsst.meas.mosaic.mosaicLib       as measMosaic
-import lsst.meas.astrom.astrom          as measAstrom
+import lsst.meas.algorithms             as measAlg
+import lsst.meas.astrom                 as measAstrom
 import lsst.afw.image                   as afwImage
 import lsst.afw.geom                    as afwGeom
 import lsst.afw.table                   as afwTable
 
-from lsst.meas.photocal.colorterms import ColortermLibraryConfig
+from lsst.pipe.tasks.colorterms import ColortermLibrary
 
 class CheckMosaicConfig(MosaicConfig):
     maxMag = pexConfig.Field(
@@ -424,7 +425,7 @@ def check(self, dataRefList, ct=None, debug=False, verbose=False):
 
         sourceSet = measMosaic.SourceGroup()
         matchList = measMosaic.SourceMatchGroup()
-        astrom = measAstrom.Astrometry(self.config.astrom)
+        astrom = measAstrom.ANetBasicAstrometryTask(self.config.astrom)
         ssVisit = dict()
         mlVisit = dict()
         dataRefListUsed = list()
@@ -467,25 +468,41 @@ def check(self, dataRefList, ct=None, debug=False, verbose=False):
                               flags=afwTable.SOURCE_IO_NO_FOOTPRINTS,
                               immediate=True)
                 packedMatches = dataRef.get('icMatch')
-                matches = astrom.joinMatchListWithCatalog(packedMatches, icSrces, True)
+                matches = astrom.joinMatchListWithCatalog(packedMatches, icSrces)
 
                 matches = [m for m in matches if m.first != None]
-                if ct != None and len(matches) != 0:
+
+                if matches:
                     refSchema = matches[0].first.schema
-                    key_p = refSchema.find(ct.primary).key
-                    key_s = refSchema.find(ct.secondary).key
-                    key_f = refSchema.find("flux").key
-                    for m in matches:
-                        refFlux1 = m.first.get(key_p)
-                        refFlux2 = m.first.get(key_s)
-                        refMag1 = -2.5*math.log10(refFlux1)
-                        refMag2 = -2.5*math.log10(refFlux2)
+                    if ct:
+                        # Add a "flux" field to the match records which contains the
+                        # colorterm-corrected reference flux. The field name is hard-coded in
+                        # lsst::meas::mosaic::Source.
+                        mapper = afwTable.SchemaMapper(refSchema)
+                        for key, field in refSchema:
+                            mapper.addMapping(key)
+                        key_f = mapper.editOutputSchema().addField("flux", type=float, doc="Reference flux")
+                        table = afwTable.SimpleTable.make(mapper.getOutputSchema())
+                        table.preallocate(len(matches))
+                        for match in matches:
+                            newMatch = table.makeRecord()
+                            newMatch.assign(match.first, mapper)
+                            match.first = newMatch
+
+                        key_p = refSchema.find(refSchema.join(ct.primary, "flux")).key
+                        key_s = refSchema.find(refSchema.join(ct.secondary, "flux")).key
+                        refFlux1 = numpy.array([m.first.get(key_p) for m in matches])
+                        refFlux2 = numpy.array([m.first.get(key_s) for m in matches])
+                        refMag1 = -2.5*numpy.log10(refFlux1)
+                        refMag2 = -2.5*numpy.log10(refFlux2)
                         refMag = ct.transformMags(refMag1, refMag2)
-                        refFlux = math.pow(10.0, -0.4*refMag)
-                        if refFlux == refFlux:
-                            m.first.set(key_f, refFlux)
-                        else:
-                            m.first = None
+                        refFlux = numpy.power(10.0, -0.4*refMag)
+                        matches = [setCatFlux(m, f, key_f) for m, f in zip(matches, refFlux) if f == f]
+                    else:
+                        # No colorterm; we can get away with aliasing the reference flux.
+                        refFluxField = measAlg.getRefFluxField(refSchema, filterName)
+                        refSchema.getAliasMap().set("flux", refFluxField)
+
                 sources = self.selectStars(sources)
                 matches = self.selectStars(matches, True)
             except Exception, e:
@@ -528,9 +545,9 @@ def check(self, dataRefList, ct=None, debug=False, verbose=False):
         self.plotPosAsMag(m0_s, dx_s, dy_s)
         self.writeCatalog(allSource, wcsDic, calibDic, ffpDic)
 
-    def run(self, camera, butler, dataRefList, debug):
+    def run(self, camera, butler, tract, dataRefList, debug):
 
-        colorterms = ColortermLibraryConfig()
+        colorterms = ColortermLibrary()
         name = os.path.join(os.environ["OBS_SUBARU_DIR"], "config", camera, "colorterms.py")
         if os.path.exists(name):
             colorterms.load(name)
@@ -544,7 +561,7 @@ def run(self, camera, butler, dataRefList, debug):
             self.log.fatal("There are %d filters in input frames" % len(filters))
             return None
 
-        ct = colorterms.selectColorTerm(butler.mapper.filters[filters[0]])
+        ct = colorterms.getColorterm(butler.mapper.filters[filters[0]], self.config.photoCatName)
 
         return self.check(dataRefList, ct, debug)