diff --git a/swiss_locator/core/profiles/profile_generator.py b/swiss_locator/core/profiles/profile_generator.py
index fde4265..29b2749 100644
--- a/swiss_locator/core/profiles/profile_generator.py
+++ b/swiss_locator/core/profiles/profile_generator.py
@@ -11,6 +11,7 @@
     QgsCsException,
     QgsFeedback,
     QgsGeometry,
+    QgsGeometryUtils,
     QgsMessageLog,
     QgsNetworkAccessManager,
     QgsPoint,
@@ -96,16 +97,29 @@ def generateProfile(self, context):  # QgsProfileGenerationContext
             QgsMessageLog.logMessage(result["error"], "Swiss locator", Qgis.Critical)
             return False
 
+        cartesian_d = 0
         for point in result:
             if self.__feedback.isCanceled():
                 return False
 
+            # Note: d is ellipsoidal from the API
             x, y, z, d = self.__parse_response_point(point)
             point_z = QgsPoint(x, y, z)
             point_z.transform(self.__transformation, Qgis.TransformDirection.Reverse)
 
+            self.__results.ellipsoidal_distance_to_height[d] = z
+            self.__results.ellipsoidal_cross_section_geometries.append(QgsGeometry(QgsPoint(d, z)))
+
+            if d != 0:
+                # QGIS elevation profile won't calculate distances
+                # using 3d, so let's stick to 2d to avoid getting
+                # displaced markers or lines in the profile canvas
+                cartesian_d += QgsGeometryUtils.distance2D(point_z, self.__results.raw_points[-1])
+
             self.__results.raw_points.append(point_z)
-            self.__results.distance_to_height[d] = z
+            self.__results.cartesian_distance_to_height[cartesian_d] = z
+            self.__results.cartesian_cross_section_geometries.append(QgsGeometry(QgsPoint(cartesian_d, z)))
+
             if z < self.__results.min_z:
                 self.__results.min_z = z
 
@@ -113,7 +127,6 @@ def generateProfile(self, context):  # QgsProfileGenerationContext
                 self.__results.max_z = z
 
             self.__results.geometries.append(QgsGeometry(point_z))
-            self.__results.cross_section_geometries.append(QgsGeometry(QgsPoint(d, z)))
 
         return not self.__feedback.isCanceled()
 
diff --git a/swiss_locator/core/profiles/profile_results.py b/swiss_locator/core/profiles/profile_results.py
index 5679415..2790b50 100644
--- a/swiss_locator/core/profiles/profile_results.py
+++ b/swiss_locator/core/profiles/profile_results.py
@@ -22,9 +22,11 @@ def __init__(self):
 
         self.__profile_curve = None
         self.raw_points = []  # QgsPointSequence
-        self.distance_to_height = {}
+        self.cartesian_distance_to_height = {}
+        self.ellipsoidal_distance_to_height = {}
         self.geometries = []
-        self.cross_section_geometries = []
+        self.cartesian_cross_section_geometries = []
+        self.ellipsoidal_cross_section_geometries = []
         self.min_z = 4500
         self.max_z = -100
 
@@ -55,7 +57,7 @@ def asFeatures(self, type, feedback):
                 result.append(feature)
 
         elif type == Qgis.ProfileExportType.Profile2D:
-            for geom in self.cross_section_geometries:
+            for geom in self.cartesian_cross_section_geometries:
                 feature = QgsAbstractProfileResults.Feature()
                 feature.geometry = geom
                 feature.layerIdentifier = self.type()
@@ -70,7 +72,7 @@ def asFeatures(self, type, feedback):
                 # Since we've got distance/elevation pairs as
                 # x,y for cross-section geometries, and since
                 # both point arrays have the same length:
-                p = self.cross_section_geometries[i].asPoint()
+                p = self.cartesian_cross_section_geometries[i].asPoint()
                 feature.attributes = {"distance": p.x(), "elevation": p.y()}
                 result.append(feature)
 
@@ -93,7 +95,7 @@ def snapPoint(self, point, context):
 
         prev_distance = float('inf')
         prev_elevation = 0
-        for k, v in self.distance_to_height.items():
+        for k, v in self.cartesian_distance_to_height.items():
             # find segment which corresponds to the given distance along curve
             if k != 0 and prev_distance <= point.distance() <= k:
                 dx = k - prev_distance
@@ -159,7 +161,7 @@ def check_line(
         current_line = QPolygonF()
         prev_distance = None
         current_part_start_distance = 0
-        for k, v in self.distance_to_height.items():
+        for k, v in self.cartesian_distance_to_height.items():
             if not len(current_line):  # new part
                 if not v:  # skip emptiness
                     continue
@@ -200,7 +202,7 @@ def __render_markers(self, context):
 
         self.marker_symbol.startRender(context.renderContext())
 
-        for k, v in self.distance_to_height.items():
+        for k, v in self.cartesian_distance_to_height.items():
             if not v:
                 continue