diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index 428e3685f05443dd6812ff48ad34f6ffc31bb646..6ee8c65266ff756eb9af4d3a179b4c6de2ffa68e 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -2,6 +2,7 @@
 
 import types
 from copy import deepcopy
+from math import nextafter
 from typing import Literal, Optional, Union
 
 import numpy as np
@@ -74,31 +75,27 @@ def get_cmap_norm(
 ) -> tuple[mcolors.Colormap, mcolors.Normalize]:
     """Construct a discrete colormap and norm for the property field."""
     vmin, vmax = (levels[0], levels[-1])
-    bilinear = property.bilinear_cmap and vmin <= 0.0 <= vmax
-    cmap_str = setup.cmap_str(property)
-    if property.is_mask():
-        conti_cmap = mcolors.ListedColormap(cmap_str)
-    elif isinstance(cmap_str, list):
-        conti_cmap = [colormaps[c] for c in cmap_str]
-    else:
-        conti_cmap = colormaps[cmap_str]
-    if property.data_name == "temperature":
-        cool_colors = conti_cmap[0](np.linspace(0, 0.75, 128 * (vmin < 0)))
-        warm_colors = conti_cmap[1](np.linspace(0, 1, 128 * (vmax >= 0)))
-        conti_cmap = mcolors.LinearSegmentedColormap.from_list(
-            "temperature_cmap", np.vstack((cool_colors, warm_colors))
-        )
-    if bilinear:
-        vmin, vmax = np.max(np.abs([vmin, vmax])) * np.array([-1.0, 1.0])
     if property.categoric:
         vmin += 0.5
         vmax += 0.5
-    conti_norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
+
+    if property.bilinear_cmap:
+        if vmin <= 0.0 <= vmax:
+            vcenter = 0.0
+            vmin, vmax = np.max(np.abs([vmin, vmax])) * np.array([-1.0, 1.0])
+        else:
+            # only use one half of the diverging colormap
+            vcenter = nextafter(*sorted([vmin, vmax], reverse=(vmin <= 0.0)))
+        conti_norm = mcolors.TwoSlopeNorm(vcenter, vmin, vmax)
+    else:
+        conti_norm = mcolors.Normalize(vmin, vmax)
+    if isinstance(property.cmap, str):
+        continuous_cmap = colormaps[property.cmap]
+    else:
+        continuous_cmap = property.cmap
     mid_levels = np.append((levels[:-1] + levels[1:]) * 0.5, levels[-1])
-    colors = [conti_cmap(conti_norm(m_l)) for m_l in mid_levels]
+    colors = [continuous_cmap(conti_norm(m_l)) for m_l in mid_levels]
     cmap = mcolors.ListedColormap(colors, name="custom")
-    if setup.custom_cmap is not None:
-        cmap = setup.custom_cmap
     boundaries = get_level_boundaries(levels) if property.categoric else levels
     norm = mcolors.BoundaryNorm(
         boundaries=boundaries, ncolors=len(boundaries), clip=True
diff --git a/ogstools/meshplotlib/plot_setup.py b/ogstools/meshplotlib/plot_setup.py
index c50115be9626a5c769a691fba4c80d1189432015..d36e7f327061557b55cabcaa6cdfdbb35a62acac 100644
--- a/ogstools/meshplotlib/plot_setup.py
+++ b/ogstools/meshplotlib/plot_setup.py
@@ -3,9 +3,7 @@
 from dataclasses import dataclass
 from typing import Union
 
-from matplotlib.colors import Colormap
-
-from ogstools.propertylib.property import Property, Scalar
+from ogstools.propertylib.property import Scalar
 
 from .plot_setup_defaults import setup_dict
 
@@ -21,16 +19,6 @@ class PlotSetup:
 
     combined_colorbar: bool
     "True if all subplots share on colorbar, else each has its own colorbar."
-    custom_cmap: Colormap
-    "If provided, this colormap will be used for any plot."
-    cmap_dict_if_bilinear: dict
-    "A dictionary that maps bilinear colormaps to properties."
-    cmap_dict: dict
-    "A dictionary that maps colormaps to properties."
-    cmap_if_mask: list
-    "A list of colors corresponding to [True, False] values of masks."
-    default_cmap: str
-    "The default colormap to use."
     dpi: int
     "The resolution (dots per inch) for the figure."
     fig_scale: float
@@ -80,17 +68,6 @@ class PlotSetup:
     embedded_region_names_color: str
     "Color of the embedded region names inside the plot."
 
-    def cmap_str(self, property: Property) -> Union[str, list]:
-        """Get the colormap string for a given property."""
-        if property.is_mask():
-            return self.cmap_if_mask
-        if property.bilinear_cmap:
-            if property.data_name in self.cmap_dict_if_bilinear:
-                return self.cmap_dict_if_bilinear[property.data_name]
-        elif property.data_name in self.cmap_dict:
-            return self.cmap_dict[property.data_name]
-        return self.default_cmap
-
     @property
     def rcParams_scaled(self) -> dict:
         """Get the scaled rcParams values."""
@@ -126,11 +103,6 @@ class PlotSetup:
             length=Scalar("", obj["length"][0], obj["length"][1], ""),
             material_names=obj["material_names"],
             combined_colorbar=obj["combined_colorbar"],
-            custom_cmap=obj["custom_cmap"],
-            cmap_dict=obj["cmap_dict"],
-            cmap_dict_if_bilinear=obj["cmap_dict_if_bilinear"],
-            cmap_if_mask=obj["cmap_if_mask"],
-            default_cmap=obj["default_cmap"],
             rcParams=obj["rcParams"],
         )
 
diff --git a/ogstools/meshplotlib/plot_setup_defaults.py b/ogstools/meshplotlib/plot_setup_defaults.py
index 29810fa350d783dfc11ac55f1a018040c8d492c1..35f8e03cedacc604bc4c7dbb26ecb88ee97462bc 100644
--- a/ogstools/meshplotlib/plot_setup_defaults.py
+++ b/ogstools/meshplotlib/plot_setup_defaults.py
@@ -31,19 +31,6 @@ setup_dict = {
     "y_label": None,
     "log_scaled": False,
     "combined_colorbar": True,
-    "custom_cmap": None,
-    "cmap_dict": {
-        "displacement": "Greens",
-        "temperature": "plasma",
-        "pressure": "Blues",
-        "velocity": "coolwarm",
-        "MaterialIDs": "tab20",
-    },
-    "cmap_dict_if_bilinear": {
-        "displacement": "PRGn",
-        "temperature": ["Blues", "plasma"],
-    },
-    "cmap_if_mask": ["lightgrey", "green"],
     "rcParams": {
         "font.weight": "normal",
         "font.family": "sans-serif",
diff --git a/ogstools/propertylib/custom_colormaps.py b/ogstools/propertylib/custom_colormaps.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd101559fe1fed313f68af186a1384dfccafcf97
--- /dev/null
+++ b/ogstools/propertylib/custom_colormaps.py
@@ -0,0 +1,14 @@
+"Definitions of custom colormaps."
+
+import numpy as np
+from matplotlib import colormaps
+from matplotlib.colors import LinearSegmentedColormap as LSC
+
+mask_cmap = LSC.from_list("mask_cmap", ["lightgrey", "g"])
+temperature_cmap = LSC.from_list(
+    "temperature_cmap",
+    np.r_[
+        colormaps["Blues"](np.linspace(0, 0.75, 128, endpoint=True)),
+        colormaps["plasma"](np.linspace(0, 1, 128, endpoint=True)),
+    ],
+)
diff --git a/ogstools/propertylib/presets.py b/ogstools/propertylib/presets.py
index 93bc44b6bdf9a0a6dc56e7fe0b968415237a83a3..6bd7127284fc9317951a4b35ef1d58b2b562ba96 100644
--- a/ogstools/propertylib/presets.py
+++ b/ogstools/propertylib/presets.py
@@ -6,23 +6,25 @@ from typing import Optional
 from . import vector2scalar as v2s
 from .property import Matrix, Property, Scalar, Vector
 
+from .custom_colormaps import temperature_cmap
+
 T_mask = "temperature_active"
 H_mask = "pressure_active"
 M_mask = "displacement_active"
 
 # fmt: off
-material_id = Scalar("MaterialIDs", categoric=True)
-displacement = Vector("displacement", "m", "m", mask=M_mask)
+material_id = Scalar("MaterialIDs", categoric=True, cmap="tab20")
+displacement = Vector("displacement", "m", "m", mask=M_mask, cmap="PRGn")
 effective_pressure = Scalar("sigma", "Pa", "MPa", "effective_pressure", M_mask, v2s.effective_pressure)
 heatflowrate = Scalar("HeatFlowRate", mask=T_mask)
 massflowrate = Scalar("MassFlowRate", mask=H_mask)
 nodal_forces = Vector("NodalForces", mask=M_mask)
-pressure = Scalar("pressure", "Pa", "MPa", "pore_pressure", H_mask)
-hydraulic_height = Scalar("pressure", "m", "m", "hydraulic_height", H_mask)
+pressure = Scalar("pressure", "Pa", "MPa", "pore_pressure", H_mask, cmap="Blues")
+hydraulic_height = Scalar("pressure", "m", "m", "hydraulic_height", H_mask, cmap="Blues")
 qp_ratio = Scalar("sigma", "Pa", "percent", "QP_ratio", M_mask, v2s.qp_ratio)
 strain = Matrix("epsilon", "", "percent", "strain", M_mask)
 stress = Matrix("sigma", "Pa", "MPa", "stress", M_mask)
-temperature = Scalar("temperature", "K", "°C", mask=T_mask, bilinear_cmap=True)
+temperature = Scalar("temperature", "K", "°C", mask=T_mask, cmap=temperature_cmap, bilinear_cmap=True)
 velocity = Vector("velocity", "m/s", "m/s", "darcy_velocity", H_mask)
 von_mises_stress = Scalar("sigma", "Pa", "MPa", "von_Mises_stress", M_mask, v2s.von_mises)
 # fmt: on
diff --git a/ogstools/propertylib/property.py b/ogstools/propertylib/property.py
index 53d467db7a8bfa57b584d579bc1f96b84837d595..711c06da2d11e07fc439e0ee8cb0845e8c182008 100644
--- a/ogstools/propertylib/property.py
+++ b/ogstools/propertylib/property.py
@@ -9,8 +9,10 @@ from dataclasses import dataclass, replace
 from typing import Any, Callable, Union
 
 import numpy as np
+from matplotlib.colors import Colormap
 from pint.facets.plain import PlainQuantity
 
+from .custom_colormaps import mask_cmap
 from .unit_registry import u_reg
 from .utils import identity, sym_tensor_to_mat
 from .vector2scalar import trace
@@ -39,6 +41,8 @@ class Property:
         Callable[[Any], Any],
     ] = identity
     """The function to be applied on the data."""
+    cmap: Union[Colormap, str] = "coolwarm"
+    """Colormap to use for plotting."""
     bilinear_cmap: bool = False
     """Should this property be displayed with a bilinear cmap?"""
     categoric: bool = False
@@ -113,7 +117,9 @@ class Property:
         """
         :returns: A property representing this properties mask.
         """
-        return Property(data_name=self.mask, mask=self.mask, categoric=True)
+        return Property(
+            data_name=self.mask, mask=self.mask, categoric=True, cmap=mask_cmap
+        )
 
     @property
     def magnitude(self) -> "Property":
@@ -149,6 +155,7 @@ class Vector(Property):
             output_name=self.output_name + f"_{suffix[0 <= index <= 2]}",
             mask=self.mask,
             func=lambda x: np.array(x)[..., index],
+            cmap=self.cmap,
             bilinear_cmap=True,
         )
 
@@ -162,6 +169,7 @@ class Vector(Property):
             output_name=self.output_name + "_magnitude",
             mask=self.mask,
             func=lambda x: np.linalg.norm(x, axis=-1),
+            cmap=self.cmap,
         )
 
     @property
@@ -172,6 +180,7 @@ class Vector(Property):
             output_name=self.output_name + "_log10",
             mask=self.mask,
             func=lambda x: np.log10(np.linalg.norm(x, axis=-1)),
+            cmap=self.cmap,
         )
 
 
@@ -200,6 +209,7 @@ class Matrix(Property):
             mask=self.mask,
             func=lambda x: np.array(x)[..., index],
             bilinear_cmap=True,
+            cmap=self.cmap,
         )
 
     @property
@@ -212,6 +222,7 @@ class Matrix(Property):
             output_name=self.output_name + "_magnitude",
             mask=self.mask,
             func=lambda x: np.linalg.norm(sym_tensor_to_mat(x), axis=(-2, -1)),
+            cmap=self.cmap,
         )
 
     @property