DM-52677: Add functors for trailLengthErr and trailAngleErr

python/lsst/pipe/tasks/functors.py

@@ -19,7 +19,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 __all__ = ["init_fromDict", "Functor", "CompositeFunctor", "mag_aware_eval",
            "CustomFunctor", "Column", "Index", "CoordColumn", "RAColumn",
            "DecColumn", "SinglePrecisionFloatColumn", "HtmIndex20", "fluxName", "fluxErrName", "Mag",
            "MagErr", "MagDiff", "Color", "DeconvolvedMoments", "SdssTraceSize",
@@ -28,6 +28,7 @@
            "ComputePixelScale", "ConvertPixelToArcseconds",
            "ConvertPixelSqToArcsecondsSq",
            "ConvertDetectorAngleToPositionAngle",
+           "ConvertDetectorAngleErrToPositionAngleErr,"
            "ReferenceBand", "Photometry",
            "NanoJansky", "NanoJanskyErr", "LocalPhotometry", "LocalNanojansky",
            "LocalNanojanskyErr", "LocalDipoleMeanFlux",
@@ -1394,8 +1395,69 @@
         # Position angle of vector from (RA1, Dec1) to (RA2, Dec2)
         return np.rad2deg(self.computePositionAngle(ra1, dec1, ra2, dec2))
 
+    def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11,
+                                                cd12, cd21, cd22):
+        """Compute position angle error (E of N) from detector angle error.
+
+        Parameters
+        ----------
+        theta : `float`
+            detector angle [radian]
+        theta_err : `float`
+            detector angle err [radian]
+        cd11 : `float`
+            [1, 1] element of the local Wcs affine transform.
+        cd12 : `float`
+            [1, 2] element of the local Wcs affine transform.
+        cd21 : `float`
+            [2, 1] element of the local Wcs affine transform.
+        cd22 : `float`
+            [2, 2] element of the local Wcs affine transform.
+        Returns
+        -------
+        Position Angle Error: `~pandas.Series`
+            Position angle error in degrees
+        """
+        # Need to compute abs(dPA/dtheta)*theta_Err to get propogated errors
+        # Unit vector in (x, y) along da
+        dx = np.cos(theta)
+        dy = np.sin(theta)
+        u = dx * cd11 + dy * cd12
+        v = dx * cd21 + dy * cd22
+
+        # Derivative of position angle wrt detector angle
+        du_dtheta = -np.sin(theta) * cd11 + np.cos(theta) * cd12
+        dv_dtheta = -np.sin(theta) * cd21 + np.cos(theta) * cd22
+
+        # Precomputing sin/cos values
+        sin_u = np.sin(u)
+        cos_u = np.cos(u)
+        sin_v = np.sin(v)
+        cos_v = np.cos(v)
+
+        # PA is atan2(Y, X), for great circle
+        pa_y = sin_u * cos_v
+        pa_x = sin_v
+
+        # We need dPA/dtheta for error propogation.
+        # X' = cos(v) * v'
+        dX_dtheta = cos_v * dv_dtheta
+        # Y' = (cos(u)*cos(v))*u' + (sin(u)*(-sin(v)))*v'
+        dY_dtheta = (cos_u * cos_v) * du_dtheta - (sin_u * sin_v) * dv_dtheta
+
+        denom = pa_x * pa_x + pa_y * pa_y
+
+        dPA_dtheta = (pa_x * dY_dtheta - pa_y * dX_dtheta) / denom
+        dPA_dtheta = np.where(denom == 0, np.nan, dPA_dtheta)
+
+        pa_err = np.rad2deg(np.abs(dPA_dtheta) * theta_err)
+
+        logging.info("PA Error: %s" % pa_err)
+        logging.info("theta_err: %s" % theta_err)
+
+        return pa_err
 
 class ComputePixelScale(LocalWcs):
     """Compute the local pixel scale from the stored CDMatrix.
     """
     name = "PixelScale"
@@ -1554,6 +1616,53 @@
         )
 
 
+class ConvertDetectorAngleErrToPositionAngleErr(LocalWcs):
+    """Compute a position angle error from a detector angle error and the
+    stored CDMatrix.
+
+    Returns
+    -------
+    position angle error : degrees
+    """
+
+    name = "PositionAngleErr"
+
+    def __init__(
+        self,
+        theta_col,
+        theta_err_col,
+        colCD_1_1,
+        colCD_1_2,
+        colCD_2_1,
+        colCD_2_2,
+        **kwargs
+    ):
+        self.theta_col = theta_col
+        self.theta_err_col = theta_err_col
+        super().__init__(colCD_1_1, colCD_1_2, colCD_2_1, colCD_2_2, **kwargs)
+
+    @property
+    def columns(self):
+        return [
+            self.theta_col,
+            self.theta_err_col,
+            self.colCD_1_1,
+            self.colCD_1_2,
+            self.colCD_2_1,
+            self.colCD_2_2
+        ]
+
+    def _func(self, df):
+        return self.getPositionAngleErrFromDetectorAngleErr(
+            df[self.theta_col],
+            df[self.theta_err_col],
+            df[self.colCD_1_1],
+            df[self.colCD_1_2],
+            df[self.colCD_2_1],
+            df[self.colCD_2_2]
+        )
+
+
 class ReferenceBand(Functor):
     """Return the band used to seed multiband forced photometry.
 
@@ -1992,14 +2101,7 @@
 
 
 class MomentsBase(Functor):
-    """Base class for functors that use shape moments and localWCS
-
-    Attributes
-    ----------
-    is_covariance : bool
-        Whether the shape columns are terms of a covariance matrix. If False,
-        they will be assumed to be terms of a correlation matrix instead.
-    """
+    """Base class for functors that use shape moments and localWCS"""
 
     is_covariance: bool = True
 
@@ -2115,8 +2217,7 @@
         CD_2_2 = df[self.colCD_2_2]
         return ((CD_1_1 * i_xx + CD_1_2 * i_xy) * CD_2_1
                 + (CD_1_1 * i_xy + CD_1_2 * i_yy) * CD_2_2)
-
     def get_g1(self, df):
         """
         Calculate shear-type ellipticity parameter G1.
         """

tests/test_functors.py

@@ -36,7 +36,7 @@
 import lsst.meas.base as measBase
 import lsst.utils.tests
 from lsst.pipe.base import InMemoryDatasetHandle
 from lsst.pipe.tasks.functors import (CompositeFunctor, CustomFunctor, Column, RAColumn,
                                       DecColumn, Mag, MagDiff, Color,
                                       DeconvolvedMoments, SdssTraceSize, PsfSdssTraceSizeDiff,
                                       E1, E2, RadiusFromQuadrupole,
@@ -53,8 +53,9 @@
                                       MomentsG1Sky, MomentsG2Sky, MomentsTraceSky,
                                       CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
                                       SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
-                                      PositionAngleFromCorrelation
-                                      )
+                                      PositionAngleFromCorrelation,
+                                      PositionAngleFromMoments,
+                                      ConvertDetectorAngleErrToPositionAngleErr)
 
 ROOT = os.path.abspath(os.path.dirname(__file__))
 
@@ -860,6 +861,109 @@
                                             atol=1e-16,
                                             rtol=1e-16))
 
+    def testConvertDetectorAngleErrToPositionAngleErr(self):
+        """Test conversion of position angle err in detector degrees to
+        position angle err on sky.
+
+        Requires a similar setup to testConvertDetectorAngleToPositionAngle.
+        """
+        dipoleSep = 10
+        ixx = 10
+        testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
+
+
+        for dec in np.linspace(-80, 80, 10):
+            for theta in (-35, 0, 90):
+                for x, y in zip(
+                        np.random.uniform(2 * 1109.99981456774, size=10),
+                        np.random.uniform(2 * 560.018167811613, size=10)):
+                    wcs = self._makeWcs(dec=dec, theta=theta)
+                    cd_matrix = wcs.getCdMatrix()
+
+                    self.dataDict["dipoleSep"] = np.full(self.nRecords,
+                                                         dipoleSep)
+                    self.dataDict["ixx"] = np.full(self.nRecords, ixx)
+                    self.dataDict["slot_Centroid_x"] = np.full(self.nRecords,
+                                                               x)
+                    self.dataDict["slot_Centroid_y"] = np.full(self.nRecords,
+                                                               y)
+                    self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
+                    self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
+                    self.dataDict["orientation"] = np.arctan2(
+                        self.dataDict["someCentroid_y"] - self.dataDict[
+                            "slot_Centroid_y"],
+                        self.dataDict["someCentroid_x"] - self.dataDict[
+                            "slot_Centroid_x"],
+                    )
+                    self.dataDict["orientation_err"] = np.arctan2(
+                        self.dataDict["someCentroid_y"] - self.dataDict[
+                            "slot_Centroid_y"],
+                        self.dataDict["someCentroid_x"] - self.dataDict[
+                            "slot_Centroid_x"],
+                    ) * 0.001
+
+                    self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(
+                        self.nRecords,
+                        cd_matrix[0, 0])
+                    self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(
+                        self.nRecords,
+                        cd_matrix[0, 1])
+                    self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(
+                        self.nRecords,
+                        cd_matrix[1, 0])
+                    self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(
+                        self.nRecords,
+                        cd_matrix[1, 1])
+                    df = self.getMultiIndexDataFrame(self.dataDict)
+
+                    # Test detector angle err to position angle err conversion
+                    func = ConvertDetectorAngleErrToPositionAngleErr(
+                        "orientation",
+                        "orientation_err",
+                        "base_LocalWcs_CDMatrix_1_1",
+                        "base_LocalWcs_CDMatrix_1_2",
+                        "base_LocalWcs_CDMatrix_2_1",
+                        "base_LocalWcs_CDMatrix_2_2"
+                    )
+                    val = self._funcVal(func, df)
+
+                    # Numerical derivative of the same PA function so a
+                    # compariosn can be made
+                    step = self.dataDict["orientation_err"] # radians. Bigger? Smaller?
+                    pa_plus_deg = func.getPositionAngleFromDetectorAngle(
+                        df[("meas", "g", "orientation")] + step,
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")],
+                    )
+                    pa_minus_deg = func.getPositionAngleFromDetectorAngle(
+                        df[("meas", "g", "orientation")] - step,
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
+                        df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")],
+                    )
+
+                    pa_plus = np.deg2rad(pa_plus_deg.to_numpy())
+                    pa_minus = np.deg2rad(pa_minus_deg.to_numpy())
+
+                    # From example: smallest signed angular difference in
+                    # (-pi, +pi]
+                    dpa = np.angle(np.exp(1j * (pa_plus - pa_minus)))
+                    dpa_dtheta = dpa / (2.0 * step)
+
+                    theta_err = df[("meas", "g", "orientation_err")].to_numpy()
+                    expected_pa_err_deg = np.rad2deg(
+                        np.abs(dpa_dtheta) * theta_err)
+
+                    np.testing.assert_allclose(
+                        val.to_numpy(),
+                        expected_pa_err_deg,
+                        rtol=0,
+                        atol=1e-6,
+                    )
+
     def _makeWcs(self, dec=53.1595451514076, theta=0):
         """Create a wcs from real CFHT values, rotated by an optional theta.