Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Initial implementation of CircleSkyRegion.to_pixel in full mode #91

Merged
merged 4 commits into from
Dec 1, 2016
Merged
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
1 change: 1 addition & 0 deletions regions/_geometry/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 *
Expand Down
40 changes: 40 additions & 0 deletions regions/_geometry/rotate_polygon.py
Original file line number Diff line number Diff line change
@@ -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
62 changes: 51 additions & 11 deletions regions/shapes/circle.py
Original file line number Diff line number Diff line change
Expand Up @@ -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']

Expand Down Expand Up @@ -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.
Expand All @@ -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')
14 changes: 11 additions & 3 deletions regions/shapes/polygon.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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()