diff options
Diffstat (limited to 'lib/skyline.py')
| -rw-r--r-- | lib/skyline.py | 551 |
1 files changed, 0 insertions, 551 deletions
diff --git a/lib/skyline.py b/lib/skyline.py deleted file mode 100644 index eb58a2f..0000000 --- a/lib/skyline.py +++ /dev/null @@ -1,551 +0,0 @@ -# -*- coding: utf-8 -*- - -# Copyright (C) 2011 Association of Universities for Research in -# Astronomy (AURA) -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# 1. Redistributions of source code must retain the above -# copyright notice, this list of conditions and the following -# disclaimer. -# -# 2. Redistributions in binary form must reproduce the above -# copyright notice, this list of conditions and the following -# disclaimer in the documentation and/or other materials -# provided with the distribution. -# -# 3. The name of AURA and its representatives may not be used to -# endorse or promote products derived from this software without -# specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR -# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -""" -This module provides support for working with footprints -on the sky. Primary use case would use the following -generalized steps: - - #. Initialize `SkyLine` objects for each input image. - This object would be the union of all the input - image's individual chips WCS footprints. - - #. Determine overlap between all images. The - determination would employ a recursive operation - to return the extended list of all overlap values - computed as [img1 vs [img2,img3,...,imgN],img2 vs - [img3,...,imgN],...] - - #. Select the pair with the largest overlap, or the - pair which produces the largest overlap with the - first input image. This defines the initial - reference `SkyLine` object. - - #. Perform some operation on the 2 images: for example, - match sky in intersecting regions, or aligning - second image with the first (reference) image. - - #. Update the second image, either apply the sky value - or correct the WCS, then generate a new `SkyLine` - object for that image. - - #. Create a new reference `SkyLine` object as the union - of the initial reference object and the newly - updated `SkyLine` object. - - #. Repeat Steps 2-6 for all remaining input images. - -This process will work reasonably fast as most operations -are performed using the `SkyLine` objects and WCS information -solely, not image data itself. - -""" -from __future__ import division, print_function, absolute_import - -# STDLIB -from copy import copy, deepcopy - -# THIRD-PARTY -import pyfits -from stwcs import wcsutil -from stwcs.distortion.utils import output_wcs - -# LOCAL -from .polygon import SphericalPolygon - -# DEBUG -SKYLINE_DEBUG = True - -__all__ = ['SkyLineMember', 'SkyLine'] -__version__ = '0.4a' -__vdate__ = '12-Jul-2012' - -class SkyLineMember(object): - """ - Container for `SkyLine` members with these attributes: - - * `fname`: Image name (with path if given) - * `ext`: Tuple of extensions read - * `wcs`: `HSTWCS` object the composite data - * `polygon`: `~sphere.polygon.SphericalPolygon` object of the composite data - - """ - def __init__(self, fname, extname): - """ - Parameters - ---------- - fname : str - FITS image. - - extname : str - EXTNAME to use. SCI is recommended for normal - HST images. PRIMARY if image is single ext. - - """ - extname = extname.upper() - ext_list = [] - wcs_list = [] - - with pyfits.open(fname) as pf: - for i,ext in enumerate(pf): - if ext.name.upper() == extname: - ext_list.append(i) - wcs_list.append(wcsutil.HSTWCS(fname, ext=i)) - - # By combining WCS first before polygon, will remove chip gaps - - n_wcs = len(wcs_list) - if n_wcs > 1: - self._wcs = output_wcs(wcs_list) - elif n_wcs == 1: - self._wcs = wcs_list[0] - else: - raise ValueError('%s has no WCS' % fname) - - self._fname = fname - self._ext = tuple(ext_list) - self._polygon = SphericalPolygon.from_wcs(self.wcs) - - def __repr__(self): - return '%s(%r, %r, %r, %r)' % (self.__class__.__name__, self.fname, - self.ext, self.wcs, self.polygon) - - @property - def fname(self): - return self._fname - - @property - def ext(self): - return self._ext - - @property - def wcs(self): - return self._wcs - - @property - def polygon(self): - return self._polygon - -class SkyLine(object): - """ - Manage outlines on the sky. - - Each `SkyLine` has a list of `~SkyLine.members` and - a composite `~SkyLine.polygon` with all the - functionalities of `~sphere.polygon.SphericalPolygon`. - - """ - def __init__(self, fname, extname='SCI'): - """ - Parameters - ---------- - fname : str - FITS image. `None` to create empty `SkyLine`. - - extname : str - EXTNAME to use. SCI is recommended for normal - HST images. PRIMARY if image is single ext. - - """ - # Convert SCI data to SkyLineMember - if fname is not None: - self.members = [SkyLineMember(fname, extname)] - else: - self.members = [] - - # Put mosaic of all the chips in SkyLine - n = len(self.members) - if n == 0: - self.polygon = SphericalPolygon([]) - elif n == 1: - self.polygon = copy(self.members[0].polygon) - else: - raise ValueError('%s cannot initialize polygon with ' - 'multiple members' % self.__class__.__name__) - - def __getattr__(self, what): - """Control attribute access to `~sphere.polygon.SphericalPolygon`.""" - if what in ('from_radec', 'from_cone', 'from_wcs', - 'multi_union', 'multi_intersection', - '_find_new_inside',): - raise AttributeError('\'%s\' object has no attribute \'%s\'' % - (self.__class__.__name__, what)) - else: - return getattr(self.polygon, what) - - def __copy__(self): - return deepcopy(self) - - def __repr__(self): - return '%s(%r, %r)' % (self.__class__.__name__, - self.polygon, self.members) - - @property - def polygon(self): - """ - `~sphere.polygon.SphericalPolygon` portion of `SkyLine` - that contains the composite skyline from `members` - belonging to *self*. - - """ - return self._polygon - - @polygon.setter - def polygon(self, value): - assert isinstance(value, SphericalPolygon) - self._polygon = copy(value) # Deep copy - - @property - def members(self): - """ - List of `SkyLineMember` objects that belong to *self*. - Duplicate members are discarded. Members are kept in - the order of their additions to *self*. - - """ - return self._members - - @members.setter - def members(self, values): - self._members = [] - - # Not using set to preserve order - for v in values: - # Report corrupted members list instead of skipping - assert isinstance(v, SkyLineMember) - - if v not in self._members: - self._members.append(v) - - def _indv_mem_wcslist(self): - """List of original HSTWCS from each EXT in each member.""" - wcs_list = [] - - for m in self.members: - for i in m.ext: - wcs_list.append(wcsutil.HSTWCS(m.fname, ext=i)) - - return wcs_list - - def to_wcs(self): - """ - Combine `HSTWCS` objects from all `members` and return - a new `HSTWCS` object. If no `members`, return `None`. - - .. warning:: This cannot return WCS of intersection. - - """ - wcs_list = self._indv_mem_wcslist() - - n = len(wcs_list) - if n > 1: - wcs = output_wcs(wcs_list) - elif n == 1: - wcs = wcs_list[0] - else: - wcs = None - - return wcs - - def _rough_id(self): - """Filename of first member.""" - if len(self.members) > 0: - return self.members[0].fname - else: - return None - - def _draw_members(self, map, **kwargs): - """ - Draw individual extensions in members. - Useful for debugging. - - Parameters - ---------- - map : Basemap axes object - - **kwargs : Any plot arguments to pass to basemap - - """ - wcs_list = self._indv_mem_wcslist() - - for wcs in wcs_list: - poly = SphericalPolygon.from_wcs(wcs) - poly.draw(map, **kwargs) - - def _find_members(self, given_members): - """ - Find `SkyLineMember` in *given_members* that is in - *self*. This is used for intersection. - - Parameters - ---------- - self : obj - `SkyLine` instance. - - given_members : list - List of `SkyLineMember` to consider. - - Returns - ------- - new_members : list - List of `SkyLineMember` belonging to *self*. - - """ - if len(self.points) > 0: - out_mem = [m for m in given_members if - self.intersects_poly(m.polygon)] - else: - out_mem = [] - return out_mem - - def add_image(self, other): - """ - Return a new `SkyLine` that is the union of *self* - and *other*. - - .. warning:: `SkyLine.union` only returns `polygon` - without `members`. - - Parameters - ---------- - other : `SkyLine` object - - Examples - -------- - >>> s1 = SkyLine('image1.fits') - >>> s2 = SkyLine('image2.fits') - >>> s3 = s1.add_image(s2) - - """ - newcls = self.__class__(None) - newcls.polygon = self.union(other) - newcls.members = self.members + other.members - return newcls - - def find_intersection(self, other): - """ - Return a new `SkyLine` that is the intersection of - *self* and *other*. - - .. warning:: `SkyLine.intersection` only returns - `polygon` without `members`. - - Parameters - ---------- - other : `SkyLine` object - - Examples - -------- - >>> s1 = SkyLine('image1.fits') - >>> s2 = SkyLine('image2.fits') - >>> s3 = s1.find_intersection(s2) - - """ - newcls = self.__class__(None) - newcls.polygon = self.intersection(other) - newcls.members = newcls._find_members(self.members + other.members) - return newcls - - def find_max_overlap(self, skylines): - """ - Find `SkyLine` from a list of *skylines* that overlaps - the most with *self*. - - Parameters - ---------- - skylines : list - A list of `SkyLine` instances. - - Returns - ------- - max_skyline : `SkyLine` instance or `None` - `SkyLine` that overlaps the most or `None` if no - overlap found. This is *not* a copy. - - max_overlap_area : float - Area of intersection. - - """ - max_skyline = None - max_overlap_area = 0.0 - - for next_s in skylines: - try: - overlap_area = self.intersection(next_s).area() - except (ValueError, AssertionError): - if SKYLINE_DEBUG: - print('WARNING: Intersection failed for %s and %s. ' - 'Ignoring %s...' % (self._rough_id(), - next_s._rough_id(), - next_s._rough_id())) - overlap_area = 0.0 - else: - raise - - if overlap_area > max_overlap_area: - max_overlap_area = overlap_area - max_skyline = next_s - - return max_skyline, max_overlap_area - - @staticmethod - def max_overlap_pair(skylines): - """ - Find a pair of skylines with maximum overlap. - - Parameters - ---------- - skylines : list - A list of `SkyLine` instances. - - Returns - ------- - max_pair : tuple - Pair of `SkyLine` objects with max overlap - among given *skylines*. If no overlap found, - return `None`. These are *not* copies. - - """ - max_pair = None - max_overlap_area = 0.0 - - for i in xrange(len(skylines) - 1): - curr_s = skylines[i] - next_s, i_area = curr_s.find_max_overlap(skylines[i+1:]) - - if i_area > max_overlap_area: - max_overlap_area = i_area - max_pair = (curr_s, next_s) - - return max_pair - - @classmethod - def mosaic(cls, skylines, verbose=True): - """ - Mosaic all overlapping *skylines*. - - A pair of skylines with the most overlap is used as - a starting point. Then a skyline that overlaps the - most with the mosaic is used, and so forth until no - overlapping skyline is found. - - Parameters - ---------- - skylines : list - A list of `SkyLine` objects. - - verbose : bool - Print info to screen. - - Returns - ------- - mosaic : `SkyLine` instance or `None` - Union of all overlapping *skylines*, or `None` if - no overlap found. - - included : list - List of image names added to mosaic in the order - of addition. - - excluded : list - List of image names excluded because they do not - overlap with mosaic. - - """ - out_order = [] - excluded = [] - - if verbose: - print('***** SKYLINE MOSAIC *****') - - starting_pair = cls.max_overlap_pair(skylines) - if starting_pair is None: - if verbose: - print(' Cannot find any overlapping skylines. Aborting...') - return starting_pair, out_order, excluded - - remaining = list(skylines) - - s1, s2 = starting_pair - if verbose: - print(' Starting pair: %s, %s' % - (s1._rough_id(), s2._rough_id())) - - mosaic = s1.add_image(s2) - out_order = [s1._rough_id(), s2._rough_id()] - remaining.remove(s1) - remaining.remove(s2) - - while len(remaining) > 0: - next_skyline, i_area = mosaic.find_max_overlap(remaining) - - if next_skyline is None: - for r in remaining: - if verbose: - print(' No overlap: Excluding %s...' % r._rough_id()) - excluded.append(r._rough_id()) - break - - try: - new_mos = mosaic.add_image(next_skyline) - except (ValueError, AssertionError): - if SKYLINE_DEBUG: - print('WARNING: Cannot add %s to mosaic. Skipping it...' % - next_skyline._rough_id()) - excluded.append(next_skyline._rough_id()) - else: - raise - else: - print(' Adding %s to mosaic...' % next_skyline._rough_id()) - mosaic = new_mos - out_order.append(next_skyline._rough_id()) - finally: - remaining.remove(next_skyline) - - return mosaic, out_order, excluded - - @classmethod - def _find_frosty(cls, show=False): - s1 = SphericalPolygon.from_cone(0, 35, 8) - s2 = SphericalPolygon.from_cone(0, 20, 13) - s3 = SphericalPolygon.from_cone(0, 0, 20) - ss = SphericalPolygon.multi_union([s1,s2,s3]) - frosty = cls(None) - frosty.polygon = ss - if show: - from mpl_toolkits.basemap import Basemap - map = Basemap() - frosty.draw(map) - print('Frosty the Snowman says Brrr so cold...') - return frosty |