diff --git a/python/lsst/pipe/tasks/functors.py b/python/lsst/pipe/tasks/functors.py index d1c1cfd69..f411f546a 100644 --- a/python/lsst/pipe/tasks/functors.py +++ b/python/lsst/pipe/tasks/functors.py @@ -28,6 +28,7 @@ "ComputePixelScale", "ConvertPixelToArcseconds", "ConvertPixelSqToArcsecondsSq", "ConvertDetectorAngleToPositionAngle", + "ConvertDetectorAngleErrToPositionAngleErr," "ReferenceBand", "Photometry", "NanoJansky", "NanoJanskyErr", "LocalPhotometry", "LocalNanojansky", "LocalNanojanskyErr", "LocalDipoleMeanFlux", @@ -1394,6 +1395,67 @@ 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 + # 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. @@ -1554,6 +1616,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. @@ -1992,14 +2101,7 @@ def _func(self, df): 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,7 +2217,6 @@ def sky_uv(self, df): 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. diff --git a/tests/test_functors.py b/tests/test_functors.py index a99907954..b8e2a93ca 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -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 @@ def testConvertPixelToArcseconds(self): 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.