diff --git a/src/ipyaladin/region_converter.py b/src/ipyaladin/region_converter.py index ee344553..1dd2b03a 100644 --- a/src/ipyaladin/region_converter.py +++ b/src/ipyaladin/region_converter.py @@ -1,8 +1,7 @@ -import math - import numpy as np -from astropy.coordinates import SkyCoord +from astropy.coordinates import SkyCoord, CartesianRepresentation, Angle +from astropy.coordinates.matrix_utilities import rotation_matrix try: from regions import ( @@ -22,108 +21,10 @@ LineSkyRegion = None from typing import Union -TWICE_PI = 2 * math.pi - - -class RefToLocalRotMatrix: - """A rotation matrix from the reference frame to the local frame. - - The reference frame is the ICRS frame. The local frame is - defined by the center of the local reference frame. - - Attributes - ---------- - rotation_matrix : np.ndarray[float] - The rotation matrix. - - """ - - def __init__( - self, - rotation_matrix: np.ndarray[float], - ) -> None: - self.rotation_matrix = rotation_matrix - - @classmethod - def from_center(cls: any, lon: float, lat: float) -> "RefToLocalRotMatrix": - """Create a rotation matrix from the center of the local reference frame. - - Parameters - ---------- - lon : float - The longitude of the center. - lat : float - The latitude of the center. - - Returns - ------- - RefToLocalRotMatrix - The rotation matrix. - - """ - ca, sa = math.cos(lon), math.sin(lon) - cd, sd = math.cos(lat), math.sin(lat) - return cls( - np.array([[ca * cd, sa * cd, sd], [-sa, ca, 0.0], [-ca * sd, -sa * sd, cd]]) - ) - - def to_global_xyz(self, x: float, y: float, z: float) -> tuple[float, float, float]: - """Convert local cartesian coordinates to global cartesian coordinates. - - Parameters - ---------- - x : float - The x coordinate. - y : float - The y coordinate. - z : float - The z coordinate. - - Returns - ------- - tuple[float, float, float] - The x, y, and z coordinates. - - """ - rotated_matrix = np.sum( - self.rotation_matrix * np.array([[x], [y], [z]]), axis=0 - ) - return rotated_matrix[0], rotated_matrix[1], rotated_matrix[2] - - def to_global_coo(self, x: float, y: float, z: float) -> tuple[float, float]: - """Convert local cartesian coordinates to global spherical coordinates. - - Parameters - ---------- - x : float - The x coordinate. - y : float - The y coordinate. - z : float - The z coordinate. - - Returns - ------- - tuple[float, float] - The longitude and latitude. - - """ - x, y, z = self.to_global_xyz(x, y, z) - # Convert cartesian coordinates to spherical coordinates. - r2 = x * x + y * y - lat = math.atan2(z, math.sqrt(r2)) - lon = math.atan2(y, x) - if lon < 0: - lon += TWICE_PI - return lon, lat - def rectangle_to_polygon_region(region: RectangleSkyRegion) -> PolygonSkyRegion: """Convert a RectangleSkyRegion to a PolygonSkyRegion. - Converted code from cdshealpix: - https://github.com/cds-astro/cds-healpix-rust/blob/master/src/nested/mod.rs#L3687-3742 - Parameters ---------- region : RectangleSkyRegion @@ -135,47 +36,50 @@ def rectangle_to_polygon_region(region: RectangleSkyRegion) -> PolygonSkyRegion: The converted region. """ - longitude, latitude = region.center.icrs.ra.rad, region.center.icrs.dec.rad - width = region.width.rad - height = region.height.rad - position_angle = region.angle.rad - - if not 0.0 <= longitude < TWICE_PI: - raise ValueError(f"Expected: lon in [0, 2pi[. Actual: {longitude}") - if not -math.pi / 2 <= latitude <= math.pi / 2: - raise ValueError(f"Expected: lat in [-pi/2, pi/2]. Actual: {latitude}") - if not 0.0 < width <= math.pi / 2: - raise ValueError(f"Expected: a in ]0, pi/2]. Actual: {width}") - if not 0.0 < height <= width: - raise ValueError(f"Expected: b in ]0, a]. Actual: {height}") - if not 0.0 <= position_angle < math.pi: - raise ValueError(f"Expected: pa in [0, pi[. Actual: {position_angle}") - - frame_rotation = RefToLocalRotMatrix.from_center(longitude, latitude) - - sin_lon, cos_lon = math.sin(width), math.cos(width) - latitude = math.atan(cos_lon * math.tan(height)) - sin_lat, cos_lat = math.sin(latitude), math.cos(latitude) - sin_pa, cos_pa = math.sin(position_angle), math.cos(position_angle) - - x1, y1, z1 = cos_lon * cos_lat, sin_lon * cos_lat, sin_lat + bl_x, bl_y, bl_z = ( + region.center.spherical_offsets_by(-region.width / 2, -region.height / 2) + .represent_as("cartesian") + .xyz.value + ) + br_x, br_y, br_z = ( + region.center.spherical_offsets_by(-region.width / 2, region.height / 2) + .represent_as("cartesian") + .xyz.value + ) + tr_x, tr_y, tr_z = ( + region.center.spherical_offsets_by(region.width / 2, region.height / 2) + .represent_as("cartesian") + .xyz.value + ) + rl_x, tly, tl_z = ( + region.center.spherical_offsets_by(region.width / 2, -region.height / 2) + .represent_as("cartesian") + .xyz.value + ) + c_x, c_y, c_z = region.center.represent_as("cartesian").xyz.value - pa_matrix = np.array([[y1 * sin_pa, z1 * sin_pa], [y1 * cos_pa, z1 * sin_pa]]) + rot_mat = rotation_matrix(Angle(-region.angle.deg, unit="deg"), (c_x, c_y, c_z)) - multipliers = [ - [[1, -1], [1, 1]], - [[1, 1], [1, -1]], - [[-1, 1], [-1, -1]], - [[-1, -1], [-1, 1]], + corners = np.array( + [ + [bl_x, bl_y, bl_z], + [br_x, br_y, br_z], + [tr_x, tr_y, tr_z], + [rl_x, tly, tl_z], + ] + ) + rotated_corners = np.dot(corners, rot_mat) + vertices = [ + SkyCoord(CartesianRepresentation(*corner), frame="icrs") + for corner in rotated_corners ] - vertices = [] - - for multiplier in multipliers: - y2, z2 = np.sum(np.array(multiplier) * pa_matrix, axis=1) - vertices.append(frame_rotation.to_global_coo(x1, y2, z2)) return PolygonSkyRegion( - vertices=SkyCoord(vertices, unit="rad", frame="icrs"), + vertices=SkyCoord( + [coord.icrs.ra.deg for coord in vertices], + [coord.icrs.dec.deg for coord in vertices], + unit="deg", + ), visual=region.visual, meta=region.meta, )