diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2afadc4..4cf7fa0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,20 @@
 ### Fixed
 
 
+## 0.93.4 (07-03-2022)
+
+### Changed
+- Update BCC slip system file and add a separate FCC file with same ordering as in DAMASK
+- Refactor boundary lines in `ebsd.map` class and add methods for warping lines to a DIC map
+- Refactor `linkEbsdMap` method and pass all arguments to transform estimate method
+- Remove IPython and jupyter as requirements
+- Move slip systems to `Phase` class and load automatically based on crystal stucture
+- Make Oxford bonary loader tolerate of unknown data fields
+
+### Fixed
+- Fix ebsd grain linker so it works again
+
+
 ## 0.93.3 (23-08-2021)
 
 ### Added
diff --git a/defdap/__init__.py b/defdap/__init__.py
index f8b4079..a3579b3 100644
--- a/defdap/__init__.py
+++ b/defdap/__init__.py
@@ -11,4 +11,9 @@
     'find_grain_report_freq': 100,
     # How to find grain in a HRDIC map, either 'floodfill' or 'warp'
     'hrdic_grain_finding_method': 'floodfill',
+    'slip_system_file': {
+        'FCC': 'cubic_fcc',
+        'BCC': 'cubic_bcc',
+        'HCP': 'hexagonal_noca',
+    },
 }
diff --git a/defdap/_version.py b/defdap/_version.py
index 327520f..c5d4fbf 100644
--- a/defdap/_version.py
+++ b/defdap/_version.py
@@ -1 +1 @@
-__version__ = '0.93.3'
+__version__ = '0.93.4'
diff --git a/defdap/base.py b/defdap/base.py
index 7e6fa68..81d83d0 100755
--- a/defdap/base.py
+++ b/defdap/base.py
@@ -105,8 +105,8 @@ def plotGrainNumbers(self, dilateBoundaries=False, ax=None, **kwargs):
             Set to true to dilate boundaries.
         ax : matplotlib.axes.Axes, optional
             axis to plot on, if not provided the current active axis is used.
-        kwargs : dict
-            Keyword arguments to pass to matplotlib.
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.plotting.MapPlot.addGrainNumbers`
 
         Returns
         -------
@@ -129,6 +129,8 @@ def locateGrainID(self, clickEvent=None, displaySelected=False, **kwargs):
             Click handler to use.
         displaySelected : bool, optional
             If true, plot slip traces for grain selected by click.
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.base.Map.plotDefault`
 
         """
         # Check that grains have been detected in the map
@@ -157,33 +159,41 @@ def clickGrainID(self, event, plot, displaySelected):
         ----------
         event :
             Click event.
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to capture clicks from.
         displaySelected : bool
             If true, plot the selected grain alone in pop-out window.
 
         """
-        if event.inaxes is plot.ax:
-            # grain id of selected grain
-            self.currGrainId = int(self.grains[int(event.ydata), int(event.xdata)] - 1)
-            print("Grain ID: {}".format(self.currGrainId))
-
-            # update the grain highlights layer in the plot
-            plot.addGrainHighlights([self.currGrainId], alpha=self.highlightAlpha)
-
-            if displaySelected:
-                currGrain = self[self.currGrainId]
-                if self.grainPlot is None or not self.grainPlot.exists:
-                    self.grainPlot = currGrain.plotDefault(makeInteractive=True)
-                else:
-                    self.grainPlot.clear()
-                    self.grainPlot.callingGrain = currGrain
-                    currGrain.plotDefault(plot=self.grainPlot)
-                    self.grainPlot.draw()
+        # check if click was on the map
+        if event.inaxes is not plot.ax:
+            return
+
+        # grain id of selected grain
+        self.currGrainId = int(self.grains[int(event.ydata), int(event.xdata)] - 1)
+        print("Grain ID: {}".format(self.currGrainId))
+
+        # update the grain highlights layer in the plot
+        plot.addGrainHighlights([self.currGrainId], alpha=self.highlightAlpha)
+
+        if displaySelected:
+            currGrain = self[self.currGrainId]
+            if self.grainPlot is None or not self.grainPlot.exists:
+                self.grainPlot = currGrain.plotDefault(makeInteractive=True)
+            else:
+                self.grainPlot.clear()
+                self.grainPlot.callingGrain = currGrain
+                currGrain.plotDefault(plot=self.grainPlot)
+                self.grainPlot.draw()
 
     def drawLineProfile(self, **kwargs):
         """Interactive plot for drawing a line profile of data.
 
+        Parameters
+        ----------
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.base.Map.plotDefault`
+
         """
         plot = self.plotDefault(makeInteractive=True, **kwargs)
 
@@ -198,21 +208,25 @@ def calcLineProfile(self, plot, startEnd, **kwargs):
 
         Parameters
         ----------
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to calculate the line profile for.
         startEnd : array_like
             Selected points (x0, y0, x1, y1).
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`matplotlib.pyplot.plot`
 
         """
-
         x0, y0 = startEnd[0:2]
         x1, y1 = startEnd[2:4]
         profile_length = np.sqrt((y1 - y0) ** 2 + (x1 - x0) ** 2)
 
         # Extract the values along the line
-        zi = profile_line(plot.imgLayers[0].get_array(),
-                          (startEnd[1], startEnd[0]), (startEnd[3], startEnd[2]),
-                          mode='nearest', **kwargs)
+        zi = profile_line(
+            plot.imgLayers[0].get_array(),
+            (startEnd[1], startEnd[0]),
+            (startEnd[3], startEnd[2]),
+            mode='nearest'
+        )
         xi = np.linspace(0, profile_length, len(zi))
 
         if self.profilePlot is None or not self.profilePlot.exists:
@@ -237,7 +251,7 @@ def setHomogPoint(self, binSize=1, points=None, **kwargs):
         points : numpy.ndarray, optional
             Array of (x,y) homologous points to set explicitly.
         kwargs : dict, optional
-            Keyword arguments for matplotlib.
+            Keyword arguments passed to :func:`defdap.base.Map.plotHomog`
 
         """
         if points is None:
@@ -269,17 +283,20 @@ def clickHomog(self, event, plot):
         ----------
         event :
             Click event.
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to monitor.
 
         """
-        if event.inaxes is plot.ax:
-            # right mouse click or shift + left mouse click
-            # shift click doesn't work in osx backend
-            if (event.button == 3 or
-                    (event.button == 1 and event.key == 'shift')):
-                plot.addPoints([int(event.xdata)], [int(event.ydata)],
-                               updateLayer=1)
+        # check if click was on the map
+        if event.inaxes is not plot.ax:
+            return
+
+        # right mouse click or shift + left mouse click
+        # shift click doesn't work in osx backend
+        if (event.button == 3 or
+                (event.button == 1 and event.key == 'shift')):
+            plot.addPoints([int(event.xdata)], [int(event.ydata)],
+                           updateLayer=1)
 
     def keyHomog(self, event, plot):
         """Event handler for moving position using keyboard after clicking on a map.
@@ -288,7 +305,7 @@ def keyHomog(self, event, plot):
         ----------
         event :
             Keypress event.
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to monitor.
 
         """
@@ -324,7 +341,7 @@ def clickSaveHomog(self, event, plot, binSize):
         ----------
         event :
             Button click event.
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to monitor.
         binSize : int, optional
             Binning applied to image, if applicable.
@@ -457,41 +474,43 @@ def clickGrainNeighbours(self, event, plot):
         ----------
         event :
             Click event.
-        plot : defdap.plotting.Plot
+        plot : defdap.plotting.MapPlot
             Plot to monitor.
 
         """
-        if event.inaxes is plot.ax:
-            # grain id of selected grain
-            grainId = int(self.grains[int(event.ydata), int(event.xdata)] - 1)
-            if grainId < 0:
-                return
-            self.currGrainId = grainId
-            grain = self[grainId]
+        # check if click was on the map
+        if event.inaxes is not plot.ax:
+            return
 
-            # find first and second nearest neighbours
-            firstNeighbours = list(self.neighbourNetwork.neighbors(grain))
-            highlightGrains = [grain] + firstNeighbours
-
-            secondNeighbours = []
-            for firstNeighbour in firstNeighbours:
-                trialSecondNeighbours = list(
-                    self.neighbourNetwork.neighbors(firstNeighbour)
-                )
-                for secondNeighbour in trialSecondNeighbours:
-                    if (secondNeighbour not in highlightGrains and
-                            secondNeighbour not in secondNeighbours):
-                        secondNeighbours.append(secondNeighbour)
-            highlightGrains.extend(secondNeighbours)
-
-            highlightGrains = [grain.grainID for grain in highlightGrains]
-            highlightColours = ['white']
-            highlightColours.extend(['yellow'] * len(firstNeighbours))
-            highlightColours.append('green')
-
-            # update the grain highlights layer in the plot
-            plot.addGrainHighlights(highlightGrains,
-                                    grainColours=highlightColours)
+        # grain id of selected grain
+        grainId = int(self.grains[int(event.ydata), int(event.xdata)] - 1)
+        if grainId < 0:
+            return
+        self.currGrainId = grainId
+        grain = self[grainId]
+
+        # find first and second nearest neighbours
+        firstNeighbours = list(self.neighbourNetwork.neighbors(grain))
+        highlightGrains = [grain] + firstNeighbours
+
+        secondNeighbours = []
+        for firstNeighbour in firstNeighbours:
+            trialSecondNeighbours = list(
+                self.neighbourNetwork.neighbors(firstNeighbour)
+            )
+            for secondNeighbour in trialSecondNeighbours:
+                if (secondNeighbour not in highlightGrains and
+                        secondNeighbour not in secondNeighbours):
+                    secondNeighbours.append(secondNeighbour)
+        highlightGrains.extend(secondNeighbours)
+
+        highlightGrains = [grain.grainID for grain in highlightGrains]
+        highlightColours = ['white']
+        highlightColours.extend(['yellow'] * len(firstNeighbours))
+        highlightColours.append('green')
+
+        # update the grain highlights layer in the plot
+        plot.addGrainHighlights(highlightGrains, grainColours=highlightColours)
 
     @property
     def proxigram(self):
@@ -673,8 +692,8 @@ def plotGrainDataMap(
             IDs of grains to plot for. Use -1 for all grains in the map.
         bg: int or real, optional
             Value to fill the background with.
-        kwargs:
-            Other parameters are passed to :func:`defdap.plotting.MapPlot.create`
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.plotting.MapPlot.create`
 
         Returns
         -------
@@ -710,18 +729,18 @@ def plotGrainDataIPF(
 
         Parameters
         ----------
-        mapData : numpy.ndarray
-            Array of map data to grain average. This must be cropped!
         direction : numpy.ndarray
             Vector of reference direction for the IPF.
-        plotColourBar : bool, optional
-            Set to False to exclude the colour bar from the plot.
-        vmin : float, optional
-            Minimum value of colour scale.
-        vmax : float, optional
-            Maximum value for colour scale.
-        cLabel : str, optional
-            Colour bar label text.
+        mapData : numpy.ndarray
+            Array of map data. This must be cropped! Either mapData or
+            grainData must be supplied.
+        grainData : list or np.array, optional
+            Grain values. This an be a single value per grain or RGB
+            values. You must supply either mapData or grainData.
+        grainIds: list of int or int, optional
+            IDs of grains to plot for. Use -1 for all grains in the map.
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.quat.Quat.plotIPF`
 
         """
         # Set default plot parameters then update with any input
@@ -759,6 +778,16 @@ def plotGrainDataIPF(
 class Grain(object):
     """
     Base class for a grain.
+
+    Attributes
+    ----------
+    grainID : int
+
+    ownerMap : defdap.base.Map
+
+    coordList : list of tuples
+
+
     """
     def __init__(self, grainID, ownerMap):
         # list of coords stored as tuples (x, y). These are coords in a
@@ -861,7 +890,7 @@ def plotOutline(self, ax=None, plotScaleBar=False, **kwargs):
         plotScaleBar : bool
             plots the scale bar on the grain if true.
         kwargs : dict
-            keyword arguments to pass to :func:`defdap.plotting.GrainPlot.addMap`.
+            keyword arguments passed to :func:`defdap.plotting.GrainPlot.addMap`
 
         Returns
         -------
@@ -1014,14 +1043,8 @@ def plotGrainData(self, mapData=None, grainData=None, **kwargs):
         grainData : numpy.ndarray
             List of data at each point in the grain. Either this or
             'mapData' must be supplied and 'grainData' takes precedence.
-        vmin : float, optional
-            Minimum value of colour scale.
-        vmax : float, optional
-            Minimum value of colour scale.
-        cLabel : str, optional
-            Colour bar label text.
-        cmap : str, optional
-            Colour map to use, default is viridis.
+        kwargs : dict, optional
+            Keyword arguments passed to :func:`defdap.plotting.GrainPlot.create`
 
         """
         # Set default plot parameters then update with any input
diff --git a/defdap/crystal.py b/defdap/crystal.py
index a7614a8..467a7b4 100755
--- a/defdap/crystal.py
+++ b/defdap/crystal.py
@@ -22,24 +22,25 @@
 
 
 class Phase(object):
-    def __init__(self, name, laueGroup, latticeParams,
-                 slipSystems=None):
+    def __init__(self, name, laueGroup, spaceGroup, latticeParams):
         """
         Parameters
         ----------
         name : str
             Name of the phase
+        laueGroup : int
+            Laue group
+        spaceGroup : int
+            Space group
         latticeParams : tuple
             Lattice parameters in order (a,b,c,alpha,beta,gamma)
-        crystalStructure : defdap.crystal.CrystalStructure
-            Crystal structure of this phase
-        slipSystems : collection of defdap.crystal.SlipSystem
-            Slip systems available in the phase
+
         """
         self.name = name
         self.laueGroup = laueGroup
+        self.spaceGroup = spaceGroup
         self.latticeParams = latticeParams
-        self.slipSystems = slipSystems
+
         try:
             self.crystalStructure = {
                 9: crystalStructures['hexagonal'],
@@ -48,12 +49,32 @@ def __init__(self, name, laueGroup, latticeParams,
         except KeyError:
             raise ValueError(f"Unknown Laue group key: {laueGroup}")
 
+        if self.crystalStructure is crystalStructures['hexagonal']:
+            self.ss_file = defaults['slip_system_file']['HCP']
+        else:
+            try:
+                self.ss_file = defaults['slip_system_file'][
+                    {225: 'FCC', 229: 'BCC'}[spaceGroup]
+                ]
+            except KeyError:
+                self.ss_file = None
+
+        if self.ss_file is None:
+            self.slipSystems = None
+            self.slipTraceColours = None
+        else:
+            self.slipSystems, self.slipTraceColours = SlipSystem.load(
+                self.ss_file, self.crystalStructure, cOverA=self.cOverA
+            )
+
     def __str__(self):
-        text = "Phase: {:}\n  Crystal structure: {:}\n  Lattice params: " \
-               "({:.2f}, {:.2f}, {:.2f}, {:.0f}, {:.0f}, {:.0f})"
+        text = ("Phase: {:}\n  Crystal structure: {:}\n  Lattice params: "
+                "({:.2f}, {:.2f}, {:.2f}, {:.0f}, {:.0f}, {:.0f})\n"
+                "  Slip systems: {:}")
         return text.format(self.name, self.crystalStructure.name,
                            *self.latticeParams[:3],
-                           *np.array(self.latticeParams[3:])*180/np.pi)
+                           *np.array(self.latticeParams[3:])*180/np.pi,
+                           self.ss_file)
 
     @property
     def cOverA(self):
@@ -61,6 +82,19 @@ def cOverA(self):
             return self.latticeParams[2] / self.latticeParams[0]
         return None
 
+    def printSlipSystems(self):
+        """Print a list of slip planes (with colours) and slip directions.
+
+        """
+        # TODO: this should be moved to static method of the SlipSystem class
+        for i, (ssGroup, colour) in enumerate(zip(self.slipSystems,
+                                                  self.slipTraceColours)):
+            print('Plane {0}: {1}\tColour: {2}'.format(
+                i, ssGroup[0].slipPlaneLabel, colour
+            ))
+            for j, ss in enumerate(ssGroup):
+                print('  Direction {0}: {1}'.format(j, ss.slipDirLabel))
+
 
 class CrystalStructure(object):
     def __init__(self, name, symmetries, vertices, faces):
@@ -256,7 +290,7 @@ class SlipSystem(object):
     """Class used for defining and performing operations on a slip system.
 
     """
-    def __init__(self, slipPlane, slipDir, crystalSym, cOverA=None):
+    def __init__(self, slipPlane, slipDir, crystalStructure, cOverA=None):
         """Initialise a slip system object.
 
         Parameters
@@ -265,24 +299,24 @@ def __init__(self, slipPlane, slipDir, crystalSym, cOverA=None):
             Slip plane.
         slipDir: numpy.ndarray
             Slip direction.
-        crystalSym : str
-            The crystal symmetry ("cubic" or "hexagonal").
+        crystalStructure : defdap.crystal.CrystalStructure
+            Crystal structure of the slip system.
         cOverA : float, optional
             C over a ratio for hexagonal crystals.
 
         """
-        self.crystalSym = crystalSym    # symmetry of material
+        self.crystalStructure = crystalStructure
 
         # Stored as Miller indices (Miller-Bravais for hexagonal)
         self.planeIdc = tuple(slipPlane)
         self.dirIdc = tuple(slipDir)
 
         # Stored as vectors in a cartesian basis
-        if crystalSym == "cubic":
+        if self.crystalStructure.name == "cubic":
             self.slipPlane = slipPlane / norm(slipPlane)
             self.slipDir = slipDir / norm(slipDir)
             self.cOverA = None
-        elif crystalSym == "hexagonal":
+        elif self.crystalStructure.name == "hexagonal":
             if cOverA is None:
                 raise Exception("No c over a ratio given")
             self.cOverA = cOverA
@@ -318,8 +352,9 @@ def __str__(self):
         return self.slipPlaneLabel + self.slipDirLabel
 
     def __repr__(self):
-        return f"SlipSystem(slipPlane={self.slipPlaneLabel}, " \
-               f"slipDir={self.slipDirLabel}, crystalSym={self.crystalSym})"
+        return (f"SlipSystem(slipPlane={self.slipPlaneLabel}, "
+                f"slipDir={self.slipDirLabel}, "
+                f"symmetry={self.crystalStructure.name})")
 
     @property
     def slipPlaneLabel(self):
@@ -355,14 +390,14 @@ def generateFamily(self):
 
         """
         #
-        symms = Quat.symEqv(self.crystalSym)
+        symms = self.crystalStructure.symmetries
 
         ss_family = set()  # will not preserve order
 
         plane = self.planeIdc
         dir = self.dirIdc
 
-        if self.crystalSym == 'hexagonal':
+        if self.crystalStructure.name == 'hexagonal':
             # Transformation from crystal to orthonormal coords
             lMatrix = CrystalStructure.lMatrix(
                 1, 1, self.cOverA, np.pi / 2, np.pi / 2, np.pi * 2 / 3
@@ -380,7 +415,7 @@ def generateFamily(self):
             plane_symm = symm.transformVector(plane)
             dir_symm = symm.transformVector(dir)
 
-            if self.crystalSym == 'hexagonal':
+            if self.crystalStructure.name == 'hexagonal':
                 # qMatrix inverse is equal to lMatrix transposed and vice-versa
                 plane_symm = reduceIdc(convertIdc(
                     'm', plane=safeIntCast(np.matmul(lMatrix.T, plane_symm))
@@ -392,13 +427,13 @@ def generateFamily(self):
             ss_family.add(SlipSystem(
                 posIdc(safeIntCast(plane_symm)),
                 posIdc(safeIntCast(dir_symm)),
-                self.crystalSym, cOverA=self.cOverA
+                self.crystalStructure, cOverA=self.cOverA
             ))
 
         return ss_family
 
     @staticmethod
-    def loadSlipSystems(name, crystalSym, cOverA=None, groupBy='plane'):
+    def load(name, crystalStructure, cOverA=None, groupBy='plane'):
         """
         Load in slip systems from file. 3 integers for slip plane
         normal and 3 for slip direction. Returns a list of list of slip
@@ -409,8 +444,8 @@ def loadSlipSystems(name, crystalSym, cOverA=None, groupBy='plane'):
         name : str
             Name of the slip system file (without file extension)
             stored in the defdap install dir or path to a file.
-        crystalSym : str
-            The crystal symmetry ("cubic" or "hexagonal").
+        crystalStructure : defdap.crystal.CrystalStructure
+            Crystal structure of the slip systems.
         cOverA : float, optional
             C over a ratio for hexagonal crystals.
         groupBy : str, optional
@@ -450,7 +485,7 @@ def loadSlipSystems(name, crystalSym, cOverA=None, groupBy='plane'):
         slipTraceColours = slipSystemFile.readline().strip().split(',')
         slipSystemFile.close()
 
-        if crystalSym == "hexagonal":
+        if crystalStructure.name == "hexagonal":
             vectSize = 4
         else:
             vectSize = 3
@@ -465,17 +500,17 @@ def loadSlipSystems(name, crystalSym, cOverA=None, groupBy='plane'):
         for row in ssData:
             slipSystems.append(SlipSystem(
                 row[0:vectSize], row[vectSize:2 * vectSize],
-                crystalSym, cOverA=cOverA
+                crystalStructure, cOverA=cOverA
             ))
 
         # Group slip systems is required
         if groupBy is not None:
-            slipSystems = SlipSystem.groupSlipSystems(slipSystems, groupBy)
+            slipSystems = SlipSystem.group(slipSystems, groupBy)
 
         return slipSystems, slipTraceColours
 
     @staticmethod
-    def groupSlipSystems(slipSystems, groupBy):
+    def group(slipSystems, groupBy):
         """
         Groups slip systems by their slip plane.
 
diff --git a/defdap/ebsd.py b/defdap/ebsd.py
index 87aee6e..831c030 100755
--- a/defdap/ebsd.py
+++ b/defdap/ebsd.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import numpy as np
-from matplotlib.widgets import Button
 from skimage import morphology as mph
 import networkx as nx
 
@@ -51,8 +50,6 @@ class Map(base.Map):
         Band contrast for each point of map. Shape (yDim, xDim).
     quatArray : numpy.ndarray of defdap.quat.Quat
         Quaterions for each point of map. Shape (yDim, xDim).
-    numPhases : int
-        Number of phases.
     phaseArray : numpy.ndarray
         Map of phase ids. 1-based, 0 is non-indexed points
     phases : list of defdap.crystal.Phase
@@ -71,12 +68,8 @@ class Map(base.Map):
         Map of misorientation axis components.
     kam : numpy.ndarray
         Map of KAM.
-    slipSystems : list of list of defdap.crystal.SlipSystem
-        Slip systems grouped by slip plane.
-    slipTraceColours list(str)
-        Colours used when plotting slip traces.
     origin : tuple(int)
-        Map origin (y, x). Used by linker class where origin is a
+        Map origin (x, y). Used by linker class where origin is a
         homologue point of the maps.
     GND : numpy.ndarray
         GND scalar map.
@@ -84,7 +77,6 @@ class Map(base.Map):
         3x3 Nye tensor at each point.
 
     """
-
     def __init__(self, fileName, dataType=None):
         """
         Initialise class and load EBSD data.
@@ -123,8 +115,6 @@ def __init__(self, fileName, dataType=None):
         self.origin = (0, 0)
         self.GND = None
         self.Nye = None
-        self.slipSystems = None
-        self.slipTraceColours = None
 
         # Phase used for the maps crystal structure and cOverA. So old
         # functions still work for the 'main' phase in the map. 0-based
@@ -328,7 +318,8 @@ def plotEulerMap(self, phases=None, **kwargs):
 
         return MapPlot.create(self, map_colours, **plot_params)
 
-    def plotIPFMap(self, direction, backgroundColour = [0., 0., 0.], phases=None, **kwargs):
+    def plotIPFMap(self, direction, backgroundColour=None, phases=None,
+                   **kwargs):
         """
         Plot a map with points coloured in IPF colouring,
         with respect to a given sample direction.
@@ -360,6 +351,9 @@ def plotIPFMap(self, direction, backgroundColour = [0., 0., 0.], phases=None, **
             phase_ids = phases
             phases = [self.phases[i] for i in phase_ids]
 
+        if backgroundColour is None:
+            backgroundColour = [0., 0., 0.]
+
         map_colours = np.tile(np.array(backgroundColour), self.shape + (1,))
 
         for phase, phase_id in zip(phases, phase_ids):
@@ -628,21 +622,18 @@ def checkDataLoaded(self):
         return True
 
     @reportProgress("building quaternion array")
-    def buildQuatArray(self, force = False):
+    def buildQuatArray(self, force=False):
         """Build quaternion array
 
         Parameters
         ----------
-        force, optional
+        force : bool, optional
             If true, re-build quaternion array
         """
         self.checkDataLoaded()
 
-        if force == False:
-            if self.quatArray is None:
-                # create the array of quat objects
-                self.quatArray = Quat.createManyQuats(self.eulerAngleArray)
-        if force == True:
+        if force or self.quatArray is None:
+            # create the array of quat objects
             self.quatArray = Quat.createManyQuats(self.eulerAngleArray)
 
         yield 1.
@@ -817,48 +808,18 @@ def findBoundaries(self, boundDef=10):
         self.boundaries = np.logical_or(self.boundariesX, self.boundariesY)
         self.boundaries = -self.boundaries.astype(int)
 
-        _, _, self.boundaryLines = Map.create_boundary_lines(
-            boundaries_x=self.boundariesX,
-            boundaries_y=self.boundariesY
+        _, _, self.boundaryLines = BoundarySegment.boundary_points_to_lines(
+            boundary_points_x=zip(*self.boundariesX.transpose().nonzero()),
+            boundary_points_y=zip(*self.boundariesY.transpose().nonzero())
         )
-        _, _, self.phaseBoundaryLines = Map.create_boundary_lines(
-            boundaries_x=self.phaseBoundariesX,
-            boundaries_y=self.phaseBoundariesY
+
+        _, _, self.phaseBoundaryLines = BoundarySegment.boundary_points_to_lines(
+            boundary_points_x=zip(*self.phaseBoundariesX.transpose().nonzero()),
+            boundary_points_y=zip(*self.phaseBoundariesY.transpose().nonzero())
         )
 
         yield 1.
 
-    @staticmethod
-    def create_boundary_lines(*, boundaries_x=None, boundaries_y=None):
-        boundary_data = {}
-        if boundaries_x is not None:
-            boundary_data['x'] = boundaries_x
-        if boundaries_y is not None:
-            boundary_data['y'] = boundaries_y
-        if not boundary_data:
-            raise ValueError("No boundaries provided.")
-
-        deltas = {
-            'x': (0.5, -0.5, 0.5, 0.5),
-            'y': (-0.5, 0.5, 0.5, 0.5)
-        }
-        all_lines = []
-        for mode, boundaries in boundary_data.items():
-            points = np.where(boundaries)
-            lines = []
-            for i, j in zip(*points):
-                lines.append((
-                    (j + deltas[mode][0], i + deltas[mode][1]),
-                    (j + deltas[mode][2], i + deltas[mode][3])
-                ))
-            all_lines.append(lines)
-
-        if len(all_lines) == 2:
-            all_lines.append(all_lines[0] + all_lines[1])
-            return tuple(all_lines)
-        else:
-            return all_lines[0]
-
     @reportProgress("constructing neighbour network")
     def buildNeighbourNetwork(self):
         # create network
@@ -907,7 +868,6 @@ def buildNeighbourNetwork(self):
         self.neighbourNetwork = nn
 
     @reportProgress("finding phase boundaries")
-
     def plotPhaseBoundaryMap(self, dilate=False, **kwargs):
         """Plot phase boundary map.
 
@@ -973,7 +933,6 @@ def findGrains(self, minGrainSize=10):
             Minimum grain area in pixels.
 
         """
-        # TODO: grains need to be assigned a phase
         # Initialise the grain map
         # TODO: Look at grain map compared to boundary map
         # self.grains = np.copy(self.boundaries)
@@ -1222,38 +1181,6 @@ def plotMisOriMap(self, component=0, **kwargs):
 
         return plot
 
-    def loadSlipSystems(self, name):
-        """Load slip system definitions from file.
-
-        Parameters
-        ----------
-        name : str
-            name of the slip system file (without file extension)
-            stored in the defdap install dir or path to a file.
-
-        """
-        # TODO: should be loaded into the phases of the map
-        self.slipSystems, self.slipTraceColours = SlipSystem.loadSlipSystems(
-            name, self.crystalSym, cOverA=self.cOverA
-        )
-
-        if self.checkGrainsDetected(raiseExc=False):
-            for grain in self:
-                grain.slipSystems = self.slipSystems
-
-    def printSlipSystems(self):
-        """Print a list of slip planes (with colours) and slip directions.
-
-        """
-        # TODO: this should be moved to static method of the SlipSystem class
-        for i, (ssGroup, colour) in enumerate(zip(self.slipSystems,
-                                                  self.slipTraceColours)):
-            print('Plane {0}: {1}\tColour: {2}'.format(
-                i, ssGroup[0].slipPlaneLabel, colour
-            ))
-            for j, ss in enumerate(ssGroup):
-                print('  Direction {0}: {1}'.format(j, ss.slipDirLabel))
-
     @reportProgress("calculating grain average Schmid factors")
     def calcAverageGrainSchmidFactors(self, loadVector, slipSystems=None):
         """
@@ -1265,8 +1192,8 @@ def calcAverageGrainSchmidFactors(self, loadVector, slipSystems=None):
         loadVector :
             Loading vector, e.g. [1, 0, 0].
         slipSystems : list, optional
-            Slip planes to calculate Schmid factor for,
-            maximum of all planes calculated if not given.
+            Slip planes to calculate Schmid factor for, maximum of all
+            planes calculated if not given.
 
         """
         # Check that grains have been detected in the map
@@ -1350,14 +1277,14 @@ class Grain(base.Grain):
 
     Attributes
     ----------
-    crystalSym : str
-        Symmetry of material e.g. "cubic", "hexagonal"
-    slipSystems : list(list(defdap.crystal.SlipSystem))
-        Slip systems
     ebsdMap : defdap.ebsd.Map
         EBSD map this grain is a member of.
     ownerMap : defdap.ebsd.Map
         EBSD map this grain is a member of.
+    phaseID : int
+
+    phase : defdap.crystal.Phase
+
     quatList : list
         List of quats.
     misOriList : list
@@ -1376,15 +1303,10 @@ class Grain(base.Grain):
          Angle between slip plane and screen plane.
 
     """
-
-    # TODO: each grain should be assigned a phase and slip systems
-    # slip systems accessed from the phase
     def __init__(self, grainID, ebsdMap):
         # Call base class constructor
         super(Grain, self).__init__(grainID, ebsdMap)
 
-        self.crystalSym = ebsdMap.crystalSym    # symmetry of material e.g. "cubic", "hexagonal"
-        self.slipSystems = ebsdMap.slipSystems
         self.ebsdMap = self.ownerMap            # ebsd map this grain is a member of
         self.quatList = []                      # list of quats
         self.misOriList = None                  # list of misOri at each point in grain
@@ -1402,6 +1324,11 @@ def plotDefault(self):
             *args, **kwargs
         )
 
+    @property
+    def crystalSym(self):
+        """Temporary"""
+        return self.phase.crystalStructure.name
+
     def addPoint(self, coord, quat):
         """Append a coordinate and a quat to a grain.
 
@@ -1592,7 +1519,7 @@ def calcAverageSchmidFactors(self, loadVector, slipSystems=None):
 
         """
         if slipSystems is None:
-            slipSystems = self.slipSystems
+            slipSystems = self.phase.slipSystems
         if self.refOri is None:
             self.calcAverageOri()
 
@@ -1636,14 +1563,17 @@ def printSlipTraces(self):
         """Print a list of slip planes (with colours) and slip directions
 
         """
-
         self.calcSlipTraces()
 
         if self.averageSchmidFactors is None:
             raise Exception("Run 'calcAverageGrainSchmidFactors' on the EBSD map first")
 
-        for ssGroup, colour, sfGroup, slipTrace in zip(self.slipSystems, self.ebsdMap.slipTraceColours,
-                                                       self.averageSchmidFactors, self.slipTraces):
+        for ssGroup, colour, sfGroup, slipTrace in zip(
+            self.phase.slipSystems,
+            self.phase.slipTraceColours,
+            self.averageSchmidFactors,
+            self.slipTraces
+        ):
             print('{0}\tColour: {1}\tAngle: {2:.2f}'.format(ssGroup[0].slipPlaneLabel, colour, slipTrace * 180 / np.pi))
             for ss, sf in zip(ssGroup, sfGroup):
                 print('  {0}   SF: {1:.3f}'.format(ss.slipDirLabel, sf))
@@ -1657,7 +1587,7 @@ def calcSlipTraces(self, slipSystems=None):
 
         """
         if slipSystems is None:
-            slipSystems = self.slipSystems
+            slipSystems = self.phase.slipSystems
         if self.refOri is None:
             self.calcAverageOri()
 
@@ -1774,6 +1704,15 @@ def boundaryPointPairsY(self):
         """
         return self.boundaryPointPairs(1)
 
+    @property
+    def boundaryLines(self):
+        """Return line points along this boundary segment"""
+        _, _, lines = self.boundary_points_to_lines(
+            boundary_points_x=self.boundaryPointsX,
+            boundary_points_y=self.boundaryPointsY
+        )
+        return lines
+
     def misorientation(self):
         misOri, minSymm = self.grain1.refOri.misOri(
             self.grain2.refOri, self.ebsdMap.crystalSym, returnQuat=2
@@ -1793,126 +1732,182 @@ def misorientation(self):
         #                                                      compVector))))
         # print(deviation * 180 / np.pi)
 
+    @staticmethod
+    def boundary_points_to_lines(*, boundary_points_x=None,
+                                 boundary_points_y=None):
+        boundary_data = {}
+        if boundary_points_x is not None:
+            boundary_data['x'] = boundary_points_x
+        if boundary_points_y is not None:
+            boundary_data['y'] = boundary_points_y
+        if not boundary_data:
+            raise ValueError("No boundaries provided.")
+
+        deltas = {
+            'x': (0.5, -0.5, 0.5, 0.5),
+            'y': (-0.5, 0.5, 0.5, 0.5)
+        }
+        all_lines = []
+        for mode, points in boundary_data.items():
+            lines = []
+            for i, j in points:
+                lines.append((
+                    (i + deltas[mode][0], j + deltas[mode][1]),
+                    (i + deltas[mode][2], j + deltas[mode][3])
+                ))
+            all_lines.append(lines)
+
+        if len(all_lines) == 2:
+            all_lines.append(all_lines[0] + all_lines[1])
+            return tuple(all_lines)
+        else:
+            return all_lines[0]
+
 
 class Linker(object):
     """Class for linking multiple EBSD maps of the same region for analysis of deformation.
 
-    Parameters
+    Attributes
     ----------
-    ebsdMaps : list(ebsd.Map)
-        List of ebsd.Map objects that are linked.
-    links : list
+    ebsd_maps : list(ebsd.Map)
+        List of `ebsd.Map` objects that are linked.
+    links : list(tuple(int))
         List of grain link. Each link is stored as a tuple of
         grain IDs (one from each map stored in same order of maps).
-    numMaps : int
-        Number of linked maps.
+    plots : list(plotting.MapPlot)
+        List of last opened plot of each map.
 
     """
+    def __init__(self, ebsd_maps):
+        """Initialise linker and set ebsd maps
 
-    def __init__(self, maps):
-        self.ebsdMaps = maps
-        self.numMaps = len(maps)
+        Parameters
+        ----------
+        ebsd_maps : list(ebsd.Map)
+            List of `ebsd.Map` objects that are linked.
+
+        """
+        self.ebsd_maps = ebsd_maps
         self.links = []
-        return
+        self.plots = None
 
-    def setOrigin(self):
+    def set_origin(self, **kwargs):
         """Interacive tool to set origin of each EBSD map.
 
-        """
-        for ebsdMap in self.ebsdMaps:
-            ebsdMap.locateGrainID(clickEvent=self.clickSetOrigin)
+        Parameters
+        ----------
+        kwargs
+            Keyword arguments passed to :func:`defdap.ebsd.Map.plotDefault`
 
-    def clickSetOrigin(self, event, currentEbsdMap):
+        """
+        self.plots = []
+        for ebsd_map in self.ebsd_maps:
+            plot = ebsd_map.plotDefault(makeInteractive=True, **kwargs)
+            plot.addEventHandler('button_press_event', self.click_set_origin)
+            plot.addPoints([ebsd_map.origin[0]], [ebsd_map.origin[1]],
+                           c='w', s=60, marker='x')
+            self.plots.append(plot)
+
+    def click_set_origin(self, event, plot):
         """Event handler for clicking to set origin of map.
 
         Parameters
         ----------
         event
             Click event.
-        currentEbsdMap : defdap.ebsd.Map
-            EBSD map to set origin for.
+        plot : defdap.plotting.MapPlot
+            Plot to capture clicks from.
 
         """
-        currentEbsdMap.origin = (int(event.ydata), int(event.xdata))
-        print("Origin set to ({:}, {:})".format(currentEbsdMap.origin[0], currentEbsdMap.origin[1]))
+        # check if click was on the map
+        if event.inaxes is not plot.ax:
+            return
+
+        origin = (int(event.xdata), int(event.ydata))
+        plot.callingMap.origin = origin
+        plot.addPoints([origin[0]], [origin[1]], updateLayer=0)
+        print(f"Origin set to ({origin[0]}, {origin[1]})")
 
-    def startLinking(self):
+    def start_linking(self):
         """Start interactive grain linking process of each EBSD map.
 
         """
-        for ebsdMap in self.ebsdMaps:
-            ebsdMap.locateGrainID(clickEvent=self.clickGrainGuess)
+        self.plots = []
+        for ebsd_map in self.ebsd_maps:
+            plot = ebsd_map.locateGrainID(clickEvent=self.click_grain_guess)
 
             # Add make link button to axes
-            btnAx = ebsdMap.fig.add_axes([0.8, 0.0, 0.1, 0.07])
-            Button(btnAx, 'Make link', color='0.85', hovercolor='0.95')
+            plot.addButton('Make link', self.make_link,
+                           color='0.85', hovercolor='0.95')
+
+            self.plots.append(plot)
 
-    def clickGrainGuess(self, event, currentEbsdMap):
+    def click_grain_guess(self, event, plot):
         """Guesses grain position in other maps, given click on one.
 
         Parameters
         ----------
         event
             Click handler.
-        currentEbsdMap : defdap.ebsd.Map
-            EBSD map that is clicked on.
+        plot : defdap.plotting.Plot
+            Plot to capture clicks from.
 
         """
-        # self is current linker instance even if run as click event handler from map class
-        if event.inaxes is currentEbsdMap.fig.axes[0]:
-            # axis 0 then is a click on the map
-
-            if currentEbsdMap is self.ebsdMaps[0]:
-                # clicked on 'master' map so highlight and guess grain on other maps
-                for ebsdMap in self.ebsdMaps:
-                    if ebsdMap is currentEbsdMap:
-                        # set current grain in ebsd map that clicked
-                        ebsdMap.clickGrainID(event)
-                    else:
-                        # Guess at grain in other maps
-                        # Calculated position relative to set origin of the map, scaled from step size of maps
-                        y0m = currentEbsdMap.origin[0]
-                        x0m = currentEbsdMap.origin[1]
-                        y0 = ebsdMap.origin[0]
-                        x0 = ebsdMap.origin[1]
-                        scaling = currentEbsdMap.stepSize / ebsdMap.stepSize
-
-                        x = int((event.xdata - x0m) * scaling + x0)
-                        y = int((event.ydata - y0m) * scaling + y0)
-
-                        ebsdMap.currGrainId = int(ebsdMap.grains[y, x]) - 1
-                        print(ebsdMap.currGrainId)
-
-                        # clear current axis and redraw euler map with highlighted grain overlay
-                        ebsdMap.ax.clear()
-                        ebsdMap.plotEulerMap(updateCurrent=True, highlightGrains=[ebsdMap.currGrainId])
-                        ebsdMap.fig.canvas.draw()
-            else:
-                # clicked on other map so correct guessed selected grain
-                currentEbsdMap.clickGrainID(event)
+        # check if click was on the map
+        if event.inaxes is not plot.ax:
+            return
+
+        curr_ebsd_map = plot.callingMap
+
+        if curr_ebsd_map is self.ebsd_maps[0]:
+            # clicked on 'master' map so highlight and guess grain on others
+
+            # set current grain in 'master' ebsd map
+            self.ebsd_maps[0].clickGrainID(event, plot, False)
+
+            # guess at grain in other maps
+            for ebsd_map, plot in zip(self.ebsd_maps[1:], self.plots[1:]):
+                # calculated position relative to set origin of the
+                # map, scaled from step size of maps
+                x0m = curr_ebsd_map.origin[0]
+                y0m = curr_ebsd_map.origin[1]
+                x0 = ebsd_map.origin[0]
+                y0 = ebsd_map.origin[1]
+                scaling = curr_ebsd_map.stepSize / ebsd_map.stepSize
 
-        elif event.inaxes is currentEbsdMap.fig.axes[1]:
-            # axis 1 then is a click on the button
-            self.makeLink()
+                x = int((event.xdata - x0m) * scaling + x0)
+                y = int((event.ydata - y0m) * scaling + y0)
 
-    def makeLink(self):
+                ebsd_map.currGrainId = int(ebsd_map.grains[y, x]) - 1
+                print(ebsd_map.currGrainId)
+
+                # update the grain highlights layer in the plot
+                plot.addGrainHighlights([ebsd_map.currGrainId],
+                                        alpha=ebsd_map.highlightAlpha)
+
+        else:
+            # clicked on other map so correct guessed selected grain
+            curr_ebsd_map.clickGrainID(event, plot, False)
+
+    def make_link(self, event, plot):
         """Make a link between the EBSD maps after clicking.
 
         """
         # create empty list for link
-        currLink = []
+        curr_link = []
 
-        for i, ebsdMap in enumerate(self.ebsdMaps):
-            if ebsdMap.currGrainId is not None:
-                currLink.append(ebsdMap.currGrainId)
+        for i, ebsd_map in enumerate(self.ebsd_maps):
+            if ebsd_map.currGrainId is not None:
+                curr_link.append(ebsd_map.currGrainId)
             else:
-                raise Exception("No grain setected in map {:d}.".format(i + 1))
-
-        self.links.append(tuple(currLink))
+                raise Exception(f"No grain setected in map {i + 1}.")
 
-        print("Link added " + str(tuple(currLink)))
+        curr_link = tuple(curr_link)
+        if curr_link not in self.links:
+            self.links.append(curr_link)
+            print("Link added " + str(curr_link))
 
-    def resetLinks(self):
+    def reset_links(self):
         """Reset links.
 
         """
@@ -1920,30 +1915,26 @@ def resetLinks(self):
 
 #   Analysis routines
 
-    def setAvOriFromInitial(self):
+    def set_ref_ori_from_master(self):
         """Loop over each map (not first/reference) and each link.
         Sets refOri of linked grains to refOri of grain in first map.
 
         """
-        masterMap = self.ebsdMaps[0]
-
-        for i, ebsdMap in enumerate(self.ebsdMaps[1:], start=1):
+        for i, ebsd_map in enumerate(self.ebsd_maps[1:], start=1):
             for link in self.links:
-                ebsdMap.grainList[link[i]].refOri = copy.deepcopy(masterMap.grainList[link[0]].refOri)
-
-        return
+                ebsd_map.grainList[link[i]].refOri = copy.deepcopy(
+                    self.ebsd_maps[0].grainList[link[0]].refOri
+                )
 
-    def updateMisOri(self, calcAxis=False):
+    def update_misori(self, calc_axis=False):
         """Recalculate misorientation for linked grain (not for first map)
 
         Parameters
         ----------
-        calcAxis : bool
+        calc_axis : bool
             Calculate the misorientation axis if True.
 
         """
-        for i, ebsdMap in enumerate(self.ebsdMaps[1:], start=1):
+        for i, ebsdMap in enumerate(self.ebsd_maps[1:], start=1):
             for link in self.links:
-                ebsdMap.grainList[link[i]].buildMisOriList(calcAxis=calcAxis)
-
-        return
+                ebsdMap.grainList[link[i]].buildMisOriList(calcAxis=calc_axis)
diff --git a/defdap/file_readers.py b/defdap/file_readers.py
index cc9edce..4cc2315 100644
--- a/defdap/file_readers.py
+++ b/defdap/file_readers.py
@@ -112,6 +112,7 @@ def parsePhase() -> Phase:
             phase = Phase(
                 lineSplit[2],
                 int(lineSplit[3]),
+                int(lineSplit[4]),
                 latticeParams
             )
             return phase
@@ -300,6 +301,7 @@ def parseLine(line: str, groupDict: Dict) -> None:
             self.loadedMetadata['phases'].append(Phase(
                 phaseMetadata['StructureName'],
                 int(phaseMetadata['LaueGroup']),
+                int(phaseMetadata['SpaceGroup']),
                 (
                     round(float(phaseMetadata['a']), 3),
                     round(float(phaseMetadata['b']), 3),
@@ -324,6 +326,7 @@ def parseLine(line: str, groupDict: Dict) -> None:
         self.checkMetadata()
 
         # Construct binary data format from listed fields
+        unknown_field_count = 0
         dataFormat = [('phase', 'uint8')]
         fieldLookup = {
             3: ('ph1', 'float32'),
@@ -342,7 +345,10 @@ def parseLine(line: str, groupDict: Dict) -> None:
                 fieldID = int(metadata['Fields']['Field{:}'.format(i + 1)])
                 dataFormat.append(fieldLookup[fieldID])
         except KeyError:
-            raise TypeError("Unknown data in EBSD file.")
+            print(f'\nUnknown field in file with key {fieldID}. '
+                  f'Assumming float32 data.')
+            unknown_field_count += 1
+            dataFormat.append((f'unknown_{unknown_field_count}', 'float32'))
 
         self.dataFormat = np.dtype(dataFormat)
 
diff --git a/defdap/hrdic.py b/defdap/hrdic.py
index c0a1aff..61d5c24 100755
--- a/defdap/hrdic.py
+++ b/defdap/hrdic.py
@@ -431,7 +431,7 @@ def setHomogPoint(self, points=None, display=None, **kwargs):
             # Call set homog points from base class setting the bin size
             super(type(self), self).setHomogPoint(binSize=binSize, points=points, **kwargs)
 
-    def linkEbsdMap(self, ebsdMap, transformType="affine", order=2):
+    def linkEbsdMap(self, ebsdMap, transformType="affine", **kwargs):
         """Calculates the transformation required to align EBSD dataset to DIC.
 
         Parameters
@@ -440,55 +440,39 @@ def linkEbsdMap(self, ebsdMap, transformType="affine", order=2):
             EBSD map object to link.
         transformType : str, optional
             affine, piecewiseAffine or polynomial.
-        order : int, optional
-            Order of polynomial transform to apply.
+        kwargs
+            All arguments are passed to `estimate` method of the transform.
 
         """
         self.ebsdMap = ebsdMap
+        calc_inv = False
         if transformType.lower() == "piecewiseaffine":
             self.ebsdTransform = tf.PiecewiseAffineTransform()
-            self.ebsdTransformInv = self.ebsdTransform.inverse
         elif transformType.lower() == "projective":
             self.ebsdTransform = tf.ProjectiveTransform()
-            self.ebsdTransformInv = self.ebsdTransform.inverse
         elif transformType.lower() == "polynomial":
+            calc_inv = True
             self.ebsdTransform = tf.PolynomialTransform()
-            # You can't calculate the inverse of a polynomial transform
-            # so have to estimate by swapping source and destination
-            # homog points
             self.ebsdTransformInv = tf.PolynomialTransform()
-            self.ebsdTransformInv.estimate(
-                np.array(self.ebsdMap.homogPoints),
-                np.array(self.homogPoints),
-                order=order
-            )
-            # calculate transform from EBSD to DIC frame
-            self.ebsdTransform.estimate(
-                np.array(self.homogPoints),
-                np.array(self.ebsdMap.homogPoints),
-                order=order
-            )
-            return
         else:
             # default to using affine
             self.ebsdTransform = tf.AffineTransform()
-            self.ebsdTransformInv = self.ebsdTransform.inverse
 
         # calculate transform from EBSD to DIC frame
         self.ebsdTransform.estimate(
             np.array(self.homogPoints),
-            np.array(self.ebsdMap.homogPoints)
+            np.array(self.ebsdMap.homogPoints),
+            **kwargs
         )
-
-        # Transform the EBSD boundaryLines to DIC reference frame
-        boundaryLineList = np.array(self.ebsdMap.boundaryLines).reshape(-1, 2)       # Flatten to coord list
-        boundaryLines = self.ebsdTransformInv(boundaryLineList).reshape(-1, 2, 2)    # Transform & reshape back
-        self.boundaryLines = np.round(boundaryLines - 0.5) + 0.5                         # Round to nearest
-
-        # Transform the EBSD phaseBoundaryLines to DIC reference frame
-        phaseBoundaryLineList = np.array(self.ebsdMap.phaseBoundaryLines).reshape(-1, 2)     # Flatten to coord list
-        phaseBoundaryLines = self.ebsdTransformInv(phaseBoundaryLineList).reshape(-1, 2, 2)  # Transform & reshape back
-        self.phaseBoundaryLines = np.round(phaseBoundaryLines - 0.5) + 0.5                   # Round to nearest
+        # Calculate inverse if required
+        if calc_inv:
+            self.ebsdTransformInv.estimate(
+                np.array(self.ebsdMap.homogPoints),
+                np.array(self.homogPoints),
+                **kwargs
+            )
+        else:
+            self.ebsdTransformInv = self.ebsdTransform.inverse
 
     def checkEbsdLinked(self):
         """Check if an EBSD map has been linked.
@@ -524,8 +508,8 @@ def warpToDicFrame(self, mapData, cropImage=True, order=1, preserve_range=False)
 
         Returns
         ----------
-        warpedMap
-            Map (i.e. EBSD map) warped to the DIC frame.
+        numpy.ndarray
+            Map (i.e. EBSD map data) warped to the DIC frame.
 
         """
         # Check a EBSD map is linked
@@ -557,11 +541,77 @@ def warpToDicFrame(self, mapData, cropImage=True, order=1, preserve_range=False)
                 order=order, preserve_range=preserve_range
             )
 
-        # return map
         return warpedMap
 
+    def warp_lines_to_dic_frame(self, lines):
+        """Warp a set of lines to the DIC reference frame.
+
+        Parameters
+        ----------
+        lines : list of tuples
+            Lines to warp. Each line is represented as a tuple of start
+            and end coordinates (x, y).
+
+        Returns
+        -------
+        list of tuples
+            List of warped lines with same representation as input.
+
+        """
+        # Flatten to coord list
+        lines = np.array(lines).reshape(-1, 2)
+        # Transform & reshape back
+        lines = self.ebsdTransformInv(lines_list).reshape(-1, 2, 2)
+        # Round to nearest
+        lines = np.round(lines - 0.5) + 0.5
+        lines = [(tuple(l[0]), tuple(l[1])) for l in lines]
+
+        return lines
+
+    @property
+    def boundaries(self):
+        """Returns EBSD map grain boundaries warped to DIC frame.
+
+        """
+        # Check a EBSD map is linked
+        self.checkEbsdLinked()
+
+        # image is returned cropped if a piecewise transform is being used
+        boundaries = self.ebsdMap.boundaries
+        boundaries = self.warpToDicFrame(-boundaries.astype(float),
+                                         cropImage=False)
+        boundaries = boundaries > 0.1
+
+        boundaries = mph.skeletonize(boundaries)
+        mph.remove_small_objects(boundaries, min_size=10, in_place=True,
+                                 connectivity=2)
+
+        # crop image if it is a simple affine transform
+        if type(self.ebsdTransform) is tf.AffineTransform:
+            # need to apply the translation of ebsd transform and
+            # remove 5% border
+            crop = np.copy(self.ebsdTransform.params[0:2, 2])
+            crop += 0.05 * np.array(self.ebsdMap.boundaries.shape)
+            # the crop is defined in EBSD coords so need to transform it
+            transformMatrix = np.copy(self.ebsdTransform.params[0:2, 0:2])
+            crop = np.matmul(np.linalg.inv(transformMatrix), crop)
+            crop = crop.round().astype(int)
+
+            boundaries = boundaries[crop[1]:crop[1] + self.yDim,
+                                    crop[0]:crop[0] + self.xDim]
+
+        return -boundaries.astype(int)
+
+    @property
+    def boundaryLines(self):
+        return self.warp_lines_to_dic_frame(self.ebsdMap.boundaryLines)
+
+    @property
+    def phaseBoundaryLines(self):
+        return self.warp_lines_to_dic_frame(self.ebsdMap.phaseBoundaryLines)
+
     def generateThresholdMask(self, mask, dilation=0, preview=True):
-        """ 
+        """
         Generate a dilated mask, based on a boolean array and previews the appication of
         this mask to the max shear map.
 
@@ -570,7 +620,7 @@ def generateThresholdMask(self, mask, dilation=0, preview=True):
         mask: numpy.array(bool)
             A boolean array where points to be removed are True
         dilation: int, optional
-            Number of pixels to dilate the mask by. Useful to remove anomalous points 
+            Number of pixels to dilate the mask by. Useful to remove anomalous points
             around masked values. No dilation applied if not specified.
         preview: bool
             If true, show the mask and preview the masked effective shear strain map.
@@ -581,43 +631,47 @@ def generateThresholdMask(self, mask, dilation=0, preview=True):
 
         >>> mask = dicMap.eMaxShear > 0.8
 
-        To remove data points in dicMap where e11 is above 1 or less than -1, use:        
-        
+        To remove data points in dicMap where e11 is above 1 or less than -1, use:
+
         >>> mask = (dicMap.e11 > 1) | (dicMap.e11 < -1)
 
         To remove data points in dicMap where corrVal is less than 0.4, use:
 
         >>> mask = dicMap.corrVal < 0.4
 
-        Note: correlation value data needs to be loaded seperately from the DIC map, 
+        Note: correlation value data needs to be loaded seperately from the DIC map,
         see :func:`defdap.hrdic.loadCorrValData`
 
         """
         self.mask = mask
 
-
         if dilation != 0:
             self.mask = binary_dilation(self.mask, iterations=dilation)
 
         numRemoved = np.sum(self.mask)
-        numTotal = self.xdim*self.ydim
+        numTotal = self.xdim * self.ydim
         numRemovedCrop = np.sum(self.crop(self.mask))
         numTotalCrop = self.xDim * self.yDim
 
         print('Filtering will remove {0} \ {1} ({2:.3f} %) datapoints in map'
-                    .format(numRemoved, numTotal,(numRemoved / numTotal)*100))
-        print('Filtering will remove {0} \ {1} ({2:.3f} %) datapoints in cropped map'
-                    .format(numRemovedCrop, numTotalCrop,(numRemovedCrop / numTotalCrop * 100)))
+              .format(numRemoved, numTotal, (numRemoved / numTotal) * 100))
+        print(
+            'Filtering will remove {0} \ {1} ({2:.3f} %) datapoints in cropped map'
+            .format(numRemovedCrop, numTotalCrop,
+                    (numRemovedCrop / numTotalCrop * 100)))
 
         if preview == True:
             plot1 = MapPlot.create(self, self.crop(self.mask), cmap='binary')
             plot1.setTitle('Removed datapoints in black')
-            plot2 = MapPlot.create(self, 
-                                  self.crop(np.where(self.mask == True, np.nan, self.eMaxShear)),
-                                  plotColourBar='True', 
-                                  clabel="Effective shear strain")
+            plot2 = MapPlot.create(self,
+                                   self.crop(
+                                       np.where(self.mask == True, np.nan,
+                                                self.eMaxShear)),
+                                   plotColourBar='True',
+                                   clabel="Effective shear strain")
             plot2.setTitle('Effective shear strain preview')
-        print('Use applyThresholdMask function to apply this filtering to data')
+        print(
+            'Use applyThresholdMask function to apply this filtering to data')
 
     def applyThresholdMask(self):
         """ Apply mask to all DIC map data by setting masked values to nan.
@@ -636,44 +690,11 @@ def applyThresholdMask(self):
         self.x_map = np.where(self.mask == True, np.nan, self.x_map)
         self.y_map = np.where(self.mask == True, np.nan, self.y_map)
 
-        self.component = {'f11': self.f11, 'f12': self.f12, 'f21': self.f21, 'f22': self.f22,
-                 'e11': self.e11, 'e12': self.e12, 'e22': self.e22,
-                 'eMaxShear': self.eMaxShear,
-                 'x_map': self.x_map, 'y_map': self.y_map}
-
-    @property
-    def boundaries(self):
-        """Returns EBSD map grain boundaries warped to DIC frame.
-
-        """
-        # Check a EBSD map is linked
-        self.checkEbsdLinked()
-
-        # image is returned cropped if a piecewise transform is being used
-        boundaries = self.ebsdMap.boundaries
-        boundaries = self.warpToDicFrame(-boundaries.astype(float),
-                                         cropImage=False)
-        boundaries = boundaries > 0.1
-
-        boundaries = mph.skeletonize(boundaries)
-        mph.remove_small_objects(boundaries, min_size=10, in_place=True,
-                                 connectivity=2)
-
-        # crop image if it is a simple affine transform
-        if type(self.ebsdTransform) is tf.AffineTransform:
-            # need to apply the translation of ebsd transform and
-            # remove 5% border
-            crop = np.copy(self.ebsdTransform.params[0:2, 2])
-            crop += 0.05 * np.array(self.ebsdMap.boundaries.shape)
-            # the crop is defined in EBSD coords so need to transform it
-            transformMatrix = np.copy(self.ebsdTransform.params[0:2, 0:2])
-            crop = np.matmul(np.linalg.inv(transformMatrix), crop)
-            crop = crop.round().astype(int)
-
-            boundaries = boundaries[crop[1]:crop[1] + self.yDim,
-                                    crop[0]:crop[0] + self.xDim]
-
-        return -boundaries.astype(int)
+        self.component = {'f11': self.f11, 'f12': self.f12, 'f21': self.f21,
+                          'f22': self.f22,
+                          'e11': self.e11, 'e12': self.e12, 'e22': self.e22,
+                          'eMaxShear': self.eMaxShear,
+                          'x_map': self.x_map, 'y_map': self.y_map}
 
     def setPatternPath(self, filePath, windowSize):
         """Set the path to the image of the pattern.
@@ -993,7 +1014,7 @@ def floodFill(self, x, y, grainIndex, points_left):
 
         return currentGrain
 
-    def runGrainInspector(self, vmax=0.1, corrAngle=None):
+    def runGrainInspector(self, vmax=0.1, corrAngle=0):
         """Run the grain inspector interactive tool.
 
         Parameters
diff --git a/defdap/inspector.py b/defdap/inspector.py
index 10e3a72..35babd7 100644
--- a/defdap/inspector.py
+++ b/defdap/inspector.py
@@ -35,8 +35,8 @@ class GrainInspector:
     """
     def __init__(self, 
         currMap: 'hrdic.Map', 
-        vmax: float = 0.1,
-        corrAngle: float = 0):
+        vmax: float,
+        corrAngle: float):
         # Initialise some values
         self.grainID = 0
         self.currMap = currMap
diff --git a/defdap/plotting.py b/defdap/plotting.py
index 3c66268..93e5e2d 100644
--- a/defdap/plotting.py
+++ b/defdap/plotting.py
@@ -260,17 +260,19 @@ def lineSlice(self, event, plot, action=None):
         >>> plot.addEventHandler('button_release_event', lambda e, p: lineSlice(e, p))
 
         """
+        # check if click was on the map
+        if event.inaxes is not self.ax:
+            return
 
-        if event.inaxes is self.ax:
-            if event.name == 'button_press_event':
-                self.p1 = (event.xdata, event.ydata)  # save 1st point
-            elif event.name == 'button_release_event':
-                self.p2 = (event.xdata, event.ydata)  # save 2nd point
-                self.addArrow(startEnd=(self.p1[0], self.p1[1], self.p2[0], self.p2[1]))
-                self.fig.canvas.draw_idle()
+        if event.name == 'button_press_event':
+            self.p1 = (event.xdata, event.ydata)  # save 1st point
+        elif event.name == 'button_release_event':
+            self.p2 = (event.xdata, event.ydata)  # save 2nd point
+            self.addArrow(startEnd=(self.p1[0], self.p1[1], self.p2[0], self.p2[1]))
+            self.fig.canvas.draw_idle()
 
-                if action is not None:
-                    action(plot=self, startEnd=(self.p1[0], self.p1[1], self.p2[0], self.p2[1]))
+            if action is not None:
+                action(plot=self, startEnd=(self.p1[0], self.p1[1], self.p2[0], self.p2[1]))
 
     @property
     def exists(self):
@@ -825,7 +827,7 @@ def addSlipTraces(self, topOnly=False, colours=None, pos=None, **kwargs):
         """
 
         if colours is None:
-            colours = self.callingGrain.ebsdMap.slipTraceColours
+            colours = self.callingGrain.ebsdGrain.phase.slipTraceColours
         slipTraceAngles = self.callingGrain.slipTraces
 
         self.addTraces(slipTraceAngles, colours, topOnly, pos=pos, **kwargs)
diff --git a/defdap/slip_systems/cubic_bcc.txt b/defdap/slip_systems/cubic_bcc.txt
index 58f1056..a852e6d 100644
--- a/defdap/slip_systems/cubic_bcc.txt
+++ b/defdap/slip_systems/cubic_bcc.txt
@@ -1,5 +1,5 @@
 # Plane normal then slip direction
-green, green, green, green, red, red, red, red, blue, blue, blue, blue
+green,green,green,green,red,red,red,red,blue,blue,blue,blue
 0	1	-1	1	1	1
 -1	0	1	1	1	1
 1	-1	0	1	1	1
@@ -25,25 +25,26 @@ green, green, green, green, red, red, red, red, blue, blue, blue, blue
 1	2	-1	-1	1	1
 1	-1	2	-1	1	1
 1	2	3	1	1	-1
+2	1	3	1	1	-1
+3	2	1	1	1	-1
+2	3	1	1	1	-1
+-1	3	2	1	1	-1
+3	-1	2	1	1	-1
 1	3	2	1	-1	1
+2	3	1	1	-1	1
+3	1	2	1	-1	1
+2	1	3	1	-1	1
+-1	2	3	1	-1	1
+3	2	-1	1	-1	1
 3	1	2	-1	1	1
 3	2	1	-1	1	1
-2	1	3	1	1	-1
-2	3	1	1	-1	1
+1 	2 	3	-1	1	1
+1	3	2	-1	1	1
+2	-1	3	-1	1	1
+2	3	-1	-1	1	1
 1	2	3	1	1	1
 1	3	2	1	1	1
 3	1	2	1	1	1
 2	1	3	1	1 	1
 2	3	1 	1	1	1
-1 	2 	3	-1	1	1
-1	3	2	-1	1	1
-3	1	2	1	-1	1
-3	2	1	1	1	-1
-2	1	3	1	-1	1
-2	3	1	1	1	-1
--1	2	3	1	-1	1
--1	3	2	1	1	-1
-3	-1	2	1	1	-1
-3	2	-1	1	-1	1
-2	-1	3	-1	1	1
-2	3	-1	-1	1	1
\ No newline at end of file
+3	1	2	1	1	1
\ No newline at end of file
diff --git a/defdap/slip_systems/cubic_bcc_110only.txt b/defdap/slip_systems/cubic_bcc_110only.txt
new file mode 100644
index 0000000..d2e62b6
--- /dev/null
+++ b/defdap/slip_systems/cubic_bcc_110only.txt
@@ -0,0 +1,14 @@
+# Plane normal then slip direction
+green,green,green,green
+0	1	-1	1	1	1
+-1	0	1	1	1	1
+1	-1	0	1	1	1
+0	1	1	1	1	-1
+-1	0	-1	1	1	-1
+1	-1	0	1	1	-1
+0	1	-1	-1	1	1
+1	0	1	-1	1	1
+-1	-1	0	-1	1	1
+0	-1	-1	1	-1	1
+-1	0	1	1	-1	1
+1	1	0	1	-1	1
\ No newline at end of file
diff --git a/defdap/slip_systems/cubic_fcc_damask.txt b/defdap/slip_systems/cubic_fcc_damask.txt
new file mode 100644
index 0000000..6794ea0
--- /dev/null
+++ b/defdap/slip_systems/cubic_fcc_damask.txt
@@ -0,0 +1,14 @@
+# Plane normal then slip direction
+blue,green,red,black
+ 1	 1	 1	 0	 1	-1
+ 1	 1	 1	-1	 0	 1
+ 1	 1	 1	 1	-1	 0
+-1	-1	 1	 0	-1	-1
+-1	-1	 1	 1	 0	 1
+-1	-1	 1	-1	 1	 0
+ 1	-1	-1	 0	-1	 1
+ 1	-1	-1	-1	 0	-1
+ 1	-1	-1	 1	 1	 0
+-1	 1	-1	 0	 1	 1
+-1	 1	-1	 1	 0	-1
+-1	 1	-1	-1	-1	 0
diff --git a/example_notebook.ipynb b/example_notebook.ipynb
index 6eb0631..2804257 100644
--- a/example_notebook.ipynb
+++ b/example_notebook.ipynb
@@ -127,7 +127,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "dicMap.printStatsTable(percentiles=[0, 50, 100], components = ['mss', 'e11', 'e22', 'e12'])"
+    "dicMap.printStatsTable(percentiles=[0, 50, 100], components = ['eMaxShear', 'e11', 'e22', 'e12'])"
    ]
   },
   {
@@ -153,7 +153,7 @@
    "metadata": {},
    "source": [
     "## Load in an EBSD map\n",
-    "Set the crystal symmetry (`cubic` in this case but `hexagonal` is also available) and (optially) load in a set of slip systems defined in files stored in the directory `defdap/slip_systems` by default (`cubic_fcc` is loaded here but `cubic_bcc` and `hexagonal` are other options). The orientation in the EBSD are converted to a quaternion representation so calculations can be applied later."
+    "The crystal structure of each phase is read from file and used to set the slip systems for the phase. The orientation in the EBSD are converted to a quaternion representation so calculations can be applied later."
    ]
   },
   {
@@ -165,7 +165,6 @@
     "ebsdFilePath = \"tests/data/testDataEBSD\"\n",
     "\n",
     "ebsdMap = ebsd.Map(ebsdFilePath)\n",
-    "ebsdMap.loadSlipSystems(\"cubic_fcc\")\n",
     "ebsdMap.buildQuatArray()"
    ]
   },
@@ -233,7 +232,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A list of the slip planes, colours and slip directions can be printed"
+    "A list of the slip planes, colours and slip directions can be printed for each phase in the map."
    ]
   },
   {
@@ -242,7 +241,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ebsdMap.printSlipSystems()"
+    "phase = ebsdMap.phases[0]\n",
+    "print(phase.name)\n",
+    "phase.printSlipSystems()"
    ]
   },
   {
@@ -981,7 +982,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -995,7 +996,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.9.7"
   }
  },
  "nbformat": 4,
diff --git a/setup.py b/setup.py
index 229be26..43c2542 100644
--- a/setup.py
+++ b/setup.py
@@ -61,8 +61,6 @@ def get_version():
         'peakutils',
         'matplotlib_scalebar',
         'networkx',
-        'IPython',
-        'jupyter'
     ],
     extras_require={
         'testing': ['pytest', 'coverage', 'pytest-cov', 'pytest_cases'],