diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py
index 233be569a..5a13f440d 100644
--- a/python/damask/_colormap.py
+++ b/python/damask/_colormap.py
@@ -6,6 +6,7 @@ from pathlib import Path
 from typing import Sequence, Union, TextIO
 
 import numpy as np
+import scipy.interpolate as interp
 import matplotlib as mpl
 if os.name == 'posix' and 'DISPLAY' not in os.environ:
     mpl.use('Agg')
@@ -41,7 +42,7 @@ class Colormap(mpl.colors.ListedColormap):
     https://doi.org/10.1016/j.ijplas.2012.09.012
 
     Matplotlib colormaps overview
-    https://matplotlib.org/tutorials/colors/colormaps.html
+    https://matplotlib.org/stable/tutorials/colors/colormaps.html
 
     """
 
@@ -191,6 +192,37 @@ class Colormap(mpl.colors.ListedColormap):
             return Colormap.from_range(definition['low'],definition['high'],name,N)
 
 
+    def at(self,
+           fraction : Union[float,Sequence[float]]) -> np.ndarray:
+        """
+        Interpolate color at fraction.
+
+        Parameters
+        ----------
+        fraction : float or sequence of float
+            Fractional coordinate(s) to evaluate Colormap at.
+
+        Returns
+        -------
+        color : np.ndarray, shape(...,4)
+            RGBA values of interpolated color(s).
+
+        Examples
+        --------
+        >>> import damask
+        >>> cmap = damask.Colormap.from_predefined('gray')
+        >>> cmap.at(0.5)
+        array([0.5, 0.5, 0.5, 1. ])
+        >>> 'rgb({},{},{})'.format(*cmap.at(0.5))
+        'rgb(0.5,0.5,0.5)'
+
+        """
+        return interp.interp1d(np.linspace(0,1,self.N),
+                               self.colors,
+                               axis=0,
+                               assume_sorted=True)(fraction)
+
+
     def shade(self,
               field: np.ndarray,
               bounds: Sequence[float] = None,
@@ -213,7 +245,6 @@ class Colormap(mpl.colors.ListedColormap):
             RGBA image of shaded data.
 
         """
-        N = len(self.colors)
         mask = np.logical_not(np.isnan(field) if gap is None else \
                np.logical_or (np.isnan(field), field == gap))                                       # mask NaN (and gap if present)
 
@@ -227,7 +258,7 @@ class Colormap(mpl.colors.ListedColormap):
 
         return Image.fromarray(
             (np.dstack((
-                        self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(N-1))).astype(np.uint16),:3],
+                        self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(self.N-1))).astype(np.uint16),:3],
                         mask.astype(float)
                        )
                       )*255
@@ -343,7 +374,7 @@ class Colormap(mpl.colors.ListedColormap):
         # ToDo: test in GOM
         GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \
                 +  '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \
-                + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(self.colors)}' \
+                + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {self.N}' \
                 + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
                 + '\n'
 
diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py
index ab9bcf92f..156907670 100644
--- a/python/tests/test_Colormap.py
+++ b/python/tests/test_Colormap.py
@@ -139,6 +139,16 @@ class TestColormap:
         c += c
         assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:]))
 
+    @pytest.mark.parametrize('N,cmap,at,result',[
+           (8,'gray',0.5,[0.5,0.5,0.5]),
+           (17,'gray',0.5,[0.5,0.5,0.5]),
+           (17,'gray',[0.5,0.75],[[0.5,0.5,0.5],[0.75,0.75,0.75]]),
+           ])
+    def test_at_value(self, N, cmap, at, result):
+        assert np.allclose(Colormap.from_predefined(cmap,N=N).at(at)[...,:3],
+                           result,
+                           rtol=0.005)
+
     @pytest.mark.parametrize('bounds',[None,[2,10]])
     def test_shade(self,ref_path,update,bounds):
         data = np.add(*np.indices((10, 11)))