Merge pull request #451 from IMBalENce/ctf_reader

Add CTF file reader and interactive crystal map plot example

CHANGELOG.rst

@@ -11,15 +11,20 @@ Unreleased
 
 Added
 -----
+- We can now read 2D crystal maps from Channel Text Files (CTFs) using ``io.load()``.
 
 Changed
 -------
+- Phase names in crystal maps read from .ang files with ``io.load()`` now prefer to use
+  the abbreviated "Formula" instead of "MaterialName" in the file header.
 
 Removed
 -------
 
 Deprecated
 ----------
+- ``loadang()`` and ``loadctf()`` are deprecated and will be removed in the next minor
+  release. Please use ``io.load()`` instead.
 
 Fixed
 -----

examples/plotting/interactive_xmap.py

@@ -0,0 +1,98 @@
+"""
+============================
+Interactive crystal map plot
+============================
+
+This example shows how to use
+:doc:`matplotlib event connections <matplotlib:users/explain/figure/event_handling>` to
+add an interactive click function to a :class:`~orix.crystal_map.CrystalMap` plot.
+Here, we navigate an inverse pole figure (IPF) map and retreive the phase name and
+corresponding Euler angles from the location of the click.
+
+.. note::
+
+    This example uses the interactive capabilities of Matplotlib, and this will not
+    appear in the static documentation.
+    Please run this code on your machine to see the interactivity.
+
+    You can copy and paste individual parts, or download the entire example using the
+    link at the bottom of the page.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from orix import data, plot
+from orix.crystal_map import CrystalMap
+
+xmap = data.sdss_ferrite_austenite(allow_download=True)
+print(xmap)
+
+pg_laue = xmap.phases[1].point_group.laue
+O_au = xmap["austenite"].orientations
+O_fe = xmap["ferrite"].orientations
+
+# Get IPF colors
+ipf_key = plot.IPFColorKeyTSL(pg_laue)
+rgb_au = ipf_key.orientation2color(O_au)
+rgb_fe = ipf_key.orientation2color(O_fe)
+
+# Combine IPF color arrays
+rgb_all = np.zeros((xmap.size, 3))
+phase_id_au = xmap.phases.id_from_name("austenite")
+phase_id_fe = xmap.phases.id_from_name("ferrite")
+rgb_all[xmap.phase_id == phase_id_au] = rgb_au
+rgb_all[xmap.phase_id == phase_id_fe] = rgb_fe
+
+
+def select_point(xmap: CrystalMap, rgb_all: np.ndarray) -> tuple[int, int]:
+    """Return location of interactive user click on image.
+
+    Interactive function for showing the phase name and Euler angles
+    from the click-position.
+    """
+    fig = xmap.plot(
+        rgb_all,
+        overlay="dp",
+        return_figure=True,
+        figure_kwargs={"figsize": (12, 8)},
+    )
+    ax = fig.axes[0]
+    ax.set_title("Click position")
+
+    # Extract array in the plot with IPF colors + dot product overlay
+    rgb_dp_2d = ax.images[0].get_array()
+
+    x = y = 0
+
+    def on_click(event):
+        x, y = (event.xdata, event.ydata)
+        if x is None:
+            print("Please click inside the IPF map")
+            return
+        print(x, y)
+
+        # Extract location in crystal map and extract phase name and
+        # Euler angles
+        xmap_yx = xmap[int(np.round(y)), int(np.round(x))]
+        phase_name = xmap_yx.phases_in_data[:].name
+        eu = xmap_yx.rotations.to_euler(degrees=True)[0].round(2)
+
+        # Format Euler angles
+        eu_str = "(" + ", ".join(np.array_str(eu)[1:-1].split()) + ")"
+
+        plt.clf()
+        plt.imshow(rgb_dp_2d)
+        plt.plot(x, y, "+", c="k", markersize=15, markeredgewidth=3)
+        plt.title(
+            f"Phase: {phase_name}, Euler angles: $(\phi_1, \Phi, \phi_2)$ = {eu_str}"
+        )
+        plt.draw()
+
+    fig.canvas.mpl_connect("button_press_event", on_click)
+
+    return x, y
+
+
+x, y = select_point(xmap, rgb_all)
+plt.show()

orix/crystal_map/crystal_map.py

@@ -17,7 +17,7 @@
 # along with orix.  If not, see <http://www.gnu.org/licenses/>.
 
 import copy
-from typing import Optional, Tuple, Union
+from typing import Dict, Optional, Tuple, Union
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -350,12 +350,12 @@ def y(self) -> Union[None, np.ndarray]:
     @property
     def dx(self) -> float:
         """Return the x coordinate step size."""
-        return self._step_size_from_coordinates(self._x)
+        return _step_size_from_coordinates(self._x)
 
     @property
     def dy(self) -> float:
         """Return the y coordinate step size."""
-        return self._step_size_from_coordinates(self._y)
+        return _step_size_from_coordinates(self._y)
 
     @property
     def row(self) -> Union[None, np.ndarray]:
@@ -1039,29 +1039,9 @@ def plot(
         if return_figure:
             return fig
 
-    @staticmethod
-    def _step_size_from_coordinates(coordinates: np.ndarray) -> float:
-        """Return step size in input ``coordinates`` array.
-
-        Parameters
-        ----------
-        coordinates
-            Linear coordinate array.
-
-        Returns
-        -------
-        step_size
-            Step size in ``coordinates`` array.
-        """
-        unique_sorted = np.sort(np.unique(coordinates))
-        step_size = 0
-        if unique_sorted.size != 1:
-            step_size = unique_sorted[1] - unique_sorted[0]
-        return step_size
-
     def _data_slices_from_coordinates(self, only_is_in_data: bool = True) -> tuple:
-        """Return a tuple of slices defining the current data extent in
-        all directions.
+        """Return a slices defining the current data extent in all
+        directions.
 
         Parameters
         ----------
@@ -1072,23 +1052,14 @@ def _data_slices_from_coordinates(self, only_is_in_data: bool = True) -> tuple:
         Returns
         -------
         slices
-            Data slice in each existing dimension, in (z, y, x) order.
+            Data slice in each existing direction in (y, x) order.
         """
         if only_is_in_data:
             coordinates = self._coordinates
         else:
             coordinates = self._all_coordinates
-
-        # Loop over dimension coordinates and step sizes
-        slices = []
-        for coords, step in zip(coordinates.values(), self._step_sizes.values()):
-            if coords is not None and step != 0:
-                c_min, c_max = np.min(coords), np.max(coords)
-                i_min = int(np.around(c_min / step))
-                i_max = int(np.around((c_max / step) + 1))
-                slices.append(slice(i_min, i_max))
-
-        return tuple(slices)
+        slices = _data_slices_from_coordinates(coordinates, self._step_sizes)
+        return slices
 
     def _data_shape_from_coordinates(self, only_is_in_data: bool = True) -> tuple:
         """Return data shape based upon coordinate arrays.
@@ -1102,21 +1073,78 @@ def _data_shape_from_coordinates(self, only_is_in_data: bool = True) -> tuple:
         Returns
         -------
         data_shape
-            Shape of data in all existing dimensions, in (z, y, x) order.
+            Shape of data in each existing direction in (y, x) order.
         """
         data_shape = []
         for dim_slice in self._data_slices_from_coordinates(only_is_in_data):
             data_shape.append(dim_slice.stop - dim_slice.start)
         return tuple(data_shape)
 
 
+def _data_slices_from_coordinates(
+    coords: Dict[str, np.ndarray], steps: Union[Dict[str, float], None] = None
+) -> Tuple[slice]:
+    """Return a list of slices defining the current data extent in all
+    directions.
+
+    Parameters
+    ----------
+    coords
+        Dictionary with coordinate arrays.
+    steps
+        Dictionary with step sizes in each direction. If not given, they
+        are computed from *coords*.
+
+    Returns
+    -------
+    slices
+        Data slice in each direction.
+    """
+    if steps is None:
+        steps = {
+            "x": _step_size_from_coordinates(coords["x"]),
+            "y": _step_size_from_coordinates(coords["y"]),
+        }
+    slices = []
+    for coords, step in zip(coords.values(), steps.values()):
+        if coords is not None and step != 0:
+            c_min, c_max = np.min(coords), np.max(coords)
+            i_min = int(np.around(c_min / step))
+            i_max = int(np.around((c_max / step) + 1))
+            slices.append(slice(i_min, i_max))
+    slices = tuple(slices)
+    return slices
+
+
+def _step_size_from_coordinates(coordinates: np.ndarray) -> float:
+    """Return step size in input *coordinates* array.
+
+    Parameters
+    ----------
+    coordinates
+        Linear coordinate array.
+
+    Returns
+    -------
+    step_size
+        Step size in *coordinates* array.
+    """
+    unique_sorted = np.sort(np.unique(coordinates))
+    if unique_sorted.size != 1:
+        step_size = unique_sorted[1] - unique_sorted[0]
+    else:
+        step_size = 0
+    return step_size
+
+
 def create_coordinate_arrays(
     shape: Optional[tuple] = None, step_sizes: Optional[tuple] = None
 ) -> Tuple[dict, int]:
-    """Create flattened coordinate arrays from a given map shape and
+    """Return flattened coordinate arrays from a given map shape and
     step sizes, suitable for initializing a
-    :class:`~orix.crystal_map.CrystalMap`. Arrays for 1D or 2D maps can
-    be returned.
+    :class:`~orix.crystal_map.CrystalMap`.
+
+    Arrays for 1D or 2D maps can be returned.
 
     Parameters
     ----------
@@ -1125,13 +1153,13 @@ def create_coordinate_arrays(
         and ten columns.
     step_sizes
         Map step sizes. If not given, it is set to 1 px in each map
-        direction given by ``shape``.
+        direction given by *shape*.
 
     Returns
     -------
     d
-        Dictionary with keys ``"y"`` and ``"x"``, depending on the
-        length of ``shape``, with coordinate arrays.
+        Dictionary with keys ``"x"`` and ``"y"``, depending on the
+        length of *shape*, with coordinate arrays.
     map_size
         Number of map points.
 
@@ -1145,10 +1173,10 @@ def create_coordinate_arrays(
     >>> create_coordinate_arrays((2, 3), (1.5, 1.5))
     ({'x': array([0. , 1.5, 3. , 0. , 1.5, 3. ]), 'y': array([0. , 0. , 0. , 1.5, 1.5, 1.5])}, 6)
     """
-    if shape is None:
+    if not shape:
         shape = (5, 10)
     ndim = len(shape)
-    if step_sizes is None:
+    if not step_sizes:
         step_sizes = (1,) * ndim
 
     if ndim == 3 or len(step_sizes) > 2:

orix/io/__init__.py

@@ -37,6 +37,7 @@
 from h5py import File, is_hdf5
 import numpy as np
 
+from orix._util import deprecated
 from orix.crystal_map import CrystalMap
 from orix.io.plugins import plugin_list
 from orix.io.plugins._h5ebsd import hdf5group2dict
@@ -45,7 +46,6 @@
 extensions = [plugin.file_extensions for plugin in plugin_list if plugin.writes]
 
 
-# Lists what will be imported when calling "from orix.io import *"
 __all__ = [
     "loadang",
     "loadctf",
@@ -54,6 +54,8 @@
 ]
 
 
+# TODO: Remove after 0.13.0
+@deprecated(since="0.13", removal="0.14", alternative="io.load")
 def loadang(file_string: str) -> Rotation:
     """Load ``.ang`` files.
 
@@ -73,6 +75,8 @@ def loadang(file_string: str) -> Rotation:
     return Rotation.from_euler(euler)
 
 
+# TODO: Remove after 0.13.0
+@deprecated(since="0.13", removal="0.14", alternative="io.load")
 def loadctf(file_string: str) -> Rotation:
     """Load ``.ctf`` files.
 

orix/io/plugins/__init__.py

@@ -28,15 +28,17 @@
 
     ang
     bruker_h5ebsd
+    ctf
     emsoft_h5ebsd
     orix_hdf5
 """
 
-from orix.io.plugins import ang, bruker_h5ebsd, emsoft_h5ebsd, orix_hdf5
+from orix.io.plugins import ang, bruker_h5ebsd, ctf, emsoft_h5ebsd, orix_hdf5
 
 plugin_list = [
     ang,
     bruker_h5ebsd,
+    ctf,
     emsoft_h5ebsd,
     orix_hdf5,
 ]