Merge pull request #183 from lsst/tickets/DM-24278

DM-24278: Apply proper motion to matched sources in Jointcal

include/lsst/jointcal/Associations.h

@@ -78,12 +78,14 @@ class Associations {
      * parallelize the creation of the ccdImages.
      *
      * @param imageList A pre-built ccdImage list.
+     * @param epoch The julian epoch year to which all proper motion corrections should be made.
      */
-    Associations(CcdImageList const &imageList)
+    Associations(CcdImageList const &imageList, double epoch = 0)
             : ccdImageList(imageList),
               _commonTangentPoint(Point(std::numeric_limits<double>::quiet_NaN(),
                                         std::numeric_limits<double>::quiet_NaN())),
-              _maxMeasuredStars(0) {}
+              _maxMeasuredStars(0),
+              _epoch(epoch) {}
 
     /// No moves or copies: jointcal only ever needs one Associations object.
     Associations(Associations const &) = delete;
@@ -97,17 +99,25 @@ class Associations {
     void computeCommonTangentPoint();
 
     /**
-     * @brief      Sets a shared tangent point for all ccdImages.
+     * Shared tangent point for all ccdImages (decimal degrees).
      *
-     * @param      commonTangentPoint  The common tangent point of all input images (decimal degrees).
+     * Used to project sidereal coordinates related to the images onto a single plane.
      */
+    ///@{
     void setCommonTangentPoint(lsst::geom::Point2D const &commonTangentPoint);
-
-    //! can be used to project sidereal coordinates related to the image set on a plane.
     Point getCommonTangentPoint() const { return _commonTangentPoint; }
+    ///@}
 
     size_t getMaxMeasuredStars() const { return _maxMeasuredStars; }
 
+    /**
+     * Common epoch of all of the ccdImages as a Julian Epoch Year (e.g. 2000.0 for J2000).
+     */
+    ///@{
+    double getEpoch() const { return _epoch; }
+    void setEpoch(double epoch) { _epoch = epoch; }
+    ///@}
+
     /**
      * @brief      Create a ccdImage from an exposure catalog and metadata, and add it to the list.
      *
@@ -157,13 +167,6 @@ class Associations {
     //! track of that.
     void deprojectFittedStars();
 
-//! Set the color field of FittedStar 's from a colored catalog.
-/* If Color is "g-i", then the color is assigned from columns "g" and "i" of the colored catalog. */
-#ifdef TODO
-    void setFittedStarColors(std::string const &dicStarListName, std::string const &color,
-                             double matchCutArcSec);
-#endif
-
     /**
      * Prepare the fittedStar list by making quality cuts and normalizing measurements.
      *
@@ -224,12 +227,18 @@ class Associations {
      */
     void normalizeFittedStars();
 
+    // Common tangent point on-sky of all of the ccdImages, typically determined by computeCommonTangentPoint.
+    // (decimal degrees)
     Point _commonTangentPoint;
 
     // The number of MeasuredStars at the start of fitting, before any outliers are removed.
     // This is used to reserve space in vectors for e.g. outlier removal, but is not updated during outlier
     // removal or cleanup, so should only be used as an upper bound on the number of MeasuredStars.
     size_t _maxMeasuredStars;
+
+    // Julian Epoch Year (e.g. 2000.0 for J2000)
+    // Common epoch of all of the ccdImages, typically computed externally via astropy and then set.
+    double _epoch;
 };
 
 }  // namespace jointcal

include/lsst/jointcal/AstrometryFit.h

@@ -91,7 +91,7 @@ class AstrometryFit : public FitterBase {
      * Set parameters to fit and assign indices in the big matrix.
      *
      * @param[in]  whatToFit   Valid strings: zero or more of "Distortions", "Positions",
-     *                         "Refrac", "PM" which define which parameter set
+     *                         "PM" which define which parameter set
      *                         is going to be variable when computing
      *                         derivatives (leastSquareDerivatives) and minimizing
      *                         (minimize()). whatToFit="Positions Distortions"
@@ -134,16 +134,10 @@ class AstrometryFit : public FitterBase {
     void saveChi2RefContributions(std::string const &filename) const override;
 
 private:
-    bool _fittingDistortions, _fittingPos, _fittingRefrac, _fittingPM;
+    bool _fittingDistortions, _fittingPos, _fittingPM;
     std::shared_ptr<AstrometryModel> _astrometryModel;
-    double _referenceColor, _sigCol;  // average and r.m.s color
-    double _refractionCoefficient;    // fit parameter
-    Eigen::Index _refracPosInMatrix;  // where it stands
-    double _JDRef;                    // average Julian date
-
-    // counts in parameter subsets.
-    std::size_t _nParRefrac;
 
+    double _epoch;     // epoch to correct proper motion/parallax to (Julian Epoch year, e.g. J2000.0)
     double _posError;  // constant term on error on position (in pixel unit)
 
     void leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList,
@@ -159,8 +153,16 @@ class AstrometryFit : public FitterBase {
 
     void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const override;
 
+    /**
+     * Transform the positions of a FittedStar into the frame of a MeasuredStar.
+     *
+     * @param fittedStar The star to transform.
+     * @param sky2TP Transformation from sky coordinates to CcdImage tangent plane.
+     * @param deltaYears Difference in years between FittedStar and MeasuredStar epochs.
+     * @return Corrected position of FittedStar.
+     */
     Point transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const &sky2TP,
-                              Point const &refractionVector, double refractionCoeff, double mjd) const;
+                              double deltaYears) const;
 
     /// Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) from one CcdImage.
     void accumulateStatImage(CcdImage const &ccdImage, Chi2Accumulator &accum) const;

include/lsst/jointcal/BaseStar.h

@@ -26,6 +26,7 @@
 #define LSST_JOINTCAL_BASE_STAR_H
 
 #include <iostream>
+#include <iomanip>
 #include <cstdio>
 #include <string>
 #include <sstream>

include/lsst/jointcal/CcdImage.h

@@ -154,8 +154,7 @@ class CcdImage {
     //!  Airmass
     double getAirMass() const { return _airMass; }
 
-    //! Julian Date
-    double getMjd() const { return _mjd; }
+    double getEpoch() const { return _epoch; }
 
     //! Return the exposure's photometric calibration
     std::shared_ptr<afw::image::PhotoCalib> getPhotoCalib() const { return _photoCalib; }
@@ -216,7 +215,7 @@ class CcdImage {
 
     lsst::geom::SpherePoint _boresightRaDec;
     double _airMass;  // airmass value.
-    double _mjd;      // modified julian date
+    double _epoch;    // julian epoch year (e.g. 2000.0 for J2000)
     std::shared_ptr<afw::image::PhotoCalib> _photoCalib;
     std::shared_ptr<afw::cameraGeom::Detector> _detector;
     // refraction

include/lsst/jointcal/FittedStar.h

@@ -29,6 +29,7 @@
 #include <fstream>
 
 #include "Eigen/Core"
+#include "lsst/jointcal/ProperMotion.h"
 #include "lsst/jointcal/BaseStar.h"
 #include "lsst/jointcal/StarList.h"
 
@@ -39,26 +40,14 @@ class MeasuredStar;
 class RefStar;
 class AstrometryTransform;
 
-/*! \file */
-
-//! objects whose position is going to be fitted. Coordinates in Common Tangent Plane.
-
-struct PmBlock {
-    // proper motion in x and y. Units depend on how you fit them
-    double pmx, pmy;
-    double epmx, epmy, epmxy;
-    double color;  // OK it is unrelated, but associated in practice
-    bool mightMove;
-
-    PmBlock() : pmx(0), pmy(0), epmx(0), epmy(0), epmxy(0), color(0), mightMove(false){};
-};
-
 /**
- * The objects which have been measured several times.
+ * FittedStars are objects whose position or flux is going to be fitted, and which come from the association
+ * of multiple MeasuredStars.
  *
+ * x/y Coordinates are in the Common Tangent Plane (degrees).
  * MeasuredStars from different CcdImages that represent the same on-sky object all point to one FittedStar.
  */
-class FittedStar : public BaseStar, public PmBlock {
+class FittedStar : public BaseStar {
 public:
     FittedStar() : BaseStar(), _indexInMatrix(-1), _measurementCount(0), _refStar(nullptr) {}
 

include/lsst/jointcal/MeasuredStar.h

@@ -41,8 +41,13 @@ double instMagFromInstFlux(double instFlux) { return -2.5 * std::log10(instFlux)
 
 class CcdImage;
 
-//! objects measured on actual images. Coordinates and uncertainties are expressed in pixel image frame. Flux
-//! expressed in ADU/s.
+/**
+ * Sources measured on images.
+ *
+ * x/y positions and uncertainties (from BaseStar) are in pixels in the image frame.
+ * instFlux[Err] are in counts.
+ * flux[Err] (from BaseStar) are nJy.
+ */
 class MeasuredStar : public BaseStar {
 public:
     MeasuredStar()

include/lsst/jointcal/PhotometryFit.h

@@ -109,11 +109,6 @@ class PhotometryFit : public FitterBase {
     /// Compute the derivatives of the reference terms
     void leastSquareDerivativesReference(FittedStarList const &fittedStarList, TripletList &tripletList,
                                          Eigen::VectorXd &grad) const override;
-
-#ifdef STORAGE
-    Point transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const *sky2TP,
-                              Point const &refractionVector, double refractionCoeff, double mjd) const;
-#endif
 };
 }  // namespace jointcal
 }  // namespace lsst

include/lsst/jointcal/ProperMotion.h

@@ -0,0 +1,83 @@
+// -*- LSST-C++ -*-
+/*
+ * This file is part of jointcal.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * 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 GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef LSST_JOINTCAL_PROPER_MOTION_H
+#define LSST_JOINTCAL_PROPER_MOTION_H
+
+#include <iostream>
+#include <memory>
+#include <math.h>  // atan2
+
+namespace lsst {
+namespace jointcal {
+
+class Point;
+
+/**
+ * Proper motion data for a reference star or fitted star.
+ *
+ * Whether to just use these values or fit them is determined by the RefStar and FittedStar they belong to.
+ *
+ * Units are radians/year
+ * Note: RA proper motion is `pm_ra*cos(dec)`
+ */
+class ProperMotion {
+public:
+    ProperMotion(double ra, double dec, double raErr, double decErr, double raDecCov = 0)
+            : _ra(ra), _dec(dec), _raErr(raErr), _decErr(decErr), _raDecCov(raDecCov) {
+        _offsetBearing = atan2(dec, ra);
+    };
+
+    ProperMotion(ProperMotion const &) = default;
+    ProperMotion(ProperMotion &&) = default;
+    ProperMotion &operator=(ProperMotion const &) = default;
+    ProperMotion &operator=(ProperMotion &&) = default;
+    ~ProperMotion() = default;
+
+    /**
+     * Apply proper motion correction to the input star, returning a star with PM-corrected coordinates.
+     *
+     * NOTE: This method does not apply to the coordinate errors (Point does not include uncertainty).
+     *
+     * @param star The star to correct for this proper motion.
+     * @param timeDeltaYears The difference in time from the correction epoch to correct for, in years.
+     *
+     * @return The star with corrected coordinates.
+     */
+    Point apply(Point star, double timeDeltaYears) const;
+
+    friend std::ostream &operator<<(std::ostream &stream, ProperMotion const &pm);
+
+private:
+    // Proper motion in ra and dec. Units are radians/year
+    double _ra, _dec;
+    double _raErr, _decErr, _raDecCov;
+
+    // Cache the bearing along-which to offset.
+    double _offsetBearing;
+};
+}  // namespace jointcal
+}  // namespace lsst
+
+#endif  // LSST_JOINTCAL_PROPER_MOTION_H

include/lsst/jointcal/RefStar.h

@@ -29,28 +29,60 @@
 #include <fstream>
 
 #include "lsst/jointcal/FittedStar.h"
+#include "lsst/jointcal/ProperMotion.h"
 #include "lsst/jointcal/StarList.h"
 
 namespace lsst {
 namespace jointcal {
 
-//! Objects used as position anchors, typically USNO stars. Coordinate system defined by user. The Common
-//! Tangent Plane seems a good idea.
+/**
+ * Objects used as position/flux anchors (e.g. Gaia DR2 stars). Coordinate system should match that of the
+ * fittedStars these are associated with; typically the common tangent plane.
+ *
+ * RefStars should have their proper motion and parallax corrections pre-applied, so that they are at
+ * the same epoch as is stored in Associations.
+ */
 class RefStar : public BaseStar {
 public:
     RefStar(double xx, double yy, double flux, double fluxErr) : BaseStar(xx, yy, flux, fluxErr) {}
 
-    /// No move or copy: each RefStar is unique, and should be accessed/managed via shared_ptr.
+    /// No copy: each RefStar is unique, and should be accessed/managed via shared_ptr.
     RefStar(RefStar const&) = delete;
-    RefStar(RefStar&&) = delete;
-    RefStar& operator=(RefStar const&) = default;
-    RefStar& operator=(RefStar&&) = delete;
+    RefStar(RefStar&&) = default;
+    RefStar& operator=(RefStar const&) = delete;
+    RefStar& operator=(RefStar&&) = default;
+
+    // pybind11 cannot handle unique_ptr as arguments, so provide this for python-level testing.
+    void setProperMotion(ProperMotion const& properMotion) {
+        _properMotion = std::make_unique<ProperMotion>(properMotion);
+    }
+    void setProperMotion(std::unique_ptr<ProperMotion const>&& properMotion) {
+        _properMotion = std::move(properMotion);
+    }
+
+    /**
+     * Apply proper motion correction to the input star, returning a star with PM-corrected
+     * coordinates and coordinate errors.
+     *
+     * Star is returned unchanged if no proper motion data is available.
+     *
+     * @param star The star to correct for this proper motion.
+     * @param timeDeltaYears The difference in time from the correction epoch to correct for, in
+     * years.
+     *
+     * @return The star with corrected coordinates.
+     */
+    Point applyProperMotion(Point star, double timeDeltaYears) const;
+
+private:
+    // RefStars are already PM corrected to a common epoch: this is to correct the associated FittedStar
+    // to each MeasuredStar's epoch. Not all refcats have PM data: this will be nullptr if no PM data is
+    // available for any reason.
+    std::unique_ptr<ProperMotion const> _properMotion;
 };
 
 /****** RefStarList ***********/
 
-class Frame;
-
 // typedef StarList<RefStar> RefStarList;
 class RefStarList : public StarList<RefStar> {};
 

python/lsst/jointcal/associations.cc

@@ -39,7 +39,7 @@ namespace {
 void declareAssociations(py::module &mod) {
     py::class_<Associations, std::shared_ptr<Associations>> cls(mod, "Associations");
     cls.def(py::init<>());
-    cls.def(py::init<CcdImageList const &>(), "imageList"_a);
+    cls.def(py::init<CcdImageList const &, double>(), "imageList"_a, "epoch"_a = 0);
 
     // NOTE: these could go away if the lists they wrap can be accessed directly.
     cls.def("refStarListSize", &Associations::refStarListSize);
@@ -66,6 +66,9 @@ void declareAssociations(py::module &mod) {
     cls.def("getCommonTangentPoint", &Associations::getCommonTangentPoint);
     cls.def("setCommonTangentPoint", &Associations::setCommonTangentPoint);
     cls.def("computeCommonTangentPoint", &Associations::computeCommonTangentPoint);
+
+    cls.def("getEpoch", &Associations::getEpoch);
+    cls.def("setEpoch", &Associations::setEpoch);
 }
 
 PYBIND11_MODULE(associations, mod) {