From 9e11e12673e0eb12ee50324e8343306beee6b387 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Germ=C3=A1n=20Carrillo?= <gcarrillo@linuxmail.org>
Date: Thu, 23 May 2024 10:53:21 -0500
Subject: [PATCH] Switch from ellipsoidal to 2d cartesian coordinates, because
 that's what QGIS elevation profile tool supports. By doing so, we ensure
 users won't get displaced profile graphs (too much difference between
 ellipsoidal and cartesian distances) not displaced exported features (using
 3d cartesian distances, exported features won't match in the graph and will
 be displaced with respect to 2d distances calculated by the elevation profile
 tool)

---
 .../core/profiles/profile_generator.py          | 17 +++++++++++++++--
 swiss_locator/core/profiles/profile_results.py  | 16 +++++++++-------
 2 files changed, 24 insertions(+), 9 deletions(-)

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