From 7749775a93f7cf479062029ce0e12f43b3d3c9b8 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Fri, 10 Oct 2025 14:07:18 -0700 Subject: [PATCH 1/6] Add ConvertDetectorAngleErrToPositionAngleErr Functor Test position angle error prop Fix Fix functor Test new functor Without debuging Logging for testing --- python/lsst/pipe/tasks/functors.py | 96 ++++++++++++++++++++++++++++++ tests/test_functors.py | 66 +++++++++++++++++++- 2 files changed, 161 insertions(+), 1 deletion(-) diff --git a/python/lsst/pipe/tasks/functors.py b/python/lsst/pipe/tasks/functors.py index 619c80f43..bd81a3a9c 100644 --- a/python/lsst/pipe/tasks/functors.py +++ b/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. diff --git a/tests/test_functors.py b/tests/test_functors.py index 2abd6d906..24441188d 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -54,7 +54,7 @@ CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky, SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation, PositionAngleFromCorrelation - ) + PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr) ROOT = os.path.abspath(os.path.dirname(__file__)) @@ -763,6 +763,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. From 0256eaedd75b9e25fb5ad45b2f47ad1fefccc468 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Wed, 4 Feb 2026 09:27:13 -0800 Subject: [PATCH 2/6] Test Fix Test Fixes --- python/lsst/pipe/tasks/functors.py | 47 ++++---- tests/test_functors.py | 174 ++++++++++++++++++----------- 2 files changed, 134 insertions(+), 87 deletions(-) diff --git a/python/lsst/pipe/tasks/functors.py b/python/lsst/pipe/tasks/functors.py index bd81a3a9c..bc360d4a3 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,7 +1395,8 @@ 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): + def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, + cd12, cd21, cd22): """Compute position angle error (E of N) from detector angle error. Parameters @@ -1417,24 +1419,36 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, cd12, Position angle error in degrees """ # Need to compute abs(dPA/dtheta)*theta_Err to get propogated errors - - # Get unit direction + # Unit vector in (x, y) along da 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 + # 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 - # Get derivative of tangent - d_ratio_dtheta = (v * du_dtheta - u * dv_dtheta) / v ** 2 - dPA_dtheta = (1 / (1 + ratio ** 2)) * d_ratio_dtheta + # 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) @@ -1443,7 +1457,6 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, cd12, return pa_err - class ComputePixelScale(LocalWcs): """Compute the local pixel scale from the stored CDMatrix. """ @@ -2088,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 @@ -2209,7 +2215,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 24441188d..a24703a82 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -53,8 +53,9 @@ MomentsG1Sky, MomentsG2Sky, MomentsTraceSky, CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky, SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation, - PositionAngleFromCorrelation - PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr) + PositionAngleFromCorrelation, + PositionAngleFromMoments, + ConvertDetectorAngleErrToPositionAngleErr) ROOT = os.path.abspath(os.path.dirname(__file__)) @@ -763,70 +764,6 @@ 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. @@ -924,6 +861,111 @@ 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. + """ + import pydevd_pycharm + pydevd_pycharm.settrace('localhost', port=8888, stdout_to_server=True, + stderr_to_server=True) + 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"], + ) * .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 = 1.0e-7 # 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=5e-6, + ) + def _makeWcs(self, dec=53.1595451514076, theta=0): """Create a wcs from real CFHT values, rotated by an optional theta. From a283dfc47161443f743c9d7e9d4e0af003f405b7 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Wed, 4 Feb 2026 12:13:44 -0800 Subject: [PATCH 3/6] Fix unit tests --- tests/test_functors.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/tests/test_functors.py b/tests/test_functors.py index a24703a82..1213dce02 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -867,13 +867,11 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): Requires a similar setup to testConvertDetectorAngleToPositionAngle. """ - import pydevd_pycharm - pydevd_pycharm.settrace('localhost', port=8888, stdout_to_server=True, - stderr_to_server=True) 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( @@ -902,7 +900,7 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): "slot_Centroid_y"], self.dataDict["someCentroid_x"] - self.dataDict[ "slot_Centroid_x"], - ) * .001 + ) * 0.001 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full( self.nRecords, @@ -931,7 +929,7 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): # Numerical derivative of the same PA function so a # compariosn can be made - step = 1.0e-7 # radians. Bigger? Smaller? + 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")], @@ -963,7 +961,7 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): val.to_numpy(), expected_pa_err_deg, rtol=0, - atol=5e-6, + atol=1e-6, ) def _makeWcs(self, dec=53.1595451514076, theta=0): From 2776ad2fa425bd8d4bbc1f474c2e4a4852ec4fc8 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Fri, 27 Mar 2026 11:33:15 -0700 Subject: [PATCH 4/6] Simplify position angle err --- python/lsst/pipe/tasks/functors.py | 54 +++++++----------------- tests/test_functors.py | 67 +++++++++++++++++++++++++++++- 2 files changed, 81 insertions(+), 40 deletions(-) diff --git a/python/lsst/pipe/tasks/functors.py b/python/lsst/pipe/tasks/functors.py index bc360d4a3..faee8260f 100644 --- a/python/lsst/pipe/tasks/functors.py +++ b/python/lsst/pipe/tasks/functors.py @@ -1413,49 +1413,27 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, [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 + Notes + ----- + The error is estimated by evaluating the position angle at + ``theta ± theta_err`` and taking half the absolute difference. + """ + pa_plus = self.getPositionAngleFromDetectorAngle( + theta + theta_err, cd11, cd12, cd21, cd22 + ) + pa_minus = self.getPositionAngleFromDetectorAngle( + theta - theta_err, cd11, cd12, cd21, cd22 + ) + # Wrap difference into (-180, 180] to handle the ±180 deg boundary. + delta = pa_plus - pa_minus + delta = (delta + 180) % 360 - 180 + return np.abs(delta) / 2 class ComputePixelScale(LocalWcs): """Compute the local pixel scale from the stored CDMatrix. diff --git a/tests/test_functors.py b/tests/test_functors.py index 1213dce02..71a1b7a4d 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -895,12 +895,12 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): self.dataDict["someCentroid_x"] - self.dataDict[ "slot_Centroid_x"], ) - self.dataDict["orientation_err"] = np.arctan2( + self.dataDict["orientation_err"] = np.abs(np.arctan2( self.dataDict["someCentroid_y"] - self.dataDict[ "slot_Centroid_y"], self.dataDict["someCentroid_x"] - self.dataDict[ "slot_Centroid_x"], - ) * 0.001 + )) * 0.001 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full( self.nRecords, @@ -964,6 +964,69 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): atol=1e-6, ) + def testPositionAngleErrPureRotationWcs(self): + """PA error for a pure scaling WCS is exactly rad2deg(theta_err). + + For CD = diag(s, -s) the sky direction of a pixel unit vector at + angle theta is PA(theta) = theta + pi/2 (flat-sky), so + |dPA/dtheta| = 1 everywhere and the propagated error is simply + np.rad2deg(theta_err) regardless of theta. This gives an + independent analytical expectation that does not depend on the + implementation under test. + """ + s = 5e-5 # degrees/pixel, representative of LSST + cd11, cd12, cd21, cd22 = s, 0.0, 0.0, -s + local_wcs = LocalWcs("cd11", "cd12", "cd21", "cd22") + + thetas = pd.Series([0.0, np.pi / 4, np.pi / 3, -np.pi / 6, 1.2, -0.8]) + theta_err = 0.01 # radians, small enough that linearity holds precisely + n = len(thetas) + cd11_s = pd.Series(np.full(n, cd11)) + cd12_s = pd.Series(np.full(n, cd12)) + cd21_s = pd.Series(np.full(n, cd21)) + cd22_s = pd.Series(np.full(n, cd22)) + + pa_err = local_wcs.getPositionAngleErrFromDetectorAngleErr( + thetas, theta_err, cd11_s, cd12_s, cd21_s, cd22_s + ) + + expected_pa_err_deg = np.rad2deg(theta_err) + np.testing.assert_allclose( + pa_err.to_numpy(), expected_pa_err_deg, rtol=1e-4, atol=0 + ) + + def testPositionAngleErrWrapAround(self): + """PA error is correct when PA crosses the ±180 deg boundary. + + For CD = diag(s, -s) and theta = pi/2, PA(theta) = pi deg exactly + (the wrap boundary). With a small theta_err: + pa_plus = PA(pi/2 + theta_err) wraps to -180 + rad2deg(theta_err) + pa_minus = PA(pi/2 - theta_err) = +180 - rad2deg(theta_err) + + Naive subtraction gives ~-360 deg → apparent error ~180 deg (wrong). + The wrap-corrected implementation should return rad2deg(theta_err). + """ + s = 5e-5 # degrees/pixel + cd11, cd12, cd21, cd22 = s, 0.0, 0.0, -s + local_wcs = LocalWcs("cd11", "cd12", "cd21", "cd22") + + theta_center = np.pi / 2 # PA(theta_center) = 180 deg → wrap boundary + theta_err = 0.01 # radians + + cd_s = lambda v: pd.Series([v]) # noqa: E731 + pa_err = local_wcs.getPositionAngleErrFromDetectorAngleErr( + pd.Series([theta_center]), + theta_err, + cd_s(cd11), cd_s(cd12), cd_s(cd21), cd_s(cd22), + ) + + expected_pa_err_deg = np.rad2deg(theta_err) + # Correct result is ~0.57 deg; naive (unwrapped) result would be ~180 deg. + np.testing.assert_allclose( + pa_err.to_numpy(), expected_pa_err_deg, rtol=1e-4, atol=0 + ) + self.assertLess(float(pa_err.iloc[0]), 1.0) + def _makeWcs(self, dec=53.1595451514076, theta=0): """Create a wcs from real CFHT values, rotated by an optional theta. From 2a04c3b661bd4d0a71bb41fdcea9cffd5b9a7fd9 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Fri, 27 Mar 2026 12:21:04 -0700 Subject: [PATCH 5/6] Fixes --- python/lsst/pipe/tasks/functors.py | 4 +++- tests/test_functors.py | 4 +--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/lsst/pipe/tasks/functors.py b/python/lsst/pipe/tasks/functors.py index faee8260f..179e0f4b1 100644 --- a/python/lsst/pipe/tasks/functors.py +++ b/python/lsst/pipe/tasks/functors.py @@ -28,7 +28,7 @@ "ComputePixelScale", "ConvertPixelToArcseconds", "ConvertPixelSqToArcsecondsSq", "ConvertDetectorAngleToPositionAngle", - "ConvertDetectorAngleErrToPositionAngleErr," + "ConvertDetectorAngleErrToPositionAngleErr", "ReferenceBand", "Photometry", "NanoJansky", "NanoJanskyErr", "LocalPhotometry", "LocalNanojansky", "LocalNanojanskyErr", "LocalDipoleMeanFlux", @@ -1435,6 +1435,7 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, delta = (delta + 180) % 360 - 180 return np.abs(delta) / 2 + class ComputePixelScale(LocalWcs): """Compute the local pixel scale from the stored CDMatrix. """ @@ -2193,6 +2194,7 @@ 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 71a1b7a4d..3d0a9a2df 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -54,7 +54,6 @@ CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky, SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation, PositionAngleFromCorrelation, - PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr) ROOT = os.path.abspath(os.path.dirname(__file__)) @@ -871,7 +870,6 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): 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( @@ -929,7 +927,7 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): # Numerical derivative of the same PA function so a # compariosn can be made - step = self.dataDict["orientation_err"] # radians. Bigger? Smaller? + 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")], From a3bc662cc5d734bffa7e73d7fb4c9553217872c1 Mon Sep 17 00:00:00 2001 From: Brianna Smart Date: Wed, 1 Apr 2026 14:30:49 -0700 Subject: [PATCH 6/6] Update tests --- tests/test_functors.py | 107 ++++++++++++++++++++++++++--------------- 1 file changed, 69 insertions(+), 38 deletions(-) diff --git a/tests/test_functors.py b/tests/test_functors.py index 3d0a9a2df..b33883cab 100644 --- a/tests/test_functors.py +++ b/tests/test_functors.py @@ -925,42 +925,12 @@ def testConvertDetectorAngleErrToPositionAngleErr(self): ) 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, - ) + # Errors must be non-negative and finite across all + # WCS configurations, declinations, and pixel positions. + self.assertTrue(np.all(val.to_numpy() >= 0), + "PA errors must be non-negative") + self.assertTrue(np.all(np.isfinite(val.to_numpy())), + "PA errors must be finite") def testPositionAngleErrPureRotationWcs(self): """PA error for a pure scaling WCS is exactly rad2deg(theta_err). @@ -972,12 +942,12 @@ def testPositionAngleErrPureRotationWcs(self): independent analytical expectation that does not depend on the implementation under test. """ - s = 5e-5 # degrees/pixel, representative of LSST + s = 5e-5 # degrees/pixel cd11, cd12, cd21, cd22 = s, 0.0, 0.0, -s local_wcs = LocalWcs("cd11", "cd12", "cd21", "cd22") thetas = pd.Series([0.0, np.pi / 4, np.pi / 3, -np.pi / 6, 1.2, -0.8]) - theta_err = 0.01 # radians, small enough that linearity holds precisely + theta_err = 0.01 # radians n = len(thetas) cd11_s = pd.Series(np.full(n, cd11)) cd12_s = pd.Series(np.full(n, cd12)) @@ -1025,6 +995,67 @@ def testPositionAngleErrWrapAround(self): ) self.assertLess(float(pa_err.iloc[0]), 1.0) + def testPositionAngleErrAccuracy(self): + """Validate PA error accuracy by comparing to a Monte Carlo distribution. + + Draw n_samples values of theta from N(theta_center, theta_err), compute + the position angle for each sample using getPositionAngleFromDetectorAngle, + and compare the circular standard deviation of the resulting PA + distribution to the functor output. + + Two CD matrices are used to test normal scaling vs shearing: + - Pure scaling (CD = diag(s, -s)): PA(theta) = theta + pi/2, Jacobian = 1. + - Sheared WCS: Jacobian varies with theta, exercising the non-linear path. + """ + rng = np.random.default_rng(42) + n_samples = 2000 + theta_err = 0.01 + tolerance = 0.05 # allow 5 % relative error + + s = 5.5e-5 # degrees/pixel. Do I need to be super accurate? + # Need to test against both?? Probably??? + cd_matrices = { + "pure_scaling": np.array([[s, 0.0], [0.0, -s]]), + "sheared": np.array([[s, 0.4 * s], [-0.2 * s, -s]]), + } + + for label, cd in cd_matrices.items(): + cd11, cd12, cd21, cd22 = cd[0, 0], cd[0, 1], cd[1, 0], cd[1, 1] + local_wcs = LocalWcs("cd11", "cd12", "cd21", "cd22") + + for theta_center in [0.0, 0.5, -1.2]: + theta_samples = rng.normal(theta_center, theta_err, size=n_samples) + n = len(theta_samples) + + pa_samples_deg = local_wcs.getPositionAngleFromDetectorAngle( + pd.Series(theta_samples), + pd.Series(np.full(n, cd11)), + pd.Series(np.full(n, cd12)), + pd.Series(np.full(n, cd21)), + pd.Series(np.full(n, cd22)), + ) + pa_rad = np.deg2rad(pa_samples_deg.to_numpy()) + # Make sure circular standard deviation handles the +/-180 deg + # wrap correctly. + R = np.abs(np.mean(np.exp(1j * pa_rad))) + mc_std_deg = np.rad2deg(np.sqrt(-2.0 * np.log(R))) + + pa_err = local_wcs.getPositionAngleErrFromDetectorAngleErr( + pd.Series([theta_center]), + theta_err, + pd.Series([cd11]), + pd.Series([cd12]), + pd.Series([cd21]), + pd.Series([cd22]), + ) + + self.assertAlmostEqual( + float(pa_err.iloc[0]), + mc_std_deg, + delta=mc_std_deg * tolerance, + msg=f"MC vs functor mismatch for {label}, theta_center={theta_center}", + ) + def _makeWcs(self, dec=53.1595451514076, theta=0): """Create a wcs from real CFHT values, rotated by an optional theta.