diff --git a/python/lsst/summit/utils/coordinates.py b/python/lsst/summit/utils/coordinates.py new file mode 100644 index 00000000..4f1ab6a3 --- /dev/null +++ b/python/lsst/summit/utils/coordinates.py @@ -0,0 +1,243 @@ +# This file is part of summit_utils. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import logging + +import numpy as np + +from lsst.afw.geom.ellipses import Quadrupole +from lsst.geom import LinearTransform + +__all__ = ["convertE1E2", "convertIxxIyyIxy"] + +logger = logging.getLogger(__name__) + + +def convertE1E2( + e1: list, + e2: list, + physical_rotator_angle: float, + input_coordinate_system: str, + output_coordinate_system: str, + sky_angle: float | None = None, +) -> list: + """Transform ellipticities between coordinate systems for LSSTCam images. + + Given an input e1, e2 with a specified coordinate system and rotator + angle (and sky angle if needed), convert the ellipticity between CCS, + DVCS, OCS, Alt/Az, equatorial, and North-West (NW) coordinate systems. + + See https://smtn-019.lsst.io/#on-sky-and-hardware-rotation-angles + for more details. + + Parameters + ---------- + e1 : `list` + A list of input e1 ellipticities. + e2 : `list` + A list of input e2 ellipticities. + physical_rotator_angle : `float` + Camera rotator angle in degrees. Equivalent to the rotator telescope + position (RotTelPos or rtp). + input_coordinate_system : `string` + The input coordinate system, which is the system that e1 and e2 are in. + Options are: 'altaz', 'ccs', 'dvcs', 'equatorial', 'nw' 'ocs'. + output_coordinate_system : `string` + The input coordinate system, which is the system that e1 and e2 are in. + Options are: 'altaz', 'ccs', 'dvcs', 'equatorial', 'nw', 'ocs'. + sky_angle : `float` | None + Optional: Sky angle in degrees. This is the orientation of the +đť‘ŚDVCS + axis (projected on the sky) measured east of north in the International + Celestial Reference Frame (ICRF). Also called RotSkyPos or given by the + boresightRotAngle if the rotType=Sky (which is true for all FBS + observations and many engineering observations). Only needed if the + coordinate system includes equatorial or nw. + + Returns + ------- + e1_out : `list` + A list of output e1 ellipticities. + e2_out : `list` + A list of output e2 ellipticities. + """ + + allowed_coordinates = ["altaz", "ccs", "dvcs", "ocs", "equatorial", "nw"] + + if input_coordinate_system not in allowed_coordinates: + raise ValueError(f"Input coordinate system must be in {allowed_coordinates}") + if output_coordinate_system not in allowed_coordinates: + raise ValueError(f"Output coordinate system must be in {allowed_coordinates}") + if (("nw" or "equatorial") in [input_coordinate_system, output_coordinate_system]) and ( + sky_angle is None + ): + raise ValueError("Sky angle must be provided when using nw or equatorial coordinate systems.") + + if input_coordinate_system == output_coordinate_system: + return e1, e2 + + # convert e1 and e2 into Ixx, Iyy, Ixy components + # letting the trace T (Ixx + Iyy) be 1. + Ixx = [(1 + e1_input) / 2 for e1_input in e1] + Iyy = [(1 - e1_input) / 2 for e1_input in e1] + Ixy = [(e2_input) / 2 for e2_input in e2] + + if input_coordinate_system == "nw": + # if input coordinate system is "nw", immediately convert to dvcs + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx, Iyy, Ixy, -1 * sky_angle) + if output_coordinate_system == "dvcs": + return e1_out, e2_out + else: + input_coordinate_system = "dvcs" + + if input_coordinate_system == "dvcs": + if output_coordinate_system == "ccs": + # for ccs you swap x and y relative to dvcs + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Iyy, Ixx, Ixy, 0) + return e1_out, e2_out + elif output_coordinate_system in ["ocs", "altaz"]: + # for ocs and altaz you swap x and y relative to dvcs + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Iyy, Ixx, Ixy, -1 * physical_rotator_angle) + if output_coordinate_system == "ocs": + return e1_out, e2_out + else: + return e1_out, [-1 * ellipticity for ellipticity in e2_out] + elif output_coordinate_system == "nw": + # for NW coordinate system, you don't need to transpose x and y + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx, Iyy, Ixy, sky_angle) + return e1_out, e2_out + + if input_coordinate_system == "altaz": + if output_coordinate_system == "ocs": + return e1, [-1 * ellipticity for ellipticity in e2] + elif output_coordinate_system in ["ccs", "dvcs", "nw"]: + _, _, _, e1_out, e2_out = convertIxxIyyIxy( + Ixx, Iyy, [-1 * xy for xy in Ixy], physical_rotator_angle + ) + if output_coordinate_system == "ccs": + return e1_out, e2_out + else: + e1_dvcs = [-1 * ellipticity for ellipticity in e1_out] + if output_coordinate_system == "dvcs": + return e1_dvcs, e2_out + else: + # rotate into NW using sky angle + Ixx_dvcs = [(1 + dvcs) / 2 for dvcs in e1_dvcs] + Iyy_dvcs = [(1 - dvcs) / 2 for dvcs in e1_dvcs] + Ixy_dvcs = [(e2_input) / 2 for e2_input in e2_out] + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx_dvcs, Iyy_dvcs, Ixy_dvcs, sky_angle) + return e1_out, e2_out + + elif input_coordinate_system == "ccs": + if output_coordinate_system in ["altaz", "ocs"]: + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx, Iyy, Ixy, -1 * physical_rotator_angle) + if output_coordinate_system == "ocs": + return e1_out, e2_out + else: + return e1_out, [-1 * ellipticity for ellipticity in e2_out] + elif output_coordinate_system in ["dvcs", "nw"]: + e1_dvcs = [-1 * ellipticity for ellipticity in e1_out] + if output_coordinate_system == "dvcs": + return e1_dvcs, e2_out + else: + # rotate into NW using sky angle + Ixx_dvcs = [(1 + dvcs) / 2 for dvcs in e1_dvcs] + Iyy_dvcs = [(1 - dvcs) / 2 for dvcs in e1_dvcs] + Ixy_dvcs = [(e2_input) / 2 for e2_input in e2_out] + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx_dvcs, Iyy_dvcs, Ixy_dvcs, sky_angle) + return e1_out, e2_out + + elif input_coordinate_system == "ocs": + if output_coordinate_system == "altaz": + return e1, [-1 * ellipticity for ellipticity in e2] + elif output_coordinate_system in ["ccs", "dvcs", "nw"]: + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx, Iyy, Ixy, physical_rotator_angle) + if output_coordinate_system == "ccs": + return e1_out, e2_out + else: + e1_dvcs = [-1 * ellipticity for ellipticity in e1_out] + if output_coordinate_system == "dvcs": + return e1_dvcs, e2_out + else: + # rotate into NW using sky angle + Ixx_dvcs = [(1 + dvcs) / 2 for dvcs in e1_dvcs] + Iyy_dvcs = [(1 - dvcs) / 2 for dvcs in e1_dvcs] + Ixy_dvcs = [(e2_input) / 2 for e2_input in e2_out] + _, _, _, e1_out, e2_out = convertIxxIyyIxy(Ixx_dvcs, Iyy_dvcs, Ixy_dvcs, sky_angle) + return e1_out, e2_out + + return [], [] + + +def convertIxxIyyIxy( + Ixx: list, + Iyy: list, + Ixy: list, + rotation_angle: float, +) -> list: + """Rotate Ixx, Iyy, Ixy by a given rotation angle. + + Parameters + ---------- + Ixx : `list` + A list of input Ixx quantities. + Iyy : `list` + A list of input Iyy quantities. + Ixy : `list` + A list of input Ixy quantities. + rotation_angle : `float` + Angle of rotation in degrees. + + Returns + ------- + Ixx_out : `list` + A list of output Ixx quantities. + Iyy_out : `list` + A list of output Iyy quantities. + Ixy_out : `list` + A list of output Ixy quantities. + e1_out : `list` + A list of output e1 quantities + e2_out : `list + A list of output e2 quantitites + """ + + # Build rotation matrix + rotation_angle_rad = rotation_angle * np.pi / 180 + srot, crot = np.sin(rotation_angle_rad), np.cos(rotation_angle_rad) + rot = np.array([[crot, -srot], [srot, crot]]) @ np.array([[0, 1], [1, 0]]) @ np.array([[-1, 0], [0, 1]]) + + # perform the rotations + transform = LinearTransform(rot) + rotShapes = [] + for xx, yy, xy in [*zip(Ixx, Iyy, Ixy)]: + shape = Quadrupole(xx, yy, xy) + rotShape = shape.transform(transform) + rotShapes.append(rotShape) + Ixx_out = [sh.getIxx() for sh in rotShapes] + Iyy_out = [sh.getIyy() for sh in rotShapes] + Ixy_out = [sh.getIxy() for sh in rotShapes] + + T_out = [xx + yy for xx, yy in [*zip(Ixx_out, Iyy_out)]] + + e1_out = [(xx + yy) / T for xx, yy, T in [*zip(Ixx_out, Iyy_out, T_out)]] + e2_out = [2 * xy / T for xy, T in [*zip(Ixy_out, T_out)]] + + return Ixx_out, Iyy_out, Ixy_out, e1_out, e2_out diff --git a/tests/test_coordinates.py b/tests/test_coordinates.py new file mode 100644 index 00000000..85a10dda --- /dev/null +++ b/tests/test_coordinates.py @@ -0,0 +1,26 @@ +import unittest + +import lsst.utils.tests +from lsst.summit.utils.coordinates import convertE1E2 + + +class CoordinatesTestCase(lsst.utils.tests.TestCase): + def test_convertE1E2(self) -> None: + e1 = [0, 0, 1, 0.5] + e2 = [0, 0, 0, 0.5] + e1_out, e2_out = convertE1E2(e1, e2, 0, "ccs", "ocs") + self.assertEqual(e1, e1_out) + self.assertEqual(e2, e2_out) + + +class TestMemory(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module: object) -> None: + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main()