From 8982b54f7f7c37ca00489d9ed2e9f8ee468a3c8d Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 27 Jan 2022 17:08:33 +0000 Subject: [PATCH 1/5] Remove bundled pvextractor and updated extraction tool to work properly with change to APE-14 coords objects --- glue/external/pvextractor/LICENSE | 29 -- glue/external/pvextractor/__init__.py | 7 - .../external/pvextractor/geometry/__init__.py | 3 - glue/external/pvextractor/geometry/helpers.py | 97 ---- .../pvextractor/geometry/line_slices.py | 75 --- glue/external/pvextractor/geometry/path.py | 262 ----------- .../pvextractor/geometry/poly_slices.py | 68 --- glue/external/pvextractor/geometry/polygon.py | 37 -- glue/external/pvextractor/geometry/slices.py | 47 -- .../pvextractor/geometry/tests/__init__.py | 0 glue/external/pvextractor/gui.py | 430 ------------------ glue/external/pvextractor/pvextractor.py | 119 ----- glue/external/pvextractor/pvregions.py | 178 -------- glue/external/pvextractor/utils/__init__.py | 0 .../external/pvextractor/utils/wcs_slicing.py | 43 -- glue/external/pvextractor/utils/wcs_utils.py | 71 --- glue/plugins/tools/pv_slicer/qt/pv_slicer.py | 11 +- setup.cfg | 2 +- 18 files changed, 8 insertions(+), 1471 deletions(-) delete mode 100644 glue/external/pvextractor/LICENSE delete mode 100644 glue/external/pvextractor/__init__.py delete mode 100644 glue/external/pvextractor/geometry/__init__.py delete mode 100644 glue/external/pvextractor/geometry/helpers.py delete mode 100644 glue/external/pvextractor/geometry/line_slices.py delete mode 100644 glue/external/pvextractor/geometry/path.py delete mode 100644 glue/external/pvextractor/geometry/poly_slices.py delete mode 100644 glue/external/pvextractor/geometry/polygon.py delete mode 100644 glue/external/pvextractor/geometry/slices.py delete mode 100644 glue/external/pvextractor/geometry/tests/__init__.py delete mode 100644 glue/external/pvextractor/gui.py delete mode 100644 glue/external/pvextractor/pvextractor.py delete mode 100644 glue/external/pvextractor/pvregions.py delete mode 100644 glue/external/pvextractor/utils/__init__.py delete mode 100644 glue/external/pvextractor/utils/wcs_slicing.py delete mode 100644 glue/external/pvextractor/utils/wcs_utils.py diff --git a/glue/external/pvextractor/LICENSE b/glue/external/pvextractor/LICENSE deleted file mode 100644 index 70b032806..000000000 --- a/glue/external/pvextractor/LICENSE +++ /dev/null @@ -1,29 +0,0 @@ -Copyright (c) 2014, pvextractor developers -All rights reserved. - -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. Neither the name of the copyright holder nor the names of its -contributors may be used to endorse or promote products derived from -this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT -HOLDER OR CONTRIBUTORS 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. diff --git a/glue/external/pvextractor/__init__.py b/glue/external/pvextractor/__init__.py deleted file mode 100644 index a3a9699a8..000000000 --- a/glue/external/pvextractor/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -from . import utils -from .pvextractor import extract_pv_slice -from .utils.wcs_slicing import slice_wcs -from .geometry import Path -from .pvregions import paths_from_regfile diff --git a/glue/external/pvextractor/geometry/__init__.py b/glue/external/pvextractor/geometry/__init__.py deleted file mode 100644 index 5ad14c776..000000000 --- a/glue/external/pvextractor/geometry/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .slices import extract_slice -from .path import Path -from .helpers import PathFromCenter diff --git a/glue/external/pvextractor/geometry/helpers.py b/glue/external/pvextractor/geometry/helpers.py deleted file mode 100644 index 82b8b038a..000000000 --- a/glue/external/pvextractor/geometry/helpers.py +++ /dev/null @@ -1,97 +0,0 @@ -import numpy as np - -from astropy import units as u -from astropy.coordinates import (SkyCoord, BaseCoordinateFrame, - UnitSphericalRepresentation, - CartesianRepresentation) - -from .path import Path - - -class PathFromCenter(Path): - """ - A simple path defined by a center, length, and position angle. - - Parameters - ---------- - center : `~astropy.coordinates.SkyCoord` - The center of the path - length : `~astropy.units.Quantity` - The length of the path in angular units - angle : `~astropy.units.Quantity` - The position angle of the path, counter-clockwise - sample : int - How many points to sample along the path. By default, this is 2 (the - two end points. For small fields of view, this will be a good - approximation to the path, but for larger fields of view, where - spherical distortions become important, this should be increased to - provide a smooth path. - width : None or float or :class:`~astropy.units.Quantity` - The width of the path. If ``coords`` is passed as a list of pixel - positions, the width should be given (if passed) as a floating-point - value in pixels. If ``coords`` is a coordinate object, the width - should be passed as a :class:`~astropy.units.Quantity` instance with - units of angle. If None, interpolation is used at the position of the - path. - - Notes - ----- - The orientation of the final path will be such that for a position angle of - zero, the path goes from South to North. For a position angle of 90 - degrees, the path will go from West to East. - """ - - def __init__(self, center, length=None, angle=None, sample=2, width=None): - - # Check input types - - if not isinstance(center, (SkyCoord, BaseCoordinateFrame)): - raise TypeError("The central position should be given as a SkyCoord object") - - if not isinstance(length, u.Quantity) or not length.unit.is_equivalent(u.deg): - raise TypeError("The length should be given as an angular Quantity") - - if not isinstance(angle, u.Quantity) or not angle.unit.is_equivalent(u.deg): - raise TypeError("The angle should be given as an angular Quantity") - - # We set up the path by adding and removing half the length along the - # declination axis, then rotate the resulting two points around the - # center. - - # Convert the central position to cartesian coordinates - c1, c2, c3 = center.cartesian.xyz.value - - # Find the end points of the path - clon, clat = center.spherical.lon, center.spherical.lat - try: - plat = clat + np.linspace(-length * 0.5, length * 0.5, sample) - except ValueError: # Numpy 1.10+ - plat = clat + np.linspace(-length.value * 0.5, length.value * 0.5, sample) * length.unit - - x, y, z = UnitSphericalRepresentation(clon, plat).to_cartesian().xyz.value - - # Rotate around central point - - # Because longitude increases to the left, we have to take -angle - angle = -angle - - # We rotate (x,y,z) around (c1,c2,c3) by making use of the following - # equations: - # - # x' = x cos a + (1 - cos a)(c1c1x + c1c2y + c1c3z) + (c2z - c3y)sin a - # y' = y cos a + (1 - cos a)(c2c1x + c2c2y + c2c3z) + (c3x - c1z)sin a - # z' = z cos a + (1 - cos a)(c3c1x + c3c2y + c3c3z) + (c1y - c2x)sin a - # - # Source: https://www.uwgb.edu/dutchs/MATHALGO/sphere0.htm - - cosa = np.cos(angle) - sina = np.sin(angle) - - xd = x * cosa + (1 - cosa) * (c1*c1*x + c1*c2*y + c1*c3*z) + (c2 * z - c3 * y) * sina - yd = y * cosa + (1 - cosa) * (c2*c1*x + c2*c2*y + c2*c3*z) + (c3 * x - c1 * z) * sina - zd = z * cosa + (1 - cosa) * (c3*c1*x + c3*c2*y + c3*c3*z) + (c1 * y - c2 * x) * sina - - # Construct representations - points = center.realize_frame(CartesianRepresentation(x=xd, y=yd, z=zd)) - - super(PathFromCenter, self).__init__(points, width=width) diff --git a/glue/external/pvextractor/geometry/line_slices.py b/glue/external/pvextractor/geometry/line_slices.py deleted file mode 100644 index e5c9e59ca..000000000 --- a/glue/external/pvextractor/geometry/line_slices.py +++ /dev/null @@ -1,75 +0,0 @@ -from __future__ import print_function - -import numpy as np - -from scipy.ndimage import map_coordinates - - -def extract_line_slice(cube, x, y, order=3, respect_nan=True): - """ - Given an array with shape (z, y, x), extract a (z, n) slice by - interpolating at n (x, y) points. - - All units are in *pixels*. - - .. note:: If there are NaNs in the cube, they will be treated as zeros when - using spline interpolation. - - Parameters - ---------- - cube : `~numpy.ndarray` - The data cube to extract the slice from - curve : list or tuple - A list or tuple of (x, y) pairs, with minimum length 2 - order : int, optional - Spline interpolation order. Set to ``0`` for nearest-neighbor - interpolation. - - Returns - ------- - slice : `numpy.ndarray` - The (z, d) slice - """ - - if order == 0: - - total_slice = np.zeros([cube.shape[0], len(x)]) + np.nan - - x = np.round(x) - y = np.round(y) - - ok = (x >= 0) & (y >= 0) & (x < cube.shape[2]) & (y < cube.shape[1]) - - total_slice[:,ok] = cube[:, y[ok].astype(int), x[ok].astype(int)] - - elif order > 0 and order == int(order): - - nx = len(x) - nz = cube.shape[0] - - zi = np.outer(np.arange(nz, dtype=int), np.ones(nx)) - xi = np.outer(np.ones(nz), x) - yi = np.outer(np.ones(nz), y) - - if np.any(np.isnan(cube)): - - # map_coordinates does not deal well with NaN values so we have - # to remove the NaN values then re-mask the final slice. - - total_slice = map_coordinates(np.nan_to_num(cube), [zi,yi,xi], - order=order, cval=np.nan) - - slice_bad = map_coordinates(np.nan_to_num(np.isnan(cube).astype(int)), - [zi,yi,xi], order=order) - - total_slice[np.nonzero(slice_bad)] = np.nan - - else: - - total_slice = map_coordinates(cube, [zi,yi,xi], order=order, cval=np.nan) - - else: - - raise TypeError("order should be a positive integer") - - return total_slice diff --git a/glue/external/pvextractor/geometry/path.py b/glue/external/pvextractor/geometry/path.py deleted file mode 100644 index ee56f8ce2..000000000 --- a/glue/external/pvextractor/geometry/path.py +++ /dev/null @@ -1,262 +0,0 @@ -from __future__ import print_function - -import sys - -import numpy as np - -from astropy.wcs import WCSSUB_CELESTIAL -from astropy.wcs.utils import wcs_to_celestial_frame -from astropy.coordinates import BaseCoordinateFrame - -from ..utils.wcs_utils import get_spatial_scale - - -class Polygon(object): - def __init__(self, x, y): - self.x = x - self.y = y - - -def segment_angles(x, y): - - dx = np.diff(x) - dy = np.diff(y) - - d = np.hypot(dx, dy) - - cos_theta = (-dx[:-1] * dx[1:] - dy[:-1] * dy[1:]) / (d[:-1] * d[1:]) - cos_theta = np.clip(cos_theta, -1., 1.) - - sin_theta = (-dx[:-1] * dy[1:] + dy[:-1] * dx[1:]) / (d[:-1] * d[1:]) - sin_theta = np.clip(sin_theta, -1., 1.) - - theta = np.arctan2(sin_theta, cos_theta) - - theta[0] = np.pi - theta[-1] = np.pi - - return theta - - -def get_endpoints(x, y, width): - - # Pad with same values at ends, to find slope of perpendicular end - # lines. - try: - xp = np.pad(x, 1, mode='edge') - yp = np.pad(y, 1, mode='edge') - except AttributeError: # Numpy < 1.7 - xp = np.hstack([x[0], x, x[-1]]) - yp = np.hstack([y[0], y, y[-1]]) - - dx = np.diff(xp) - dy = np.diff(yp) - - alpha = segment_angles(xp, yp) / 2. - beta = np.arctan2(dy, dx)[:-1] - beta[0] = beta[1] - gamma = -(np.pi - alpha - beta) - - dx = np.cos(gamma) - dy = np.sin(gamma) - - angles = segment_angles(xp, yp) / 2. - - # Find points offset from main curve, on bisecting lines - x1 = x - dx * width * 0.5 / np.sin(angles) - x2 = x + dx * width * 0.5 / np.sin(angles) - y1 = y - dy * width * 0.5 / np.sin(angles) - y2 = y + dy * width * 0.5 / np.sin(angles) - - return x1, y1, x2, y2 - - -class Path(object): - """ - A curved path that may have a non-zero width and is used to extract - slices from cubes. - - Parameters - ---------- - xy_or_coords : list or Astropy coordinates - The points defining the path. This can be passed as a list of (x, y) - tuples, which is interpreted as being pixel positions, or it can be - an Astropy coordinate object containing an array of 2 or more - coordinates. - width : None or float or :class:`~astropy.units.Quantity` - The width of the path. If ``coords`` is passed as a list of pixel - positions, the width should be given (if passed) as a floating-point - value in pixels. If ``coords`` is a coordinate object, the width - should be passed as a :class:`~astropy.units.Quantity` instance with - units of angle. If None, interpolation is used at the position of the - path. - """ - - def __init__(self, xy_or_coords, width=None): - if isinstance(xy_or_coords, list): - self._xy = xy_or_coords - self._coords = None - elif sys.version_info[0] > 2 and isinstance(xy_or_coords, zip): - self._xy = list(xy_or_coords) - self._coords = None - else: - self._xy = None - self._coords = xy_or_coords - self.width = width - - def add_point(self, xy_or_coord): - """ - Add a point to the path - - Parameters - ---------- - xy_or_coord : tuple or Astropy coordinate - A tuple (x, y) containing the coordinates of the point to add (if - the path is defined in pixel space), or an Astropy coordinate - object (if it is defined in world coordinates). - """ - if self._xy is not None: - if isinstance(xy_or_coord, tuple): - self._xy.append(xy_or_coord) - else: - raise TypeError("Path is defined as a list of pixel " - "coordinates, so `xy_or_coord` should be " - "a tuple of `(x,y)` pixel coordinates.") - else: - if isinstance(xy_or_coord, BaseCoordinateFrame): - raise NotImplementedError("Cannot yet append world coordinates to path") - else: - raise TypeError("Path is defined in world coordinates, " - "so `xy_or_coord` should be an Astropy " - "coordinate object.") - - def get_xy(self, wcs=None): - """ - Return the pixel coordinates of the path. - - If the path is defined in world coordinates, the appropriate WCS - transformation should be passed. - - Parameters - ---------- - wcs : :class:`~astropy.wcs.WCS` - The WCS transformation to assume in order to transform the path - to pixel coordinates. - """ - if self._xy is not None: - return self._xy - else: - if wcs is None: - raise ValueError("`wcs` is needed in order to compute " - "the pixel coordinates") - else: - - # Extract the celestial component of the WCS - wcs_sky = wcs.sub([WCSSUB_CELESTIAL]) - - # Find the astropy name for the coordinates - celestial_system = wcs_to_celestial_frame(wcs_sky) - - world_coords = self._coords.transform_to(celestial_system) - - xw, yw = world_coords.spherical.lon.degree, world_coords.spherical.lat.degree - - return list(zip(*wcs_sky.wcs_world2pix(xw, yw, 0))) - - def sample_points_edges(self, spacing, wcs=None): - - x, y = zip(*self.get_xy(wcs=wcs)) - - # Find the distance interval between all pairs of points - dx = np.diff(x) - dy = np.diff(y) - dd = np.hypot(dx, dy) - - # Find the total displacement along the broken curve - d = np.hstack([0., np.cumsum(dd)]) - - # Figure out the number of points to sample, and stop short of the - # last point. - n_points = int(np.floor(d[-1] / spacing)) - - if n_points == 0: - raise ValueError("Path is shorter than spacing") - - d_sampled = np.linspace(0., n_points * spacing, n_points + 1) - - x_sampled = np.interp(d_sampled, d, x) - y_sampled = np.interp(d_sampled, d, y) - - return d_sampled, x_sampled, y_sampled - - def sample_points(self, spacing, wcs=None): - - d_sampled, x_sampled, y_sampled = self.sample_points_edges(spacing, wcs=wcs) - - x_sampled = 0.5 * (x_sampled[:-1] + x_sampled[1:]) - y_sampled = 0.5 * (y_sampled[:-1] + y_sampled[1:]) - - return x_sampled, y_sampled - - def sample_polygons(self, spacing, wcs=None): - - x, y = zip(*self.get_xy(wcs=wcs)) - - d_sampled, x_sampled, y_sampled = self.sample_points_edges(spacing, wcs=wcs) - - # Find the distance interval between all pairs of points - dx = np.diff(x) - dy = np.diff(y) - dd = np.hypot(dx, dy) - - # Normalize to find unit vectors - dx = dx / dd - dy = dy / dd - - # Find the total displacement along the broken curve - d = np.hstack([0., np.cumsum(dd)]) - - interval = np.searchsorted(d, d_sampled) - 1 - interval[0] = 0 - - dx = dx[interval] - dy = dy[interval] - - polygons = [] - - x_beg = x_sampled - x_end = x_sampled + dx * spacing - - y_beg = y_sampled - y_end = y_sampled + dy * spacing - - if hasattr(self.width, 'unit'): - scale = get_spatial_scale(wcs) - width = (self.width / scale).decompose().value - else: - width = self.width - - x1 = x_beg - dy * width * 0.5 - y1 = y_beg + dx * width * 0.5 - - x2 = x_end - dy * width * 0.5 - y2 = y_end + dx * width * 0.5 - - x3 = x_end + dy * width * 0.5 - y3 = y_end - dx * width * 0.5 - - x4 = x_beg + dy * width * 0.5 - y4 = y_beg - dx * width * 0.5 - - for i in range(len(x_sampled) - 1): - p = Polygon([x1[i], x2[i], x3[i], x4[i]], [y1[i], y2[i], y3[i], y4[i]]) - polygons.append(p) - - return polygons - - def to_patches(self, spacing, wcs=None, **kwargs): - from matplotlib.patches import Polygon as MPLPolygon - patches = [] - for poly in self.sample_polygons(spacing, wcs=wcs): - patches.append(MPLPolygon(list(zip(poly.x, poly.y)), **kwargs)) - return patches diff --git a/glue/external/pvextractor/geometry/poly_slices.py b/glue/external/pvextractor/geometry/poly_slices.py deleted file mode 100644 index c0129d415..000000000 --- a/glue/external/pvextractor/geometry/poly_slices.py +++ /dev/null @@ -1,68 +0,0 @@ -from __future__ import print_function - -import numpy as np -from astropy.utils.console import ProgressBar - -from .polygon import square_polygon_overlap_area - - -def extract_poly_slice(cube, polygons, return_area=False): - """ - Extract the values of polygonal chunks from a data cube - - Parameters - ---------- - cube : np.ndarray - polygons : - return_area : bool - If set, return the area of each polygon and the sum over that area. - Otherwise, return the mean. - """ - - nx = len(polygons) - nz = cube.shape[0] - - total_slice = np.zeros((nz, nx)) - total_area = np.zeros((nz, nx)) - - p = ProgressBar(len(polygons)) - - for i, polygon in enumerate(polygons): - - p.update() - - # Find bounding box - bbxmin = int(round(np.min(polygon.x))-1) - bbxmax = int(round(np.max(polygon.x))+2) - bbymin = int(round(np.min(polygon.y))-1) - bbymax = int(round(np.max(polygon.y))+2) - - # Clip to cube box - bbxmin = max(bbxmin, 0) - bbxmax = min(bbxmax, cube.shape[2]) - bbymin = max(bbymin, 0) - bbymax = min(bbymax, cube.shape[1]) - - # Loop through pixels that might overlap - for xmin in np.arange(bbxmin, bbxmax): - for ymin in np.arange(bbymin, bbymax): - - area = square_polygon_overlap_area(xmin-0.5, xmin+0.5, - ymin-0.5, ymin+0.5, - polygon.x, polygon.y) - - if area > 0: - dataslice = cube[:, ymin, xmin] - good_values = np.isfinite(dataslice) - if np.any(good_values): - total_slice[good_values, i] += dataslice[good_values] * area - total_area[good_values, i] += area - - total_slice[total_area == 0.] = np.nan - if return_area: - return total_slice, total_area - total_slice[total_area > 0.] /= total_area[total_area > 0.] - - print("") - - return total_slice diff --git a/glue/external/pvextractor/geometry/polygon.py b/glue/external/pvextractor/geometry/polygon.py deleted file mode 100644 index c5db660c9..000000000 --- a/glue/external/pvextractor/geometry/polygon.py +++ /dev/null @@ -1,37 +0,0 @@ -""" -This module implements polygon-square intersection using matplotlib. It is -twice as fast as Shapely for this specific case and avoids requiring another -dependency. -""" - -import numpy as np - -from matplotlib.path import Path -from matplotlib.transforms import Bbox - - -def square_polygon_intersection(xmin, xmax, ymin, ymax, x, y): - poly = Path(list(zip(x, y))) - box = Bbox([[xmin, ymin], [xmax, ymax]]) - try: - clipped_poly = poly.clip_to_bbox(box) - except ValueError: - return [], [] - else: - return clipped_poly.vertices[:, 0], clipped_poly.vertices[:, 1] - - -def polygon_area(x, y): - x1 = x - x2 = np.roll(x, -1) - y1 = y - y2 = np.roll(y, -1) - return abs(0.5 * np.sum(x1 * y2 - x2 * y1)) - - -def square_polygon_overlap_area(xmin, xmax, ymin, ymax, x, y): - x, y = square_polygon_intersection(xmin, xmax, ymin, ymax, x, y) - if len(x) == 0: - return 0. - else: - return polygon_area(x, y) diff --git a/glue/external/pvextractor/geometry/slices.py b/glue/external/pvextractor/geometry/slices.py deleted file mode 100644 index a918a0ed7..000000000 --- a/glue/external/pvextractor/geometry/slices.py +++ /dev/null @@ -1,47 +0,0 @@ -import numpy as np - -from .line_slices import extract_line_slice -from .poly_slices import extract_poly_slice - - -def extract_slice(cube, path, spacing=1.0, order=3, respect_nan=True, - wcs=None): - """ - Given an array with shape (z, y, x), extract a (z, n) slice from a path - with ``n`` segments. - - All units are in *pixels* - - .. note:: If there are NaNs in the cube, they will be treated as zeros when - using spline interpolation. - - Parameters - ---------- - path : `Path` - The path along which to define the slice - spacing : float - The position resolution in the final slice - order : int, optional - Spline interpolation order when using line paths. Does not have any - effect for polygon paths. - respect_nan : bool, optional - If set to `False`, NaN values are changed to zero before computing - the slices. - - Returns - ------- - slice : `numpy.ndarray` - The slice - """ - - if not respect_nan: - cube = np.nan_to_num(cube) - - if path.width is None: - x, y = path.sample_points(spacing=spacing, wcs=wcs) - slice = extract_line_slice(cube, x, y, order=order) - else: - polygons = path.sample_polygons(spacing=spacing, wcs=wcs) - slice = extract_poly_slice(cube, polygons) - - return slice diff --git a/glue/external/pvextractor/geometry/tests/__init__.py b/glue/external/pvextractor/geometry/tests/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/glue/external/pvextractor/gui.py b/glue/external/pvextractor/gui.py deleted file mode 100644 index 991132882..000000000 --- a/glue/external/pvextractor/gui.py +++ /dev/null @@ -1,430 +0,0 @@ -from __future__ import print_function - -import os -import math -import warnings - -import numpy as np - -from matplotlib.collections import LineCollection -from matplotlib.transforms import Bbox -from matplotlib.patches import Polygon - -from .geometry.path import Path, get_endpoints -from . import extract_pv_slice - - -def distance(x1, y1, x2, y2, x3, y3): - """ - Find the shortest distance between a point (x3, y3) and the line passing - through the points (x1, y1) and (x2, y2). - """ - - px = x2-x1 - py = y2-y1 - - something = px * px + py * py - - u = ((x3 - x1) * px + (y3 - y1) * py) / float(something) - - x = x1 + u * px - y = y1 + u * py - - dx = x - x3 - dy = y - y3 - - dist = math.sqrt(dx*dx + dy*dy) - - return dist - - -class MovableSliceBox(object): - - def __init__(self, box, callback): - self.box = box - self.press = None - self.background = None - self.point_counter = 0 - self.callback = callback - self.mode = 0 - self.show_poly = False - self.cidpress = self.box.figure.canvas.mpl_connect('draw_event', self.draw_slicer) - - def connect(self): - self.cidpress = self.box.figure.canvas.mpl_connect('key_press_event', self.key_press) - self.cidpress = self.box.figure.canvas.mpl_connect('button_press_event', self.on_press) - self.cidmotion = self.box.figure.canvas.mpl_connect('motion_notify_event', self.on_motion) - - def draw_slicer(self, event): - - axes = self.box.axes - canvas = self.box.figure.canvas - - self.box.axes.draw_artist(self.box) - - if self.show_poly: - - path = Path(zip(self.box.x, self.box.y)) - path.width = self.box.width - - patches = path.to_patches(1, ec='green', fc='none', - transform=self.box.axes.transData, - clip_on=True, clip_box=self.box.axes.bbox) - - for patch in patches: - self.box.axes.draw_artist(patch) - - def on_press(self, event): - - if self.box.figure.canvas.toolbar.mode != '': - return - - if event.inaxes != self.box.axes: - return - - if self.mode == 1: - self.callback(self.box) - self.mode += 1 - return - - if self.mode == 2: - self.box.x = [] - self.box.y = [] - self.mode = 0 - self.point_counter = 0 - - self.press = event.xdata, event.ydata - - self.point_counter += 1 - - axes = self.box.axes - canvas = self.box.figure.canvas - - if self.point_counter == 1: # first point - - self.box.x.append(event.xdata) - self.box.x.append(event.xdata) - self.box.y.append(event.ydata) - self.box.y.append(event.ydata) - - self.box.width = 0. - - self.box.set_animated(True) - canvas.draw() - self.background = canvas.copy_from_bbox(self.box.axes.bbox) - - elif self.mode == 0: - - self.box.x.append(event.xdata) - self.box.y.append(event.ydata) - - self.box._update_segments() - - # now redraw just the lineangle - axes.draw_artist(self.box) - - def key_press(self, event): - - if self.box.figure.canvas.toolbar.mode != '': - return - - if event.key == 'enter' and self.mode == 0: - self.mode += 1 - self.box.x = self.box.x[:-1] - self.box.y = self.box.y[:-1] - - if event.key == 'y' and self.mode == 2: - self.show_poly = not self.show_poly - self.draw_slicer(event) - self.box.figure.canvas.draw() - - def on_motion(self, event): - - if self.box.figure.canvas.toolbar.mode != '': - return - - if self.point_counter == 0: - return - - if self.mode == 2: - return - - canvas = self.box.figure.canvas - axes = self.box.axes - canvas.restore_region(self.background) - - if event.inaxes != axes: - return - - if self.mode == 0: - self.box.x[-1] = event.xdata - self.box.y[-1] = event.ydata - elif self.mode == 1: - self.box.width = distance(self.box.x[-2], self.box.y[-2], self.box.x[-1], self.box.y[-1], event.xdata, event.ydata) * 2 - - self.box._update_segments() - - # redraw just the current lineangle - axes.draw_artist(self.box) - - canvas.blit() - - def disconnect(self): - self.box.figure.canvas.mpl_disconnect(self.cidpress) - self.box.figure.canvas.mpl_disconnect(self.cidmotion) - - -class SliceCurve(LineCollection): - - def __init__(self, x=[], y=[], width=None, **kwargs): - - super(SliceCurve, self).__init__([], **kwargs) - - self.x = x - self.y = y - self.width = width - - self._update_segments() - - def _update_segments(self): - - if not self.x or self.width is None or len(self.x) < 2: - return - - # Find central line - line = zip(self.x, self.y) - - if self.width: - - x1, y1, x2, y2 = get_endpoints(self.x, self.y, self.width) - - # Find bounding rectangle - rect = zip(np.hstack([x1,x2[::-1], x1[0]]), - np.hstack([y1,y2[::-1], y1[0]])) - - self.set_segments([list(line), list(rect)]) - self.set_linestyles(['solid', 'dashed']) - self.set_linewidths([2, 1]) - - else: - - self.set_segments([list(line)]) - self.set_linestyles(['solid']) - self.set_linewidths([2,]) - - -def unitless(x): - if hasattr(x, 'unit'): - return x.value - else: - return x - - -class PVSlicer(object): - - def __init__(self, filename_or_cube, backend=None, clim=None, cmap=None): - - try: - from spectral_cube import SpectralCube - if isinstance(filename_or_cube, SpectralCube): - cube = filename_or_cube - else: - cube = SpectralCube.read(filename_or_cube, format='fits') - self.cube = cube - self.array = self.cube - self.shape = cube.shape - except ImportError: - warnings.warn("spectral_cube package is not available - using astropy.io.fits directly") - from astropy.io import fits - self.array = fits.getdata(filename_or_cube) - self.shape = array.shape - self.ok_mask = np.isfinite(self.array) - if self.array.ndim != 3: - raise ValueError("dataset does not have 3 dimensions (install the spectral_cube package to avoid this error)") - - import matplotlib as mpl - - # We reset the rc parameters to the default values to make sure we use - # the default interactive backend and to make sure that the UI is - # consistent. - mpl.rcdefaults() - - if backend is not None: - mpl.use(backend) - - import matplotlib.pyplot as plt - - self.fig = plt.figure(figsize=(8, 5)) - - self.backend = mpl.get_backend() - print("Using Matplotlib backend: {0}".format(self.backend)) - - self.cmap = cmap - - self.ax1 = self.fig.add_axes([0.1, 0.1, 0.4, 0.7], aspect='equal', adjustable='datalim') - - if clim is None: - warnings.warn("clim not defined and will be determined from the data") - # To work with large arrays, sub-sample the data - # (but don't do it for small arrays) - n1 = int(np.round(max(self.shape[0] / 10, 1))) - n2 = int(np.round(max(self.shape[1] / 10, 1))) - n3 = int(np.round(max(self.shape[2] / 10, 1))) - if hasattr(self,'cube'): - sub_cube = self.cube[::n1,::n2,::n3] - cmin = sub_cube.min().value - cmax = sub_cube.max().value - else: - sub_array = self.array[::n1,::n2,::n3] - sub_mask = self.ok_mask[::n1,::n2,::n3] - cmin = sub_array[sub_mask].min() - cmax = sub_array[sub_mask].max() - crange = cmax - cmin - self._clim = (cmin - crange, cmax + crange) - else: - self._clim = clim - - self.slice = int(round(self.shape[0] / 2.)) - - from matplotlib.widgets import Slider - - self.slice_slider_ax = self.fig.add_axes([0.1, 0.95, 0.4, 0.03]) - self.slice_slider_ax.set_xticklabels("") - self.slice_slider_ax.set_yticklabels("") - self.slice_slider = Slider(self.slice_slider_ax, "3-d slice", 0, self.shape[0]-1, valinit=self.slice, valfmt="%i") - self.slice_slider.on_changed(self.update_slice) - self.slice_slider.drawon = False - - self.image = self.ax1.imshow(unitless(self.array[self.slice, :,:]), - origin='lower', interpolation='nearest', - vmin=self._clim[0], vmax=self._clim[1], - cmap=self.cmap) - - self.vmin_slider_ax = self.fig.add_axes([0.1, 0.90, 0.4, 0.03]) - self.vmin_slider_ax.set_xticklabels("") - self.vmin_slider_ax.set_yticklabels("") - self.vmin_slider = Slider(self.vmin_slider_ax, "vmin", self._clim[0], self._clim[1], valinit=self._clim[0]) - self.vmin_slider.on_changed(self.update_vmin) - self.vmin_slider.drawon = False - - self.vmax_slider_ax = self.fig.add_axes([0.1, 0.85, 0.4, 0.03]) - self.vmax_slider_ax.set_xticklabels("") - self.vmax_slider_ax.set_yticklabels("") - self.vmax_slider = Slider(self.vmax_slider_ax, "vmax", self._clim[0], self._clim[1], valinit=self._clim[1]) - self.vmax_slider.on_changed(self.update_vmax) - self.vmax_slider.drawon = False - - self.grid1 = None - self.grid2 = None - self.grid3 = None - - self.ax2 = self.fig.add_axes([0.55, 0.1, 0.4, 0.7]) - - # Add slicing box - self.box = SliceCurve(colors=(0.8, 0.0, 0.0)) - self.ax1.add_collection(self.box) - self.movable = MovableSliceBox(self.box, callback=self.update_pv_slice) - self.movable.connect() - - # Add save button - from matplotlib.widgets import Button - self.save_button_ax = self.fig.add_axes([0.65, 0.90, 0.20, 0.05]) - self.save_button = Button(self.save_button_ax, 'Save slice to FITS') - self.save_button.on_clicked(self.save_fits) - self.file_status_text = self.fig.text(0.75, 0.875, "", ha='center', va='center') - self.set_file_status(None) - - self.set_file_status(None) - self.pv_slice = None - - self.cidpress = self.fig.canvas.mpl_connect('button_press_event', self.click) - - def set_file_status(self, status, filename=None): - if status == 'instructions': - self.file_status_text.set_text('Please enter filename in terminal') - self.file_status_text.set_color('red') - elif status == 'saved': - self.file_status_text.set_text('File successfully saved to {0}'.format(filename)) - self.file_status_text.set_color('green') - else: - self.file_status_text.set_text('') - self.file_status_text.set_color('black') - self.fig.canvas.draw() - - def click(self, event): - - if event.inaxes != self.ax2: - return - - self.slice_slider.set_val(event.ydata) - - def save_fits(self, *args, **kwargs): - - if self.pv_slice is None: - return - - # When using Qt, we need to use a proper Qt save dialog since input() - # does not play nicely with the Qt event loop. - if self.backend.lower().startswith('qt'): - from qtpy.compat import getsavefilename - plot_name, _ = getsavefilename() - if not plot_name: - return - else: - self.set_file_status('instructions') - - print("Enter filename: ", end='') - plot_name = str(input()) - - from astropy.io import fits - self.pv_slice.writeto(plot_name, overwrite=True) - print("Saved file to: ", plot_name) - - self.set_file_status('saved', filename=plot_name) - - def update_pv_slice(self, box): - - path = Path(zip(box.x, box.y)) - path.width = box.width - - self.pv_slice = extract_pv_slice(self.array, path) - - self.ax2.cla() - self.ax2.imshow(self.pv_slice.data, origin='lower', aspect='auto', - interpolation='nearest', cmap=self.cmap) - - self.fig.canvas.draw() - - def show(self, block=True): - import matplotlib.pyplot as plt - plt.show(block=block) - - def update_slice(self, pos=None): - - if self.array.ndim == 2: - self.image.set_array(unitless(self.array)) - else: - self.slice = int(round(pos)) - self.image.set_array(unitless(self.array[self.slice, :, :])) - - self.fig.canvas.draw() - - def update_vmin(self, vmin): - if vmin > self._clim[1]: - self._clim = (self._clim[1], self._clim[1]) - else: - self._clim = (vmin, self._clim[1]) - self.image.set_clim(*self._clim) - self.fig.canvas.draw() - - def update_vmax(self, vmax): - if vmax < self._clim[0]: - self._clim = (self._clim[0], self._clim[0]) - else: - self._clim = (self._clim[0], vmax) - self.image.set_clim(*self._clim) - self.fig.canvas.draw() - - def close(self): - import matplotlib.pyplot as plt - plt.close(self.fig) diff --git a/glue/external/pvextractor/pvextractor.py b/glue/external/pvextractor/pvextractor.py deleted file mode 100644 index 1fa39e662..000000000 --- a/glue/external/pvextractor/pvextractor.py +++ /dev/null @@ -1,119 +0,0 @@ -from __future__ import print_function - -import numpy as np -import six -import warnings - -from astropy import units as u -from astropy.io.fits import PrimaryHDU, ImageHDU, Header - -from .utils.wcs_utils import get_spatial_scale, sanitize_wcs -from .geometry import extract_slice -from .geometry import path as paths -from .utils.wcs_slicing import slice_wcs - - -def extract_pv_slice(cube, path, wcs=None, spacing=1.0, order=3, - respect_nan=True): - """ - Given a position-position-velocity cube with dimensions (nv, ny, nx), and - a path, extract a position-velocity slice. - - Alternative implementations: - gipsy::sliceview - karma::kpvslice - casaviewer::slice - - Parameters - ---------- - cube : :class:`~numpy.ndarray` or :class:`~spectral_cube.SpectralCube` or str or HDU - The cube to extract a slice from. If this is a plain - :class:`~numpy.ndarray` instance, the WCS information can optionally - be specified with the ``wcs`` parameter. If a string, it should be - the name of a file containing a spectral cube. - path : `Path` or list of 2-tuples - The path along which to define the position-velocity slice. The path - can contain coordinates defined in pixel or world coordinates. - wcs : :class:`~astropy.wcs.WCS`, optional - The WCS information to use for the cube. This should only be - specified if the ``cube`` parameter is a plain - :class:`~numpy.ndarray` instance. - spacing : float - The position resolution in the final position-velocity slice. This - can be given in pixel coordinates or as a - :class:`~astropy.units.Quantity` instance with angle units. - order : int, optional - Spline interpolation order when using paths with zero width. Does not - have any effect for paths with a non-zero width. - respect_nan : bool, optional - If set to `False`, NaN values are changed to zero before computing - the slices. If set to `True`, in the case of line paths a second - computation is performed to ignore the NaN value while interpolating, - and set the output values of NaNs to NaN. - - Returns - ------- - slice : `PrimaryHDU` - The position-velocity slice, as a FITS HDU object - """ - - if isinstance(cube, (str, ImageHDU, PrimaryHDU)): - try: - from spectral_cube import SpectralCube - cube = SpectralCube.read(cube) - except ImportError: - raise ImportError("spectral_cube package required for working " - "with fits data. Install spectral_cube or " - "use NumPy arrays") - - if _is_spectral_cube(cube): - wcs = cube.wcs - # The fits HEADER will preserve the UNIT, but pvextractor does not care - # what the flux units are - cube = cube.filled_data[...].value - - if wcs is not None: - wcs = sanitize_wcs(wcs) - - if not isinstance(cube, np.ndarray) or wcs is not None: - scale = get_spatial_scale(wcs) - if isinstance(spacing, u.Quantity): - pixel_spacing = (spacing / scale).decompose() - world_spacing = spacing - else: - pixel_spacing = spacing - world_spacing = spacing * scale - else: - if isinstance(spacing, u.Quantity): - raise TypeError("No WCS has been specified, so spacing should be given in pixels") - else: - pixel_spacing = spacing - world_spacing = None - - # Allow path to be passed in as list of 2-tuples - if not isinstance(path, paths.Path): - path = paths.Path(path) - - pv_slice = extract_slice(cube, path, wcs=wcs, spacing=pixel_spacing, - order=order, respect_nan=respect_nan) - - # Generate output header - if wcs is None: - header = Header() - else: - header = slice_wcs(wcs, spatial_scale=world_spacing).to_header() - - # TODO: write path to BinTableHDU - - return PrimaryHDU(data=pv_slice, header=header) - - -def _is_spectral_cube(obj): - try: - from spectral_cube.spectral_cube import BaseSpectralCube - return isinstance(obj, BaseSpectralCube) - except ImportError: - if 'SpectralCube' in str(obj.__class__): - warnings.warn("Object appears to be a SpectralCube, but" - " the spectral_cube module could not be loaded") - return False diff --git a/glue/external/pvextractor/pvregions.py b/glue/external/pvextractor/pvregions.py deleted file mode 100644 index c94d4a304..000000000 --- a/glue/external/pvextractor/pvregions.py +++ /dev/null @@ -1,178 +0,0 @@ -import numpy as np -from .geometry import path -from astropy import coordinates -from astropy import units as u -import re - -csystems = {'galactic':coordinates.Galactic, - 'fk5':coordinates.FK5, - 'fk4':coordinates.FK4, - 'icrs':coordinates.ICRS} -cel_systems = ['fk5','fk4','icrs'] -# ecliptic, detector, etc. not supported (because I don't know what they mean) -# (or with ecliptic, how to deal with them) -all_systems = cel_systems+['galactic','image','physical'] - -class SimpleRegion(object): - def __init__(self, coord_list, coord_format, name): - self.name = name - self.coord_format = coord_format - self.coord_list = coord_list - - def __repr__(self): - return "Region: {0}, {1}, {2}".format(self.name, self.coord_list, - self.coord_format) - - -valid_regions = ['line', 'segment', 'vector'] -valid_region_re = [re.compile("^"+n) for n in valid_regions] - -def simple_region_parser(regionstring, coord_format): - rs = regionstring.lstrip("# ") - - rtype = None - for rt, rre in zip(valid_regions, valid_region_re): - if rre.search(rs): - rtype = rt - break - - if rtype is None: - # not a usable region - return - - coordre = re.compile("^[a-z]*\((.*)\)") - coord_list = coordre.findall(rs) - if len(coord_list) != 1: - raise ValueError("Invalid region") - - coords = coord_list[0].split(",") - - outcoords = [] - for ii,cs in enumerate(coords): - if coord_format in csystems: - if ":" in cs: - # sexagesimal - if coord_format in cel_systems and ii % 2 == 0: - # odd, celestial = RA = hours - crd = coordinates.Angle(cs, unit=u.hour) - else: - crd = coordinates.Angle(cs, unit=u.deg) - else: - try: - # if it's a float, it's in degrees - crd = float(cs) * u.deg - except ValueError: - crd = coordinates.Angle(cs) - else: - # assume pixel units - crd = float(cs) - outcoords.append(crd) - - reg = SimpleRegion(coord_list=outcoords, coord_format=coord_format, - name=rtype) - - return reg - -def load_regions_file(rfile): - with open(rfile,'r') as fh: - lines = fh.readlines() - return load_regions_stringlist(lines) - -def load_regions_stringlist(lines): - - coord_format = None - for line in lines: - if line.strip() in all_systems: - coord_format = line.strip() - break - if coord_format is None: - raise ValueError("No valid coordinate format found.") - - regions_ = [simple_region_parser(line, coord_format) for line in lines] - regions = [r for r in regions_ if r is not None] - - return regions - - -def line_to_path(region): - """ - Convert a line or segment to a path - """ - - l,b = None,None - - endpoints = [] - - for x in region.coord_list: - if l is None: - if hasattr(x,'unit'): - l = x.to(u.deg).value - else: - l = x - else: - if hasattr(x,'unit'): - b = x.to(u.deg).value - else: - b = x - if l is not None and b is not None: - if hasattr(b,'unit') or hasattr(l,'unit'): - raise TypeError("Can't work with a list of quantities") - endpoints.append((l,b)) - l,b = None,None - else: - raise ValueError("unmatched l,b") - - lbarr = np.array(endpoints) - C = csystems[region.coord_format](lbarr[:,0]*u.deg, lbarr[:,1]*u.deg) - - # TODO: add widths for projection - - p = path.Path(C) - - return p - -def vector_to_path(vector_region): - """ - Convert a vector region to a path - - # vector(48.944348,-0.36432694,485.647",124.082) vector=1 - """ - - x,y = vector_region.coord_list[:2] - length = vector_region.coord_list[2] - angle = vector_region.coord_list[3] - - C1 = csystems[vector_region.coord_format](x, y) - dx,dy = length * np.cos(angle), length * np.sin(angle) - # -dx because we're in the flippy coordsys - C2 = csystems[vector_region.coord_format](C1.spherical.lon - dx, C1.spherical.lat + dy) - - C = csystems[vector_region.coord_format]([C1.spherical.lon.deg,C2.spherical.lon.deg]*u.deg, - [C1.spherical.lat.deg,C2.spherical.lat.deg]*u.deg) - - p = path.Path(C) - - return p - -region_converters = {'line':line_to_path, 'segment':line_to_path, - 'vector':vector_to_path} - -def paths_from_regfile(regfile): - """ - Given a ds9 region file, extract pv diagrams for each: - group of points [NOT IMPLEMENTED] - panda [NOT IMPLEMENTED] - vector [NOT IMPLEMENTED] - segment [NOT IMPLEMENTED] - group of lines [NOT IMPLEMENTED] - """ - #import pyregion - #regions = pyregion.open(regfile) - regions = load_regions_file(regfile) - return paths_from_regions(regions) - -def paths_from_regions(regions): - paths = [region_converters[r.name](r) - for r in regions - if r.name in region_converters] - return paths diff --git a/glue/external/pvextractor/utils/__init__.py b/glue/external/pvextractor/utils/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/glue/external/pvextractor/utils/wcs_slicing.py b/glue/external/pvextractor/utils/wcs_slicing.py deleted file mode 100644 index e4ad4ff73..000000000 --- a/glue/external/pvextractor/utils/wcs_slicing.py +++ /dev/null @@ -1,43 +0,0 @@ -import numpy as np - -from astropy import units as u -from astropy.wcs import WCS, WCSSUB_SPECTRAL -from .wcs_utils import get_spectral_scale - - -def slice_wcs(wcs, spatial_scale): - """ - Slice a WCS header for a spectral cube to a Position-Velocity WCS, with - ctype "OFFSET" for the spatial offset direction - - Parameters - ---------- - wcs : :class:`~astropy.wcs.WCS` - The WCS of the spectral cube. This should already be sanitized and - have the spectral axis along the third dimension. - spatial_scale: :class:`~astropy.units.Quantity` - The spatial scale of the position axis - - Returns - ------- - wcs_slice :class:`~astropy.wcs.WCS` - The resulting WCS slice - """ - - # Extract spectral slice - wcs_slice = wcs.sub([0, WCSSUB_SPECTRAL]) - - # Set spatial parameters - wcs_slice.wcs.crpix[0] = 1. - wcs_slice.wcs.cdelt[0] = spatial_scale.to(u.degree).value - wcs_slice.wcs.crval[0] = 0. - wcs_slice.wcs.ctype[0] = "OFFSET" - wcs_slice.wcs.cunit[0] = 'deg' - # Not clear why this is needed, but apparently sub with 0 sets pc[1,0] = 1, - # which is incorrect - try: - wcs_slice.wcs.pc[1,0] = wcs_slice.wcs.pc[0,1] = 0 - except AttributeError: - pass - - return wcs_slice diff --git a/glue/external/pvextractor/utils/wcs_utils.py b/glue/external/pvextractor/utils/wcs_utils.py deleted file mode 100644 index a9b368a43..000000000 --- a/glue/external/pvextractor/utils/wcs_utils.py +++ /dev/null @@ -1,71 +0,0 @@ -import numpy as np -from astropy import units as u -from astropy.wcs import WCSSUB_CELESTIAL, WCSSUB_SPECTRAL - - -def get_spatial_scale(wcs, assert_square=True): - - # Code adapted from APLpy - - wcs = wcs.sub([WCSSUB_CELESTIAL]) - cdelt = np.matrix(wcs.wcs.get_cdelt()) - pc = np.matrix(wcs.wcs.get_pc()) - scale = np.array(cdelt * pc) - - if assert_square: - try: - np.testing.assert_almost_equal(abs(cdelt[0,0]), abs(cdelt[0,1])) - np.testing.assert_almost_equal(abs(pc[0,0]), abs(pc[1,1])) - np.testing.assert_almost_equal(abs(scale[0,0]), abs(scale[0,1])) - except AssertionError: - raise ValueError("Non-square pixels. Please resample data.") - - return abs(scale[0,0]) * u.Unit(wcs.wcs.cunit[0]) - - -def get_spectral_scale(wcs): - - # Code adapted from APLpy - - wcs = wcs.sub([WCSSUB_SPECTRAL]) - cdelt = np.matrix(wcs.wcs.get_cdelt()) - pc = np.matrix(wcs.wcs.get_pc()) - scale = np.array(cdelt * pc) - - return abs(scale[0,0]) * u.Unit(wcs.wcs.cunit[0]) - - -def sanitize_wcs(mywcs): - pc = np.matrix(mywcs.wcs.get_pc()) - if (pc[:,2].sum() != pc[2,2] or pc[2,:].sum() != pc[2,2]): - raise ValueError("Non-independent 3rd axis.") - axtypes = mywcs.get_axis_types() - if ((axtypes[0]['coordinate_type'] != 'celestial' or - axtypes[1]['coordinate_type'] != 'celestial' or - axtypes[2]['coordinate_type'] != 'spectral')): - cunit3 = mywcs.wcs.cunit[2] - ctype3 = mywcs.wcs.ctype[2] - if cunit3 != '': - cunit3 = u.Unit(cunit3) - if cunit3.is_equivalent(u.m/u.s): - mywcs.wcs.ctype[2] = 'VELO' - elif cunit3.is_equivalent(u.Hz): - mywcs.wcs.ctype[2] = 'FREQ' - elif cunit3.is_equivalent(u.m): - mywcs.wcs.ctype[2] = 'WAVE' - else: - raise ValueError("Could not determine type of 3rd axis.") - elif ctype3 != '': - if 'VELO' in ctype3: - mywcs.wcs.ctype[2] = 'VELO' - elif 'FELO' in ctype3: - mywcs.wcs.ctype[2] = 'VELO-F2V' - elif 'FREQ' in ctype3: - mywcs.wcs.ctype[2] = 'FREQ' - elif 'WAVE' in ctype3: - mywcs.wcs.ctype[2] = 'WAVE' - else: - raise ValueError("Could not determine type of 3rd axis.") - else: - raise ValueError("Cube axes not in expected orientation: PPV") - return mywcs diff --git a/glue/plugins/tools/pv_slicer/qt/pv_slicer.py b/glue/plugins/tools/pv_slicer/qt/pv_slicer.py index 1d802a6ea..76df4b647 100644 --- a/glue/plugins/tools/pv_slicer/qt/pv_slicer.py +++ b/glue/plugins/tools/pv_slicer/qt/pv_slicer.py @@ -199,7 +199,7 @@ def _slice_from_path(x, y, data, attribute, slc): :note: For >3D cubes, the "V-axis" of the PV slice is the longest cube axis ignoring the x/y axes of `slc` """ - from glue.external.pvextractor import Path, extract_pv_slice + from pvextractor import Path, extract_pv_slice p = Path(list(zip(x, y))) cube = data[attribute] @@ -207,7 +207,12 @@ def _slice_from_path(x, y, data, attribute, slc): s = list(slc) ind = _slice_index(data, slc) - cube_wcs = getattr(data.coords, 'wcs', None) + from astropy.wcs import WCS + + if isinstance(data.coords, WCS): + cube_wcs = data.coords + else: + cube_wcs = None # transpose cube to (z, y, x, ) def _swap(x, s, i, j): @@ -231,8 +236,6 @@ def _swap(x, s, i, j): spacing = 1 # pixel x, y = [np.round(_x).astype(int) for _x in p.sample_points(spacing)] - from astropy.wcs import WCS - try: result = extract_pv_slice(cube, path=p, wcs=cube_wcs, order=0) wcs = WCS(result.header) diff --git a/setup.cfg b/setup.cfg index 345969894..190001cf5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -42,6 +42,7 @@ install_requires = openpyxl>=3.0 h5py>=2.10 mpl-scatter-density>=0.7 + pvextractor>=0.2 [options.entry_points] glue.plugins = @@ -99,7 +100,6 @@ test = [options.package_data] * = *.png, *.ui, *.glu, *.hdf5, *.fits, *.xlsx, *.txt, *.csv, *.svg, *.vot glue.core.data_factories.tests = data/*.jpg -glue.external.pvextractor = LICENSE glue.viewers.histogram.qt.tests = data/*.glu glue.viewers.image.qt.tests = data/*.glu, baseline/*.png glue.viewers.profile.qt.tests = data/*.glu From f65031638ac6e9f7fe32ba7e2f01007e8fe2f7ff Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 27 Jan 2022 17:42:57 +0000 Subject: [PATCH 2/5] Added changelog entry --- CHANGES.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index db15be515..9193fbb55 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,6 +1,12 @@ Full changelog ============== +v1.3.0 (unreleased) +------------------- + +* Remove bundled version of pvextractor and include it as a proper + dependency. [#2252] + v1.2.4 (2022-01-27) ------------------- From ff72ed275c9af2e2dbb0de383b07791444a8086c Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 27 Jan 2022 21:17:55 +0000 Subject: [PATCH 3/5] Remove usage of nan* functions from glue.utils --- glue/core/util.py | 8 +++----- glue/utils/array.py | 3 +-- glue/viewers/histogram/python_export.py | 6 +++--- glue/viewers/profile/qt/profile_tools.py | 12 ++++++------ glue/viewers/profile/state.py | 10 +++++----- glue/viewers/scatter/layer_artist.py | 4 ++-- 6 files changed, 20 insertions(+), 23 deletions(-) diff --git a/glue/core/util.py b/glue/core/util.py index 6a59a45ea..971669745 100644 --- a/glue/core/util.py +++ b/glue/core/util.py @@ -11,8 +11,6 @@ from matplotlib.projections.polar import ThetaFormatter, ThetaLocator import matplotlib.ticker as mticker -from glue.utils import nanmin, nanmax - __all__ = ["ThetaRadianFormatter", "relim", "split_component_view", "join_component_view", "facet_subsets", "colorize_subsets", "disambiguate", 'small_view', 'small_view_array', 'visible_limits', @@ -205,9 +203,9 @@ def facet_subsets(data_collection, cid, lo=None, hi=None, steps=5, raise ValueError("Cannot infer data limits for ComponentID %s" % cid) if lo is None: - lo = nanmin(vals) + lo = np.nanmin(vals) if hi is None: - hi = nanmax(vals) + hi = np.nanmax(vals) reverse = lo > hi if log: @@ -362,7 +360,7 @@ def visible_limits(artists, axis): if data.size == 0: return - lo, hi = nanmin(data), nanmax(data) + lo, hi = np.nanmin(data), np.nanmax(data) if not np.isfinite(lo): return diff --git a/glue/utils/array.py b/glue/utils/array.py index 49275fb62..9b7846a0f 100644 --- a/glue/utils/array.py +++ b/glue/utils/array.py @@ -10,8 +10,7 @@ __all__ = ['unique', 'shape_to_string', 'view_shape', 'stack_view', 'coerce_numeric', 'check_sorted', 'broadcast_to', 'unbroadcast', - 'iterate_chunks', 'combine_slices', 'nanmean', 'nanmedian', 'nansum', - 'nanmin', 'nanmax', 'format_minimal', 'compute_statistic', + 'iterate_chunks', 'combine_slices', 'format_minimal', 'compute_statistic', 'categorical_ndarray', 'index_lookup', 'ensure_numerical', 'broadcast_arrays_minimal', 'random_views_for_dask_array'] diff --git a/glue/viewers/histogram/python_export.py b/glue/viewers/histogram/python_export.py index 38284663c..87856676f 100644 --- a/glue/viewers/histogram/python_export.py +++ b/glue/viewers/histogram/python_export.py @@ -1,9 +1,9 @@ from distutils.version import LooseVersion +import numpy as np from matplotlib import __version__ from glue.viewers.common.python_export import code, serialize_options -from glue.utils import nanmin, nanmax MATPLOTLIB_GE_30 = LooseVersion(__version__) > '3' @@ -17,8 +17,8 @@ def python_export_histogram_layer(layer, *args): imports = ["import numpy as np"] x = layer.layer[layer._viewer_state.x_att] - x_min = nanmin(x) - x_max = nanmax(x) + x_min = np.nanmin(x) + x_max = np.nanmax(x) hist_x_min = layer._viewer_state.hist_x_min hist_x_max = layer._viewer_state.hist_x_max diff --git a/glue/viewers/profile/qt/profile_tools.py b/glue/viewers/profile/qt/profile_tools.py index f0e44375c..59dd2a7bb 100644 --- a/glue/viewers/profile/qt/profile_tools.py +++ b/glue/viewers/profile/qt/profile_tools.py @@ -8,7 +8,7 @@ from qtpy.QtCore import Qt from qtpy import QtWidgets, QtGui -from glue.utils import color2hex, nanmean, nanmedian, nanmin, nanmax, nansum +from glue.utils import color2hex from glue.config import fit_plugin, viewer_tool from glue.utils.qt import load_ui, fix_tab_widget_fontsize from glue.viewers.profile.qt.mouse_mode import NavigateMouseMode, RangeMouseMode @@ -28,11 +28,11 @@ MODES = ['navigate', 'fit', 'collapse'] -COLLAPSE_FUNCS = OrderedDict([(nanmean, 'Mean'), - (nanmedian, 'Median'), - (nanmin, 'Minimum'), - (nanmax, 'Maximum'), - (nansum, 'Sum'), +COLLAPSE_FUNCS = OrderedDict([(np.nanmean, 'Mean'), + (np.nanmedian, 'Median'), + (np.nanmin, 'Minimum'), + (np.nanmax, 'Maximum'), + (np.nansum, 'Sum'), (mom1, 'Moment 1'), (mom2, 'Moment 2')]) diff --git a/glue/viewers/profile/state.py b/glue/viewers/profile/state.py index 5a2544600..37df108c1 100644 --- a/glue/viewers/profile/state.py +++ b/glue/viewers/profile/state.py @@ -9,7 +9,7 @@ DeferredDrawCallbackProperty as DDCProperty, DeferredDrawSelectionCallbackProperty as DDSCProperty) from glue.core.data_combo_helper import ManualDataComboHelper, ComponentIDComboHelper -from glue.utils import defer_draw, nanmin, nanmax +from glue.utils import defer_draw from glue.core.link_manager import is_convertible_to_single_pixel_cid from glue.core.exceptions import IncompatibleDataException @@ -132,8 +132,8 @@ def _reset_y_limits(self, *event): if profile is not None: x, y = profile if len(y) > 0: - y_min = min(y_min, nanmin(y)) - y_max = max(y_max, nanmax(y)) + y_min = min(y_min, np.nanmin(y)) + y_max = max(y_max, np.nanmax(y)) with delay_callback(self, 'y_min', 'y_max'): if y_max > y_min: self.y_min = y_min @@ -325,5 +325,5 @@ def update_limits(self, update_profile=True): if update_profile: self.update_profile(update_limits=False) if self._profile_cache is not None and len(self._profile_cache[1]) > 0: - self.v_min = nanmin(self._profile_cache[1]) - self.v_max = nanmax(self._profile_cache[1]) + self.v_min = np.nanmin(self._profile_cache[1]) + self.v_max = np.nanmax(self._profile_cache[1]) diff --git a/glue/viewers/scatter/layer_artist.py b/glue/viewers/scatter/layer_artist.py index c2bc3da7e..9c1ffa1fc 100644 --- a/glue/viewers/scatter/layer_artist.py +++ b/glue/viewers/scatter/layer_artist.py @@ -9,7 +9,7 @@ from astropy.visualization import (ImageNormalize, LinearStretch, SqrtStretch, AsinhStretch, LogStretch) -from glue.utils import defer_draw, broadcast_to, nanmax, ensure_numerical, datetime64_to_mpl +from glue.utils import defer_draw, broadcast_to, ensure_numerical, datetime64_to_mpl from glue.viewers.scatter.state import ScatterLayerState from glue.viewers.scatter.python_export import python_export_scatter_layer from glue.viewers.matplotlib.layer_artist import MatplotlibLayerArtist @@ -57,7 +57,7 @@ def min(self, array): return 0 def max(self, array): - return 10. ** (np.log10(nanmax(array)) * self.contrast) + return 10. ** (np.log10(np.nanmax(array)) * self.contrast) def set_mpl_artist_cmap(artist, values, state=None, cmap=None, vmin=None, vmax=None): From 65389353fc07b87b2f19973721089069827d356c Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 27 Jan 2022 21:20:24 +0000 Subject: [PATCH 4/5] Fix codestyle --- glue/utils/array.py | 2 +- glue/viewers/scatter/layer_artist.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/glue/utils/array.py b/glue/utils/array.py index 9b7846a0f..9c74909d2 100644 --- a/glue/utils/array.py +++ b/glue/utils/array.py @@ -6,7 +6,7 @@ import pandas as pd # For backward-compatibility with when we used to define bottleneck wrappers -from numpy import nanmin, nanmax, nanmean, nanmedian, nansum +from numpy import nanmin, nanmax, nanmean, nanmedian, nansum # noqa __all__ = ['unique', 'shape_to_string', 'view_shape', 'stack_view', 'coerce_numeric', 'check_sorted', 'broadcast_to', 'unbroadcast', diff --git a/glue/viewers/scatter/layer_artist.py b/glue/viewers/scatter/layer_artist.py index 9c1ffa1fc..12585a3b6 100644 --- a/glue/viewers/scatter/layer_artist.py +++ b/glue/viewers/scatter/layer_artist.py @@ -295,7 +295,7 @@ def _update_data(self): hw = 1 hl = 0 - vmax = nanmax(np.hypot(vx, vy)) + vmax = np.nanmax(np.hypot(vx, vy)) self.vector_artist = self.axes.quiver(x, y, vx, vy, units='width', pivot=self.state.vector_origin, From af93ccad9c92aa7bb3e6ed588be120483896e367 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 27 Jan 2022 22:22:52 +0000 Subject: [PATCH 5/5] Fix tests --- glue/core/tests/test_state_objects.py | 3 +-- glue/viewers/profile/qt/tests/test_profile_tools.py | 5 ++--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/glue/core/tests/test_state_objects.py b/glue/core/tests/test_state_objects.py index 83c98a8c7..163d64f55 100644 --- a/glue/core/tests/test_state_objects.py +++ b/glue/core/tests/test_state_objects.py @@ -4,7 +4,6 @@ from echo import CallbackProperty, ListCallbackProperty from glue.core import Data, DataCollection -from glue.utils import nanmedian from .test_state import clone from ..state_objects import (State, StateAttributeLimitsHelper, @@ -202,7 +201,7 @@ class SimpleState(State): self.state = SimpleState() self.helper = StateAttributeSingleValueHelper(self.state, attribute='comp', - function=nanmedian, value='val') + function=np.nanmedian, value='val') self.state.data = self.data diff --git a/glue/viewers/profile/qt/tests/test_profile_tools.py b/glue/viewers/profile/qt/tests/test_profile_tools.py index 4910b6805..3e3e1a5fe 100644 --- a/glue/viewers/profile/qt/tests/test_profile_tools.py +++ b/glue/viewers/profile/qt/tests/test_profile_tools.py @@ -6,7 +6,6 @@ from glue.core import Data from glue.tests.helpers import PYSIDE2_INSTALLED # noqa from glue.app.qt import GlueApplication -from glue.utils import nanmean from glue.utils.qt import process_events from glue.viewers.image.state import AggregateSlice @@ -159,7 +158,7 @@ def test_collapse(self): assert image_viewer.state.slices[0].slice.start == 1 assert image_viewer.state.slices[0].slice.stop == 15 assert image_viewer.state.slices[0].center == 0 - assert image_viewer.state.slices[0].function is nanmean + assert image_viewer.state.slices[0].function is np.nanmean # Next, try in world coordinates @@ -181,7 +180,7 @@ def test_collapse(self): assert image_viewer.state.slices[0].slice.start == 1 assert image_viewer.state.slices[0].slice.stop == 15 assert image_viewer.state.slices[0].center == 0 - assert image_viewer.state.slices[0].function is nanmean + assert image_viewer.state.slices[0].function is np.nanmean def test_collapse_reverse(self, capsys):