Add ConvertDetectorAngleErrToPositionAngleErr Functor

Test position angle error prop

Fix

Fix functor

Test new functor
Without debuging

Logging for testing

Author: bsmartradio

python/lsst/pipe/tasks/functors.py

@@ -1394,6 +1394,55 @@ def getPositionAngleFromDetectorAngle(self, theta, cd11, cd12, cd21, cd22):
         # 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
+
+        # Get unit direction
+        dx = np.cos(theta)
+        dy = np.sin(theta)
+
+        # Transform it using WCS?
+        u = dx * cd11 + dy * cd12
+        v = dx * cd21 + dy * cd22
+        # Now we are computing the tangent
+        ratio = u / v
+
+        # Get derivative of theta
+        du_dtheta = -np.sin(theta) * cd11 + np.cos(theta) * cd12
+        dv_dtheta = -np.sin(theta) * cd21 + np.cos(theta) * cd22
+
+        # Get derivative of tangent
+        d_ratio_dtheta = (v * du_dtheta - u * dv_dtheta) / v ** 2
+        dPA_dtheta = (1 / (1 + ratio ** 2)) * d_ratio_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.
@@ -1554,6 +1603,53 @@ def _func(self, df):
         )
 
 
+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.
 

tests/test_functors.py

@@ -49,11 +49,10 @@
                                       ConvertDetectorAngleToPositionAngle,
                                       HtmIndex20, Ebv, MomentsIuuSky, MomentsIvvSky, MomentsIuvSky,
                                       SemimajorAxisFromMoments, SemiminorAxisFromMoments,
-                                      PositionAngleFromMoments,
                                       CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
                                       SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
                                       PositionAngleFromCorrelation
-                                      )
+                                      PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr)
 
 ROOT = os.path.abspath(os.path.dirname(__file__))
 
@@ -762,6 +761,70 @@ def testConvertDetectorAngleToPositionAngle(self):
 
                     np.testing.assert_allclose(coord_diff, 0, rtol=0, atol=5e-6)
 
+    # Test position angle error propogation
+    def testConvertDetectorAngleErrToPositionAngleErr(self):
+        """Test conversion of position angle err in detector degrees to
+        position angle erron sky
+        """
+        dipoleSep = 10
+        ixx = 10
+        testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
+        # testConvertPixelToArcSecond uses a meas_base LocalWcsPlugin
+        # but we're using a simple WCS as our example, so this doesn't really matter
+        # and we'll just use the WCS directly
+        for dec in np.linspace(-90, 90, 10):
+            for theta in (-45, 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"],
+                    )*.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 to position angle 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"
+                    )
+
+                    func_pa = ConvertDetectorAngleToPositionAngle(
+                        "orientation",
+                        "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)
+                    val_pa = self._funcVal(func_pa, df)
+
+
     def testConvertPixelToArcseconds(self):
         """Test calculations of the pixel scale, conversions of pixel to
         arcseconds.