diff --git a/ogstools/meshlib/mesh_series.py b/ogstools/meshlib/mesh_series.py
index 055abca7145ca6b36426a2a41c20c1ae664fa2e6..dbbc86706a9bdcea3a0d30936a1d6c8743a27c2a 100644
--- a/ogstools/meshlib/mesh_series.py
+++ b/ogstools/meshlib/mesh_series.py
@@ -192,6 +192,7 @@ class MeshSeries:
     ):
         values = self.hdf5["meshes"][self.hdf5_bulk_name][data_name][:]
         geom = self.hdf5["meshes"][self.hdf5_bulk_name]["geometry"][0]
+        values = np.swapaxes(values, 0, 1)
 
         # remove flat dimensions for interpolation
         for index, axis in enumerate(geom.T):
@@ -202,11 +203,11 @@ class MeshSeries:
         if interp_method is None:
             interp_method = "linear"
         interp = {
-            "nearest": NearestNDInterpolator(geom, values.T),
-            "linear": LinearNDInterpolator(geom, values.T, np.nan),
+            "nearest": NearestNDInterpolator(geom, values),
+            "linear": LinearNDInterpolator(geom, values, np.nan),
         }[interp_method]
 
-        return interp(points).T
+        return np.swapaxes(interp(points), 0, 1)
 
     def probe(
         self,
diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index bf7bb19c6404243ad6dbdac899dda124d6105787..53dc33132a8aa34b800740af9540096dd62c926f 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -23,7 +23,7 @@ from ogstools.propertylib.unit_registry import u_reg
 from . import plot_features as pf
 from . import setup
 from .levels import get_levels
-from .utils import get_style_cycler, justified_labels
+from .utils import get_style_cycler
 
 # TODO: define default data_name for regions in setup
 
@@ -571,13 +571,13 @@ def plot_probe(
     else:
         fig = None
     ax.set_prop_cycle(get_style_cycler(len(points), colors, linestyles))
-
-    if labels is None:
-        j_labels = justified_labels(points)
-        labels = [f"{i}: {label}" for i, label in enumerate(j_labels)]
     ax.plot(times, values, label=labels, **kwargs)
-    ax.legend(facecolor="white", framealpha=1, prop={"family": "monospace"})
+    if labels is not None:
+        ax.legend(facecolor="white", framealpha=1, prop={"family": "monospace"})
     time_label = f"time / {time_unit}" if time_unit else "time"
+    ax.set_axisbelow(True)
+    ax.grid(which="major", color="lightgrey", linestyle="-")
+    ax.grid(which="minor", color="0.95", linestyle="--")
     unit_str = (
         f" / {mesh_property.get_output_unit()}"
         if mesh_property.get_output_unit()
@@ -586,7 +586,6 @@ def plot_probe(
     y_label = mesh_property.output_name.replace("_", " ") + unit_str
     ax.set_xlabel(time_label)
     ax.set_ylabel(y_label)
-    ax.grid(which="major", color="lightgrey", linestyle="-")
-    ax.grid(which="minor", color="0.95", linestyle="--")
+    ax.label_outer()
     ax.minorticks_on()
     return fig
diff --git a/ogstools/meshplotlib/utils.py b/ogstools/meshplotlib/utils.py
index 748732a3334b865efea37c5df0bac6f3fdb525f5..e31ce760c94b50f10df3f2ac01913d0ebc62718d 100644
--- a/ogstools/meshplotlib/utils.py
+++ b/ogstools/meshplotlib/utils.py
@@ -13,8 +13,9 @@ def justified_labels(points: np.ndarray) -> list[str]:
     col_lens = np.max(
         [[len(fmt(coord)) for coord in point] for point in points], axis=0
     )
+    dim = points.shape[1]
     return [
-        ",".join(fmt(point[i]).rjust(col_lens[i]) for i in range(3))
+        ",".join(fmt(point[i]).rjust(col_lens[i]) for i in range(dim))
         for point in points
     ]
 
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index 6104a16ee3c09b4a9941eb55ee5e62015842d571..d75bf6541ebd5fa4f4d4048595653033f89ab401 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -98,6 +98,11 @@ class MeshplotlibTest(unittest.TestCase):
         plot_probe(mesh_series, points, presets.temperature)
         points = mesh_series.read(0).points[[0, -1]]
         plot_probe(mesh_series, points, presets.temperature)
+        mesh_series = examples.meshseries_XDMF
+        points = mesh_series.read(0).center
+        plot_probe(mesh_series, points, presets.temperature)
+        mesh_property = presets.velocity.replace(data_name="darcy_velocity")
+        plot_probe(mesh_series, points, mesh_property)
 
     def test_animation(self):
         """Test creation of animation."""