Skip to content

Commit dfbd201

Browse files
bsmartradiobsmartradio
committed
Test
Fix Test Fixes
1 parent 687798f commit dfbd201

2 files changed

Lines changed: 134 additions & 87 deletions

File tree

‎python/lsst/pipe/tasks/functors.py‎

Lines changed: 26 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
"ComputePixelScale", "ConvertPixelToArcseconds",
2929
"ConvertPixelSqToArcsecondsSq",
3030
"ConvertDetectorAngleToPositionAngle",
31+
"ConvertDetectorAngleErrToPositionAngleErr,"
3132
"ReferenceBand", "Photometry",
3233
"NanoJansky", "NanoJanskyErr", "LocalPhotometry", "LocalNanojansky",
3334
"LocalNanojanskyErr", "LocalDipoleMeanFlux",
@@ -1394,7 +1395,8 @@ def getPositionAngleFromDetectorAngle(self, theta, cd11, cd12, cd21, cd22):
13941395
# Position angle of vector from (RA1, Dec1) to (RA2, Dec2)
13951396
return np.rad2deg(self.computePositionAngle(ra1, dec1, ra2, dec2))
13961397

1397-
def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, cd12, cd21, cd22):
1398+
def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11,
1399+
cd12, cd21, cd22):
13981400
"""Compute position angle error (E of N) from detector angle error.
13991401
14001402
Parameters
@@ -1417,24 +1419,36 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, cd12,
14171419
Position angle error in degrees
14181420
"""
14191421
# Need to compute abs(dPA/dtheta)*theta_Err to get propogated errors
1420-
1421-
# Get unit direction
1422+
# Unit vector in (x, y) along da
14221423
dx = np.cos(theta)
14231424
dy = np.sin(theta)
1424-
1425-
# Transform it using WCS?
14261425
u = dx * cd11 + dy * cd12
14271426
v = dx * cd21 + dy * cd22
1428-
# Now we are computing the tangent
1429-
ratio = u / v
14301427

1431-
# Get derivative of theta
1428+
# Derivative of position angle wrt detector angle
14321429
du_dtheta = -np.sin(theta) * cd11 + np.cos(theta) * cd12
14331430
dv_dtheta = -np.sin(theta) * cd21 + np.cos(theta) * cd22
14341431

1435-
# Get derivative of tangent
1436-
d_ratio_dtheta = (v * du_dtheta - u * dv_dtheta) / v ** 2
1437-
dPA_dtheta = (1 / (1 + ratio ** 2)) * d_ratio_dtheta
1432+
# Precomputing sin/cos values
1433+
sin_u = np.sin(u)
1434+
cos_u = np.cos(u)
1435+
sin_v = np.sin(v)
1436+
cos_v = np.cos(v)
1437+
1438+
# PA is atan2(Y, X), for great circle
1439+
pa_y = sin_u * cos_v
1440+
pa_x = sin_v
1441+
1442+
# We need dPA/dtheta for error propogation.
1443+
# X' = cos(v) * v'
1444+
dX_dtheta = cos_v * dv_dtheta
1445+
# Y' = (cos(u)*cos(v))*u' + (sin(u)*(-sin(v)))*v'
1446+
dY_dtheta = (cos_u * cos_v) * du_dtheta - (sin_u * sin_v) * dv_dtheta
1447+
1448+
denom = pa_x * pa_x + pa_y * pa_y
1449+
1450+
dPA_dtheta = (pa_x * dY_dtheta - pa_y * dX_dtheta) / denom
1451+
dPA_dtheta = np.where(denom == 0, np.nan, dPA_dtheta)
14381452

14391453
pa_err = np.rad2deg(np.abs(dPA_dtheta) * theta_err)
14401454

@@ -1443,7 +1457,6 @@ def getPositionAngleErrFromDetectorAngleErr(self, theta, theta_err, cd11, cd12,
14431457

14441458
return pa_err
14451459

1446-
14471460
class ComputePixelScale(LocalWcs):
14481461
"""Compute the local pixel scale from the stored CDMatrix.
14491462
"""
@@ -2088,14 +2101,7 @@ def _func(self, df):
20882101

20892102

20902103
class MomentsBase(Functor):
2091-
"""Base class for functors that use shape moments and localWCS
2092-
2093-
Attributes
2094-
----------
2095-
is_covariance : bool
2096-
Whether the shape columns are terms of a covariance matrix. If False,
2097-
they will be assumed to be terms of a correlation matrix instead.
2098-
"""
2104+
"""Base class for functors that use shape moments and localWCS"""
20992105

21002106
is_covariance: bool = True
21012107

@@ -2211,7 +2217,6 @@ def sky_uv(self, df):
22112217
CD_2_2 = df[self.colCD_2_2]
22122218
return ((CD_1_1 * i_xx + CD_1_2 * i_xy) * CD_2_1
22132219
+ (CD_1_1 * i_xy + CD_1_2 * i_yy) * CD_2_2)
2214-
22152220
def get_g1(self, df):
22162221
"""
22172222
Calculate shear-type ellipticity parameter G1.

‎tests/test_functors.py‎

Lines changed: 108 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -53,8 +53,9 @@
5353
MomentsG1Sky, MomentsG2Sky, MomentsTraceSky,
5454
CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
5555
SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
56-
PositionAngleFromCorrelation
57-
PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr)
56+
PositionAngleFromCorrelation,
57+
PositionAngleFromMoments,
58+
ConvertDetectorAngleErrToPositionAngleErr)
5859

5960
ROOT = os.path.abspath(os.path.dirname(__file__))
6061

@@ -763,70 +764,6 @@ def testConvertDetectorAngleToPositionAngle(self):
763764

764765
np.testing.assert_allclose(coord_diff, 0, rtol=0, atol=5e-6)
765766

766-
# Test position angle error propogation
767-
def testConvertDetectorAngleErrToPositionAngleErr(self):
768-
"""Test conversion of position angle err in detector degrees to
769-
position angle erron sky
770-
"""
771-
dipoleSep = 10
772-
ixx = 10
773-
testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
774-
# testConvertPixelToArcSecond uses a meas_base LocalWcsPlugin
775-
# but we're using a simple WCS as our example, so this doesn't really matter
776-
# and we'll just use the WCS directly
777-
for dec in np.linspace(-90, 90, 10):
778-
for theta in (-45, 0, 90):
779-
for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
780-
np.random.uniform(2 * 560.018167811613, size=10)):
781-
wcs = self._makeWcs(dec=dec, theta=theta)
782-
cd_matrix = wcs.getCdMatrix()
783-
784-
self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
785-
self.dataDict["ixx"] = np.full(self.nRecords, ixx)
786-
self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
787-
self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
788-
self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
789-
self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
790-
self.dataDict["orientation"] = np.arctan2(
791-
self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"],
792-
self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
793-
)
794-
self.dataDict["orientation_err"] = np.arctan2(
795-
self.dataDict["someCentroid_y"] - self.dataDict[ "slot_Centroid_y"],
796-
self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
797-
)*.001
798-
799-
self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
800-
cd_matrix[0, 0])
801-
self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
802-
cd_matrix[0, 1])
803-
self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
804-
cd_matrix[1, 0])
805-
self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
806-
cd_matrix[1, 1])
807-
df = self.getMultiIndexDataFrame(self.dataDict)
808-
809-
# Test detector angle to position angle conversion
810-
func = ConvertDetectorAngleErrToPositionAngleErr(
811-
"orientation",
812-
"orientation_err",
813-
"base_LocalWcs_CDMatrix_1_1",
814-
"base_LocalWcs_CDMatrix_1_2",
815-
"base_LocalWcs_CDMatrix_2_1",
816-
"base_LocalWcs_CDMatrix_2_2"
817-
)
818-
819-
func_pa = ConvertDetectorAngleToPositionAngle(
820-
"orientation",
821-
"base_LocalWcs_CDMatrix_1_1",
822-
"base_LocalWcs_CDMatrix_1_2",
823-
"base_LocalWcs_CDMatrix_2_1",
824-
"base_LocalWcs_CDMatrix_2_2"
825-
)
826-
val = self._funcVal(func, df)
827-
val_pa = self._funcVal(func_pa, df)
828-
829-
830767
def testConvertPixelToArcseconds(self):
831768
"""Test calculations of the pixel scale, conversions of pixel to
832769
arcseconds.
@@ -924,6 +861,111 @@ def testConvertPixelToArcseconds(self):
924861
atol=1e-16,
925862
rtol=1e-16))
926863

864+
def testConvertDetectorAngleErrToPositionAngleErr(self):
865+
"""Test conversion of position angle err in detector degrees to
866+
position angle err on sky.
867+
868+
Requires a similar setup to testConvertDetectorAngleToPositionAngle.
869+
"""
870+
import pydevd_pycharm
871+
pydevd_pycharm.settrace('localhost', port=8888, stdout_to_server=True,
872+
stderr_to_server=True)
873+
dipoleSep = 10
874+
ixx = 10
875+
testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
876+
877+
for dec in np.linspace(-80, 80, 10):
878+
for theta in (-35, 0, 90):
879+
for x, y in zip(
880+
np.random.uniform(2 * 1109.99981456774, size=10),
881+
np.random.uniform(2 * 560.018167811613, size=10)):
882+
wcs = self._makeWcs(dec=dec, theta=theta)
883+
cd_matrix = wcs.getCdMatrix()
884+
885+
self.dataDict["dipoleSep"] = np.full(self.nRecords,
886+
dipoleSep)
887+
self.dataDict["ixx"] = np.full(self.nRecords, ixx)
888+
self.dataDict["slot_Centroid_x"] = np.full(self.nRecords,
889+
x)
890+
self.dataDict["slot_Centroid_y"] = np.full(self.nRecords,
891+
y)
892+
self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
893+
self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
894+
self.dataDict["orientation"] = np.arctan2(
895+
self.dataDict["someCentroid_y"] - self.dataDict[
896+
"slot_Centroid_y"],
897+
self.dataDict["someCentroid_x"] - self.dataDict[
898+
"slot_Centroid_x"],
899+
)
900+
self.dataDict["orientation_err"] = np.arctan2(
901+
self.dataDict["someCentroid_y"] - self.dataDict[
902+
"slot_Centroid_y"],
903+
self.dataDict["someCentroid_x"] - self.dataDict[
904+
"slot_Centroid_x"],
905+
) * .001
906+
907+
self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(
908+
self.nRecords,
909+
cd_matrix[0, 0])
910+
self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(
911+
self.nRecords,
912+
cd_matrix[0, 1])
913+
self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(
914+
self.nRecords,
915+
cd_matrix[1, 0])
916+
self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(
917+
self.nRecords,
918+
cd_matrix[1, 1])
919+
df = self.getMultiIndexDataFrame(self.dataDict)
920+
921+
# Test detector angle err to position angle err conversion
922+
func = ConvertDetectorAngleErrToPositionAngleErr(
923+
"orientation",
924+
"orientation_err",
925+
"base_LocalWcs_CDMatrix_1_1",
926+
"base_LocalWcs_CDMatrix_1_2",
927+
"base_LocalWcs_CDMatrix_2_1",
928+
"base_LocalWcs_CDMatrix_2_2"
929+
)
930+
val = self._funcVal(func, df)
931+
932+
# Numerical derivative of the same PA function so a
933+
# compariosn can be made
934+
step = 1.0e-7 # radians. Bigger? Smaller?
935+
pa_plus_deg = func.getPositionAngleFromDetectorAngle(
936+
df[("meas", "g", "orientation")] + step,
937+
df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
938+
df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
939+
df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
940+
df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")],
941+
)
942+
pa_minus_deg = func.getPositionAngleFromDetectorAngle(
943+
df[("meas", "g", "orientation")] - step,
944+
df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
945+
df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
946+
df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
947+
df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")],
948+
)
949+
950+
pa_plus = np.deg2rad(pa_plus_deg.to_numpy())
951+
pa_minus = np.deg2rad(pa_minus_deg.to_numpy())
952+
953+
# From example: smallest signed angular difference in
954+
# (-pi, +pi]
955+
dpa = np.angle(np.exp(1j * (pa_plus - pa_minus)))
956+
dpa_dtheta = dpa / (2.0 * step)
957+
958+
theta_err = df[("meas", "g", "orientation_err")].to_numpy()
959+
expected_pa_err_deg = np.rad2deg(
960+
np.abs(dpa_dtheta) * theta_err)
961+
962+
np.testing.assert_allclose(
963+
val.to_numpy(),
964+
expected_pa_err_deg,
965+
rtol=0,
966+
atol=5e-6,
967+
)
968+
927969
def _makeWcs(self, dec=53.1595451514076, theta=0):
928970
"""Create a wcs from real CFHT values, rotated by an optional theta.
929971

0 commit comments

Comments
 (0)