Merge pull request #91 from astrofrog/circle-to-pix-full

Initial implementation of CircleSkyRegion.to_pixel in full mode

regions/_geometry/__init__.py

@@ -3,6 +3,7 @@
 Geometry subpackage for low-level geometry functions.
 """
 
+from .rotate_polygon import *
 from .circular_overlap import *
 from .elliptical_overlap import *
 from .rectangular_overlap import *

regions/_geometry/rotate_polygon.py

@@ -0,0 +1,40 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+import numpy as np
+
+from astropy import units as u
+from astropy.coordinates.representation import UnitSphericalRepresentation
+
+try:
+    from astropy.coordinates.matrix_utilities import rotation_matrix as rotation_array
+    def rotation_matrix(*args, **kwargs):
+        return np.matrix(rotation_array(*args, **kwargs))
+except ImportError:  # Astropy < 1.3
+    from astropy.coordinates.angles import rotation_matrix
+
+__all__ = ['rotate_polygon']
+
+
+def rotate_polygon(lon, lat, lon0, lat0):
+    """
+    Given a polygon with vertices defined by (lon, lat), rotate the polygon
+    such that the North pole of the spherical coordinates is now at (lon0,
+    lat0). Therefore, to end up with a polygon centered on (lon0, lat0), the
+    polygon should initially be drawn around the North pole.
+    """
+
+    # Create a representation object
+    polygon = UnitSphericalRepresentation(lon=lon, lat=lat)
+
+    # Determine rotation matrix to make it so that the circle is centered
+    # on the correct longitude/latitude.
+    m1 = rotation_matrix(-(0.5 * np.pi * u.radian - lat0), axis='y')
+    m2 = rotation_matrix(-lon0, axis='z')
+    transform_matrix = m2 * m1
+
+    # Apply 3D rotation
+    polygon = polygon.to_cartesian()
+    polygon = polygon.transform(transform_matrix)
+    polygon = UnitSphericalRepresentation.from_cartesian(polygon)
+
+    return polygon.lon, polygon.lat

regions/shapes/circle.py

@@ -3,12 +3,17 @@
 
 import math
 
-from astropy.coordinates import Angle
+import numpy as np
+
+from astropy import units as u
+from astropy.coordinates import Angle, SkyCoord
 from astropy.wcs.utils import pixel_to_skycoord
 
+from .polygon import PolygonSkyRegion
+
 from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
 from .._utils.wcs_helpers import skycoord_to_pixel_scale_angle
-from .._geometry import circular_overlap_grid
+from .._geometry import circular_overlap_grid, rotate_polygon
 
 __all__ = ['CirclePixelRegion', 'CircleSkyRegion']
 
@@ -153,7 +158,37 @@ def area(self):
     def contains(self, skycoord):
         return self.center.separation(skycoord) < self.radius
 
-    def to_pixel(self, wcs, mode='local', tolerance=None):
+    def to_polygon(self, points=100):
+        """
+        Convert the circle to a polygon.
+
+        Parameters
+        ----------
+        points : int, optional
+            Number of points in the final polygon.
+        """
+
+        # TODO: avoid converting to unit spherical or spherical if already
+        #       using a spherical representation
+
+        # Extract longitude/latitude, either from a tuple of two quantities, or
+        # a single 2-element Quantity.
+        rep = self.center.represent_as('unitspherical')
+        longitude, latitude = rep.lon, rep.lat
+
+        # Start off by generating the circle around the North pole
+        lon = np.linspace(0., 2 * np.pi, points + 1)[:-1] * u.radian
+        lat = np.repeat(0.5 * np.pi - self.radius.to(u.radian).value, points) * u.radian
+
+        # Now rotate it to the correct longitude/latitude
+        lon, lat = rotate_polygon(lon, lat, longitude, latitude)
+
+        # Make a new SkyCoord
+        vertices_sky = SkyCoord(lon, lat, frame=self.center)
+
+        return PolygonSkyRegion(vertices_sky)
+
+    def to_pixel(self, wcs, mode='local', tolerance=100):
         """
         Given a WCS, convert the circle to a best-approximation circle in pixel
         dimensions.
@@ -164,20 +199,25 @@ def to_pixel(self, wcs, mode='local', tolerance=None):
             A world coordinate system
         mode : 'local' or not
             not implemented
-        tolerance : None
+        tolerance : int
             not implemented
 
         Returns
         -------
         CirclePixelRegion
         """
 
-        if mode != 'local':
-            raise NotImplementedError
-        if tolerance is not None:
-            raise NotImplementedError
+        if mode == 'local':
+            center, scale, _ = skycoord_to_pixel_scale_angle(self.center, wcs)
+            radius = self.radius.to('deg').value * scale
+            return CirclePixelRegion(center, radius)
+
+        elif mode == 'affine':
+            raise NotImplementedError()
 
-        center, scale, _ = skycoord_to_pixel_scale_angle(self.center, wcs)
-        radius = self.radius.to('deg').value * scale
+        elif mode == 'full':
 
-        return CirclePixelRegion(center, radius)
+            return self.to_polygon(points=tolerance).to_pixel(wcs)
+
+        else:
+            raise ValueError('mode should be one of local/affine/full')

regions/shapes/polygon.py

@@ -4,7 +4,9 @@
 
 import numpy as np
 
-from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
+from astropy.wcs.utils import skycoord_to_pixel
+
+from ..core import PixelRegion, SkyRegion, Mask, BoundingBox, PixCoord
 from .._geometry import polygonal_overlap_grid
 from .._geometry.pnpoly import points_in_polygon
 
@@ -141,5 +143,11 @@ def contains(self, skycoord):
         raise NotImplementedError
 
     def to_pixel(self, wcs, mode='local', tolerance=None):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+
+        if mode == 'local':
+            x, y = skycoord_to_pixel(self.vertices, wcs)
+            vertices_pix = PixCoord(x, y)
+            return PolygonPixelRegion(vertices_pix)
+        else:
+
+            raise NotImplementedError()