diff --git a/CHANGELOG.md b/CHANGELOG.md
index f340443..2afadc4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,7 +9,19 @@
 ### Fixed
 
 
-## 0.93.1 (16-04-2020)
+## 0.93.3 (23-08-2021)
+
+### Added
+- Store EDX data as a dictionary called `EDX` accessible from the EBSD `Map` object
+- Add option to change IPF map background colour
+
+### Fixed
+- Fix bug with reading cpr EBSD file without EDX data
+- Fix issue with plotting Schmid factor maps
+- Fix bug with maps `component` not updating after masking
+
+
+## 0.93.2 (16-04-2021)
 
 ### Added
 - Reading of Channel project files that contain EDX data
@@ -19,7 +31,7 @@
 - Plotting lines with grain inspector
 
 
-## 0.93.1 (12-04-2020)
+## 0.93.1 (12-04-2021)
 
 ### Added
 - Started adding type hinting
@@ -45,7 +57,7 @@
 - Fix bug in 'warp' grain finding algorithm
 
 
-## 0.93.0 (20-02-2020)
+## 0.93.0 (20-02-2021)
 
 ### Added
 - Add EBSD file writing to CTF format.
diff --git a/defdap/_version.py b/defdap/_version.py
index ad880ac..327520f 100644
--- a/defdap/_version.py
+++ b/defdap/_version.py
@@ -1 +1 @@
-__version__ = '0.93.2'
+__version__ = '0.93.3'
diff --git a/defdap/ebsd.py b/defdap/ebsd.py
index 27ce8ea..87aee6e 100755
--- a/defdap/ebsd.py
+++ b/defdap/ebsd.py
@@ -71,8 +71,6 @@ class Map(base.Map):
         Map of misorientation axis components.
     kam : numpy.ndarray
         Map of KAM.
-    averageSchmidFactor : numpy.ndarray
-        Map of average Schmid factor.
     slipSystems : list of list of defdap.crystal.SlipSystem
         Slip systems grouped by slip plane.
     slipTraceColours list(str)
@@ -122,7 +120,6 @@ def __init__(self, fileName, dataType=None):
         self.misOri = None
         self.misOriAxis = None
         self.kam = None
-        self.averageSchmidFactor = None
         self.origin = (0, 0)
         self.GND = None
         self.Nye = None
@@ -167,6 +164,8 @@ def loadData(self, fileName, dataType=None):
         self.bandSlopeArray = dataDict['bandSlope']
         self.meanAngularDeviationArray = dataDict['meanAngularDeviation']
         self.phaseArray = dataDict['phase']
+        if int(metadataDict['EDX Windows']['Count']) > 0:
+            self.EDX = dataDict['EDXDict']
 
         # write final status
         yield "Loaded EBSD data (dimensions: {:} x {:} pixels, step " \
@@ -329,7 +328,7 @@ def plotEulerMap(self, phases=None, **kwargs):
 
         return MapPlot.create(self, map_colours, **plot_params)
 
-    def plotIPFMap(self, direction, phases=None, **kwargs):
+    def plotIPFMap(self, direction, backgroundColour = [0., 0., 0.], phases=None, **kwargs):
         """
         Plot a map with points coloured in IPF colouring,
         with respect to a given sample direction.
@@ -338,8 +337,10 @@ def plotIPFMap(self, direction, phases=None, **kwargs):
         ----------
         direction : np.array len 3
             Sample directiom.
+        backgroundColour : np.array len 3
+            Colour of background (i.e. for phases not plotted).
         phases : list of int
-            Which phases to plot for
+            Which phases to plot IPF data for.
         kwargs
             Other arguments passed to :func:`defdap.plotting.MapPlot.create`.
 
@@ -359,7 +360,7 @@ def plotIPFMap(self, direction, phases=None, **kwargs):
             phase_ids = phases
             phases = [self.phases[i] for i in phase_ids]
 
-        map_colours = np.zeros(self.shape + (3,))
+        map_colours = np.tile(np.array(backgroundColour), self.shape + (1,))
 
         for phase, phase_id in zip(phases, phase_ids):
             # calculate IPF colours for phase
@@ -1299,57 +1300,53 @@ def plotAverageGrainSchmidFactorsMap(self, planes=None, directions=None,
 
         """
         # Set default plot parameters then update with any input
-        plotParams = {
+        plot_params = {
             'vmin': 0,
             'vmax': 0.5,
             'cmap': 'gray',
             'plotColourBar': True,
             'clabel': "Schmid factor"
         }
-        plotParams.update(kwargs)
+        plot_params.update(kwargs)
 
         # Check that grains have been detected in the map
         self.checkGrainsDetected()
-        self.averageSchmidFactor = np.zeros([self.yDim, self.xDim])
 
         if self[0].averageSchmidFactors is None:
             raise Exception("Run 'calcAverageGrainSchmidFactors' first")
 
+        grains_sf = []
         for grain in self.grainList:
-            currentSchmidFactor = []
+            current_sf = []
 
             if planes is not None:
-                # Error catching
-                if np.max(planes) > len(self.slipSystems) - 1:
-                    raise Exception("Check plane IDs exists, IDs range from 0 "
-                                    "to {0}".format(len(self.slipSystems) - 1))
-
                 for plane in planes:
                     if directions is not None:
                         for direction in directions:
-                            currentSchmidFactor.append(grain.averageSchmidFactors[plane][direction])
+                            current_sf.append(
+                                grain.averageSchmidFactors[plane][direction]
+                            )
                     else:
-                        currentSchmidFactor.append(grain.averageSchmidFactors[plane])
-                # TODO: what is this doing?
-                currentSchmidFactor = [max(s) for s in zip(*currentSchmidFactor)]
+                        current_sf += grain.averageSchmidFactors[plane]
             else:
-                currentSchmidFactor = [max(s) for s in zip(*grain.averageSchmidFactors)]
+                for sf_group in grain.averageSchmidFactors:
+                    current_sf += sf_group
 
-            # Fill grain with colour
-            for coord in grain.coordList:
-                self.averageSchmidFactor[coord[1], coord[0]] = currentSchmidFactor[0]
+            grains_sf.append(current_sf)
 
-        self.averageSchmidFactor[self.averageSchmidFactor == 0] = 0.5
+        grains_sf = np.array(grains_sf)
+        grains_sf_max = np.max(grains_sf, axis=1)
 
-        plot = MapPlot.create(self, self.averageSchmidFactor, **plotParams)
+        plot = self.plotGrainDataMap(grainData=grains_sf_max, bg=0.5,
+                                     **plot_params)
 
         return plot
 
 
 class Grain(base.Grain):
     """
-    Class to encapsulate a grain in an EBSD map and useful analysis and plotting
-    methods.
+    Class to encapsulate a grain in an EBSD map and useful analysis and
+    plotting methods.
 
     Attributes
     ----------
@@ -1517,7 +1514,7 @@ def plotOriSpread(self, direction=np.array([0, 0, 1]), **kwargs):
         plotParams.update(kwargs)
         return Quat.plotIPF(self.quatList, direction, self.crystalSym,
                             **plotParams)
-                            
+
     def plotUnitCell(self, fig=None, ax=None, plot=None, **kwargs):
         """Plot an unit cell of the average grain orientation.
 
diff --git a/defdap/file_readers.py b/defdap/file_readers.py
index e3d1a0b..cc9edce 100644
--- a/defdap/file_readers.py
+++ b/defdap/file_readers.py
@@ -149,6 +149,9 @@ def parsePhase() -> Phase:
             *(np.array(acqEulers) * np.pi / 180)
         )
 
+        # TODO: Load EDX data from .ctf file, if it's accesible
+        self.loadedMetadata['EDX Windows'] = {'Count': int(0)}
+
         self.checkMetadata()
 
         # Construct data format from table header
@@ -306,15 +309,17 @@ def parseLine(line: str, groupDict: Dict) -> None:
                     round(float(phaseMetadata['gamma']), 3) * np.pi / 180
                 )
             ))
-
+ 
         # Deal with EDX data
+        edx_fields = {}
         if 'EDX Windows' in metadata:
-            edx_dict = metadata['EDX Windows']
-            num_edx = int(edx_dict['Count'])
+            self.loadedMetadata['EDX Windows'] = metadata['EDX Windows']
             edx_fields = {}
-            for i in range(1, num_edx + 1):
-                key = f"Window{i}"
-                edx_fields[100+i] = (f'EDX {edx_dict[key]}', 'float32')
+            for i in range(1, int(self.loadedMetadata['EDX Windows']['Count']) + 1):
+                name = self.loadedMetadata['EDX Windows'][f"Window{i}"]
+                edx_fields[100+i] = (f'EDX {name}', 'float32')
+        else:
+            self.loadedMetadata['EDX Windows'] = {'Count': int(0)}
 
         self.checkMetadata()
 
@@ -378,6 +383,14 @@ def loadOxfordCRC(self, fileName: str, fileDir: str = "") -> None:
         eulerAngles = np.reshape(
             binData[['ph1', 'phi', 'ph2']], (yDim, xDim)
         )
+
+        # Load EDX data into a dict
+        if int(self.loadedMetadata['EDX Windows']['Count']) > 0:
+            EDXFields = [key for key in binData.dtype.fields.keys() if key.startswith('EDX')]        
+            self.loadedData['EDXDict'] = dict(
+                [(field[4:], np.reshape(binData[field], (yDim, xDim))) for field in EDXFields]
+                )
+
         # flatten the structures so that the Euler angles are stored
         # into a normal array
         eulerAngles = np.array(eulerAngles.tolist()).transpose((2, 0, 1))
@@ -406,6 +419,7 @@ def load(self, dataDict: Dict[str, Any]) -> None:
         self.loadedMetadata['stepSize'] = dataDict['stepSize']
         assert type(dataDict['phases']) is list
         self.loadedMetadata['phases'] = dataDict['phases']
+        self.loadedMetadata['EDX Windows'] = {'Count': int(0)}
 
         self.checkMetadata()
 
diff --git a/defdap/hrdic.py b/defdap/hrdic.py
index e58c686..c0a1aff 100755
--- a/defdap/hrdic.py
+++ b/defdap/hrdic.py
@@ -175,7 +175,6 @@ def __init__(self, path, fname, dataType=None):
         # max shear component
         self.eMaxShear = np.sqrt(((self.e11 - self.e22) / 2.)**2 + self.e12**2)
 
-        # Dictionary references to all map types
         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,
@@ -185,7 +184,7 @@ def __init__(self, path, fname, dataType=None):
         self.cropDists = np.array(((0, 0), (0, 0)), dtype=int)
 
         self.plotDefault = lambda *args, **kwargs: self.plotMaxShear(plotGBs=True, *args, **kwargs)
-
+    
     @property
     def crystalSym(self):
         return self.ebsdMap.crystalSym
@@ -637,6 +636,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.