Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 96 additions & 0 deletions python/lsst/pipe/tasks/functors.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.

Expand Down
67 changes: 65 additions & 2 deletions tests/test_functors.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,10 @@
ConvertDetectorAngleToPositionAngle,
HtmIndex20, Ebv, MomentsIuuSky, MomentsIvvSky, MomentsIuvSky,
SemimajorAxisFromMoments, SemiminorAxisFromMoments,
PositionAngleFromMoments,
CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
PositionAngleFromCorrelation
)
PositionAngleFromMoments, ConvertDetectorAngleErrToPositionAngleErr)

Check failure on line 55 in tests/test_functors.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

E999

SyntaxError: invalid syntax

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

Expand Down Expand Up @@ -762,6 +761,70 @@

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.
Expand Down
Loading