Docs about custom coordinates

CHANGES.md

@@ -4,6 +4,10 @@ Full changelog
 v0.15.0 (unreleased)
 --------------------
 
+* Added a new ``glue.core.coordinates.AffineCoordinates`` class for common
+  affine transformations, and also added documentation on defining custom
+  coordinates. [#1994]
+
 * Make it possible to view only a subset of data in the table
   viewer. [#1988]
 

doc/customizing_guide/coordinates.rst

@@ -0,0 +1,122 @@
+Customizing the coordinate system of a data object
+==================================================
+
+Background
+----------
+
+Data objects represented by the :class:`~glue.core.data.Data` class can have
+a coordinate system defined, for display and/or linking purposes. This
+coordinate system is defined in the ``.coords`` attribute of data objects::
+
+    >>> from glue.core import Data
+    >>> data = Data(x=[1, 2, 3])
+    >>> data.coords
+    <glue.core.coordinates.Coordinates object at 0x7fa52f5547b8>
+
+This attribute can be used to convert pixel to so-called 'world' coordinates
+and vice-versa::
+
+    >>> data.coords.pixel2world(2)
+    (2,)
+    >>> data.coords.world2pixel(3)
+    (3,)
+
+By default the ``coords`` object for :class:`~glue.core.data.Data` objects
+created manually is an instance of :class:`~glue.core.coordinates.Coordinates`
+which is an identity transform, as can be seen above. However, it is possible to
+use other coordinate systems or define your own.
+
+Affine coordinates
+------------------
+
+The most common cases of transformation between pixel and world coordinates are
+`affine transformations <https://en.wikipedia.org/wiki/Affine_transformation>`_,
+which can represent combinations of e.g. reflections, scaling, translations,
+rotations, and shear. A common way of representing an affine transformations is
+through an `augmented
+matrix <https://en.wikipedia.org/wiki/Affine_transformation>`_, which has shape
+N+1 x N+1, where N is the number of pixel and world dimensions.
+
+Glue provides a class for representing arbitrary affine transformations::
+
+     >>> from glue.core.coordinates import AffineCoordinates
+
+To initialize it, you will need to provide an augmented matrix, and optionally
+lists of units and axis names (as strings). For example, to construct an affine
+transformation where the x and y coordinates are each doubled, you would do::
+
+     >>> import numpy as np
+     >>> matrix = np.array([[2, 0, 0], [0, 2, 0], [0, 0, 1]])
+     >>> affine_coords = AffineCoordinates(matrix, units=['m', 'm'], labels=['xw', 'yw'])
+
+To use a custom coordinate system, when creating a data object you should specify
+the coordinates object via the ``coords=`` keyword argument::
+
+   >>> data_double = Data(x=[1, 2, 3], coords=affine_coords)
+   >>> data_double.coords.pixel2world(2, 1)
+   [array(4.), array(2.)]
+   >>> data_double.coords.world2pixel(4.0, 2.0)
+   [array(2.), array(1.)]
+
+Custom coordinates
+------------------
+
+If you want to define a fully customized coordinate transformation, you will
+need to either define a :class:`~glue.core.coordinates.Coordinates` subclass
+with the following methods::
+
+    from glue.core.coordinates import Coordinates
+
+
+    class MyCoordinates(Coordinates):
+
+        def pixel2world(self, *args):
+            # This should take N arguments (where N is the number of dimensions
+            # in your dataset) and assume these are 0-based pixel coordinates,
+            # then return N world coordinates with the same shape as the input.
+
+        def world2pixel(self, *args):
+            # This should take N arguments (where N is the number of dimensions
+            # in your dataset) and assume these are 0-based pixel coordinates,
+            # then return N world coordinates with the same shape as the input.
+
+        def world_axis_unit(self, axis):
+            # For a given axis (0-based) return the units of the world
+            # coordinate as a string. This is optional and will return '' by
+            # default if not defined.
+
+        def axis_label(self, axis):
+            # For a given axis (0-based) return the name of the world
+            # coordinate as a string. This is optional and will return
+            # 'World {axis}' by default if not defined.
+
+        def dependent_axes(self, axis):
+            # This should return a tuple of all the world dimensions that are
+            # correlated with the specified pixel axis. As an example, for a
+            # 2-d coordinate system rotated compared to the pixel coordinates,
+            # both world coordinates depend on both pixel coordinates, so this
+            # should return (0, 1). If all axes are independent, then this
+            # should return (axis,) (the default implementation)
+
+For example, let's consider a coordinate system where the world coordinates are
+simply scaled by a factor of two compared to the pixel coordinates. The minimal
+class implementing this would look like::
+
+    >>> from glue.core.coordinates import Coordinates
+
+    >>> class DoubleCoordinates(Coordinates):
+    ...
+    ...     def pixel2world(self, *args):
+    ...        return tuple([2.0 * x for x in args])
+    ...
+    ...     def world2pixel(self, *args):
+    ...        return ([0.5 * x for x in args])
+
+To use a custom coordinate system, when creating a data object you should specify
+the coordinates object via the ``coords=`` keyword argument::
+
+    >>> data_double = Data(x=[1, 2, 3], coords=DoubleCoordinates())
+    >>> data_double.coords.pixel2world(2)
+    (4.0,)
+    >>> data_double.coords.world2pixel(4.0)
+    [2.0]

doc/index.rst

@@ -95,6 +95,7 @@ Customizing/Hacking Glue
    customizing_guide/configuration.rst
    customizing_guide/customization.rst
    customizing_guide/writing_plugin.rst
+   customizing_guide/coordinates.rst
    python_guide/data_viewer_options.rst
    customizing_guide/custom_viewer.rst
    python_guide/liveupdate.rst

glue/core/coordinates.py

@@ -4,10 +4,12 @@
 
 import numpy as np
 
+from astropy.wcs import WCS
+
 from glue.utils import unbroadcast, broadcast_to, axis_correlation_matrix
 
 
-__all__ = ['Coordinates', 'WCSCoordinates', 'coordinates_from_header', 'coordinates_from_wcs']
+__all__ = ['Coordinates', 'AffineCoordinates', 'WCSCoordinates', 'coordinates_from_header', 'coordinates_from_wcs']
 
 
 class Coordinates(object):
@@ -231,7 +233,6 @@ class WCSCoordinates(Coordinates):
 
     def __init__(self, header=None, wcs=None):
         super(WCSCoordinates, self).__init__()
-        from astropy.wcs import WCS
 
         if header is None and wcs is None:
             raise ValueError('Must provide either FITS header or WCS or both')
@@ -347,6 +348,67 @@ def __setgluestate__(cls, rec, context):
         return cls(fits.Header.fromstring(rec['header']))
 
 
+class AffineCoordinates(WCSCoordinates):
+    """
+    Coordinates determined via an affine transformation represented by an
+    augmented matrix of shape N+1 x N+1 matrix, where N is the number of pixel
+    and world coordinates. The last column of the matrix should be used for
+    the translation term, and the last row should be set to 0 except for the
+    last column which should be 1.
+
+    Note that the order of the dimensions in the matrix (x, y) should be the
+    opposite of the order of the order of dimensions of Numpy arrays (y, x).
+    """
+
+    # Note that for now the easiest way to implement this is to sub-class
+    # WCS, which means that this will automatically work with WCSAxes. In
+    # future it would be good to make this independent from WCS but it will
+    # require changes to WCSAxes.
+
+    def __init__(self, matrix, units=None, labels=None):
+
+        if matrix.ndim != 2:
+            raise ValueError("Affine matrix should be two-dimensional")
+
+        if matrix.shape[0] != matrix.shape[1]:
+            raise ValueError("Affine matrix should be square")
+
+        if np.any(matrix[-1, :-1] != 0) or matrix[-1, -1] != 1:
+            raise ValueError("Last row of matrix should be zeros and a one")
+
+        if units is not None and len(units) != matrix.shape[0] - 1:
+            raise ValueError("Expected {0} units, got {1}".format(matrix.shape[0] - 1, len(units)))
+
+        if labels is not None and len(labels) != matrix.shape[0] - 1:
+            raise ValueError("Expected {0} labels, got {1}".format(matrix.shape[0] - 1, len(labels)))
+
+        self.matrix = matrix
+        self.units = units
+        self.labels = labels
+
+        wcs = WCS(naxis=self.matrix.shape[0] - 1)
+        wcs.wcs.cd = self.matrix[:-1, :-1]
+        wcs.wcs.crpix = np.ones(wcs.naxis)
+        wcs.wcs.crval = self.matrix[:-1, -1]
+        if labels is not None:
+            wcs.wcs.ctype = labels
+        if units is not None:
+            wcs.wcs.cunit = units
+
+        super(AffineCoordinates, self).__init__(wcs=wcs)
+
+    def __gluestate__(self, context):
+        return dict(matrix=context.do(self.matrix),
+                    labels=self.labels,
+                    units=self.units)
+
+    @classmethod
+    def __setgluestate__(cls, rec, context):
+        return cls(context.object(rec['matrix']),
+                   units=rec['units'],
+                   labels=rec['labels'])
+
+
 def coordinates_from_header(header):
     """
     Convert a FITS header into a glue Coordinates object.

glue/core/tests/test_coordinates.py

@@ -6,10 +6,11 @@
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
+from glue.core.tests.test_state import clone
 from glue.tests.helpers import requires_astropy
 
 from ..coordinates import (coordinates_from_header,
-                           WCSCoordinates,
+                           WCSCoordinates, AffineCoordinates,
                            Coordinates,
                            header_from_string)
 
@@ -370,3 +371,86 @@ def test_pixel2world_single_axis():
     assert_allclose(coord.pixel2world_single_axis(x, y, z, axis=0), [1.21004705, 1.42012044, 1.63021455])
     assert_allclose(coord.pixel2world_single_axis(x, y, z, axis=1), [1.24999002, 1.499947, 1.74985138])
     assert_allclose(coord.pixel2world_single_axis(x, y, z, axis=2), [1.5, 1.5, 1.5])
+
+
+def test_affine():
+
+    matrix = np.array([[2, 3, -1], [1, 2, 2], [0, 0, 1]])
+
+    coords = AffineCoordinates(matrix)
+
+    assert coords.axis_label(1) == 'World 1'
+    assert coords.axis_label(0) == 'World 0'
+
+    assert coords.world_axis_unit(1) == ''
+    assert coords.world_axis_unit(0) == ''
+
+    xp = np.array([1, 2, 3])
+    yp = np.array([2, 3, 4])
+
+    xw, yw = coords.pixel2world(xp, yp)
+
+    assert_allclose(xw, 2 * xp + 3 * yp - 1)
+    assert_allclose(yw, 1 * xp + 2 * yp + 2)
+
+    coords2 = clone(coords)
+
+    xw, yw = coords2.pixel2world(xp, yp)
+
+    assert_allclose(xw, 2 * xp + 3 * yp - 1)
+    assert_allclose(yw, 1 * xp + 2 * yp + 2)
+
+
+def test_affine_labels_units():
+
+    matrix = np.array([[2, 3, -1], [1, 2, 2], [0, 0, 1]])
+
+    coords = AffineCoordinates(matrix, units=['km', 'km'], labels=['xw', 'yw'])
+
+    assert coords.axis_label(1) == 'Xw'
+    assert coords.axis_label(0) == 'Yw'
+
+    assert coords.world_axis_unit(1) == 'km'
+    assert coords.world_axis_unit(0) == 'km'
+
+    xp = np.array([1, 2, 3])
+    yp = np.array([2, 3, 4])
+
+    xw, yw = coords.pixel2world(xp, yp)
+
+    assert_allclose(xw, 2 * xp + 3 * yp - 1)
+    assert_allclose(yw, 1 * xp + 2 * yp + 2)
+
+    coords2 = clone(coords)
+
+    xw, yw = coords2.pixel2world(xp, yp)
+
+    assert_allclose(xw, 2 * xp + 3 * yp - 1)
+    assert_allclose(yw, 1 * xp + 2 * yp + 2)
+
+
+def test_affine_invalid():
+
+    matrix = np.array([[2, 3, -1], [1, 2, 2], [0, 0, 1]])
+
+    with pytest.raises(ValueError) as exc:
+        AffineCoordinates(matrix[0])
+    assert exc.value.args[0] == 'Affine matrix should be two-dimensional'
+
+    with pytest.raises(ValueError) as exc:
+        AffineCoordinates(matrix[:-1])
+    assert exc.value.args[0] == 'Affine matrix should be square'
+
+    with pytest.raises(ValueError) as exc:
+        AffineCoordinates(matrix, labels=['a', 'b', 'c'])
+    assert exc.value.args[0] == 'Expected 2 labels, got 3'
+
+    with pytest.raises(ValueError) as exc:
+        AffineCoordinates(matrix, units=['km', 'km', 'km'])
+    assert exc.value.args[0] == 'Expected 2 units, got 3'
+
+    matrix[-1] = 1
+
+    with pytest.raises(ValueError) as exc:
+        AffineCoordinates(matrix)
+    assert exc.value.args[0] == 'Last row of matrix should be zeros and a one'