Use length spheroid to full compatibility with Pseudo Mercator projection (3857) (#5038)

* refactor: use length spheroid

Author: submarcos

docs/changelog.rst

@@ -5,6 +5,10 @@ CHANGELOG
 2.120.0+dev     (XXXX-XX-XX)
 ----------------------------
 
+**Improvements**
+
+* Length are now computed with Earth spheroïd. We can now support using SRID 3857 (Pseudo mercator) to cover - almost - all Earth surface.
+
 
 2.120.0         (2025-10-22)
 ----------------------------

geotrek/altimetry/templates/altimetry/elevationinfo_fragment.html

@@ -6,7 +6,7 @@
     {# ``length`` attribute comes from AltimetryMixin and triggers #}
     <span title="{% trans "3D" %}">&#8605; {{ object.length|floatformat }} m</span>
     {# ``geom.length`` is 2D and is computed by GEOS #}
-    <span title="{% trans "2D" %}">(&#8594; {{ object.geom.length|floatformat }} m)</span></td>
+    <span title="{% trans "2D" %}">(&#8594; {{ object.length_2d|floatformat }} m)</span></td>
 </tr>
 {% endif %}
 {% if object.slope %}

geotrek/altimetry/tests/test_elevation.py

@@ -6,6 +6,7 @@
 from django.test import TestCase
 
 from geotrek.altimetry.helpers import AltimetryHelper
+from geotrek.common.tests.utils import LineStringInBounds
 from geotrek.core.models import Path, Topology
 from geotrek.core.tests.factories import TopologyFactory
 
@@ -338,7 +339,7 @@ def setUpTestData(cls):
         # Create a simple fake DEM
         with connection.cursor() as cur:
             cur.execute(
-                "INSERT INTO altimetry_dem (rast) VALUES (ST_MakeEmptyRaster(100, 125, 0, 125, 25, -25, 0, 0, %s))",
+                "INSERT INTO altimetry_dem (rast) VALUES (ST_MakeEmptyRaster(100, 125, 489353.59, 6587677.2, 25, -25, 0, 0, %s))",
                 [settings.SRID],
             )
             cur.execute("UPDATE altimetry_dem SET rast = ST_AddBand(rast, '16BSI')")
@@ -356,17 +357,17 @@ def setUpTestData(cls):
                         [x + 1, y + 1, demvalues[y][x]],
                     )
             cls.path = Path.objects.create(
-                geom=LineString((1, 101), (81, 101), (81, 99))
+                geom=LineStringInBounds((1, 101), (81, 101), (81, 99))
             )
 
     def test_2dlength_is_preserved(self):
-        self.assertEqual(self.path.geom_3d.length, self.path.geom.length)
+        self.assertAlmostEqual(self.path.geom.length, self.path.geom_3d.length)
 
     def test_3dlength(self):
         # before smoothing: (1 101 0, 21 101 0, 41 101 0, 61 101 3, 81 101 5, 81 99 15)
         # after smoothing:  (1 101 0, 21 101 0, 41 101 0, 61 101 1, 81 101 3, 81 99  9)
         # length: 20 + 20 + (20 ** 2 + 1) ** .5 + (20 ** 2 + 2 ** 2) ** .5 + (2 ** 2 + 6 ** 2) ** .5
-        self.assertEqual(round(self.path.length, 9), 83.127128724)
+        self.assertEqual(round(self.path.length, 9), 83.203851919)
 
 
 @skipIf(not settings.TREKKING_TOPOLOGY_ENABLED, "Test with dynamic segmentation only")

geotrek/api/v2/functions.py

@@ -22,7 +22,6 @@
 from rest_framework_gis.serializers import GeoFeatureModelSerializer
 
 from geotrek.api.v2.filters import get_published_filter_expression
-from geotrek.api.v2.functions import Length3D
 from geotrek.api.v2.mixins import (
     PDFSerializerMixin,
     PublishedRelatedObjectsSerializerMixin,
@@ -788,7 +787,7 @@ class PathSerializer(DynamicFieldsMixin, serializers.ModelSerializer):
         provider = serializers.SlugRelatedField(read_only=True, slug_field="name")
 
         def get_length_3d(self, obj):
-            return round(obj.length_3d_m, 1)
+            return round(obj.length, 1)
 
         class Meta:
             model = core_models.Path
@@ -936,7 +935,7 @@ def get_description_teaser(self, obj):
             )
 
         def get_length_3d(self, obj):
-            return round(obj.length_3d_m, 1)
+            return round(obj.length, 1)
 
         def get_gpx_url(self, obj):
             return build_url(
@@ -1132,7 +1131,6 @@ def get_steps(self, obj):
                 )
                 .annotate(
                     geom3d_transformed=Transform(F("geom_3d"), settings.API_SRID),
-                    length_3d_m=Length3D("geom_3d"),
                 )
             )
             FinalClass = override_serializer(

geotrek/api/v2/serializers.py

@@ -6,7 +6,6 @@
 from geotrek.api.v2 import filters as api_filters
 from geotrek.api.v2 import serializers as api_serializers
 from geotrek.api.v2 import viewsets as api_viewsets
-from geotrek.api.v2.functions import Length3D
 from geotrek.core import models as core_models
 
 
@@ -28,7 +27,6 @@ def get_queryset(self):
             .prefetch_related("usages", "networks")
             .annotate(
                 geom3d_transformed=Transform(F("geom_3d"), settings.API_SRID),
-                length_3d_m=Length3D("geom_3d"),
             )
             .order_by("pk")
         )  # Required for reliable pagination

geotrek/api/v2/views/core.py

@@ -12,7 +12,6 @@
 from geotrek.api.v2 import serializers as api_serializers
 from geotrek.api.v2 import viewsets as api_viewsets
 from geotrek.api.v2.decorators import cache_response_detail
-from geotrek.api.v2.functions import Length3D
 from geotrek.api.v2.renderers import SVGProfileRenderer
 from geotrek.common.models import AccessibilityAttachment, Attachment, HDViewPoint
 from geotrek.trekking import models as trekking_models
@@ -71,7 +70,6 @@ def get_queryset(self):
                 )
                 .annotate(
                     geom3d_transformed=Transform(F("geom_3d"), settings.API_SRID),
-                    length_3d_m=Length3D("geom_3d"),
                 )
                 .order_by("name")
             )  # Required for reliable pagination

geotrek/api/v2/views/trekking.py

@@ -9,6 +9,10 @@ class Length(GeoFunc):
     output_field = FloatField()
 
 
+class LengthSpheroid(GeoFunc):
+    output_field = FloatField()
+
+
 class SimplifyPreserveTopology(GeomOutputGeoFunc):
     """ST_SimplifyPreserveTopology postgis function"""
 

geotrek/common/functions.py

@@ -1,4 +1,6 @@
 import fiona
+from django.conf import settings
+from django.contrib.gis.geos import LineString, Point
 
 
 def update_gis(input_file_path: str, output_file_path: str, new_properties: dict):
@@ -22,3 +24,37 @@ def update_gis(input_file_path: str, output_file_path: str, new_properties: dict
                         properties={**feat.properties, **new_properties},
                     )
                 )
+
+
+def _get_y(y):
+    return 6587552.2 + y
+
+
+def _get_x(x):
+    return 489353.59 + x
+
+
+class LineStringInBounds(LineString):
+    """
+    A LineString that is guaranteed to be within the bounds of default projection.
+    Useful for tests that require geometries within certain bounds,
+    because ST_LENGTHSPHEROID can compute very bad data if coords is out of bounds.
+    """
+
+    def __init__(self, *args, **kwargs):
+        # Define a simple LineString within typical world bounds
+        kwargs.setdefault("srid", settings.SRID)  # Default SRID for French Lambert-93
+        coords = [[_get_x(point[0]), _get_y(point[1])] for point in args]
+        super().__init__(*coords, **kwargs)
+
+
+class PointInBounds(Point):
+    """
+    A Point that is guaranteed to be within the bounds of default projection.
+    Useful for tests that require geometries within certain bounds,
+    because ST_LENGTHSPHEROID can compute very bad data if coords is out of bounds.
+    """
+
+    def __init__(self, x=0, y=0, **kwargs):
+        kwargs.setdefault("srid", settings.SRID)
+        super().__init__((_get_x(x), _get_y(y)), **kwargs)

geotrek/common/tests/utils.py

@@ -1,6 +1,8 @@
 from django.contrib.gis.db import models
+from django.contrib.gis.db.models.functions import Transform
+from django.db.models import Value
 
-from geotrek.common.functions import Length
+from geotrek.common.functions import LengthSpheroid
 from geotrek.common.mixins.managers import NoDeleteManager
 
 
@@ -10,12 +12,7 @@ class PathManager(models.Manager):
 
     def get_queryset(self):
         """Hide all ``Path`` records that are not marked as visible."""
-        return (
-            super()
-            .get_queryset()
-            .filter(visible=True)
-            .annotate(length_2d=Length("geom"))
-        )
+        return super().get_queryset().filter(visible=True)
 
 
 class PathInvisibleManager(models.Manager):
@@ -30,7 +27,14 @@ class TopologyManager(NoDeleteManager):
     use_for_related_fields = True
 
     def get_queryset(self):
-        return super().get_queryset().annotate(length_2d=Length("geom"))
+        qs = super().get_queryset()
+        qs = qs.annotate(
+            length_2d=LengthSpheroid(
+                Transform("geom", 4326),
+                Value('SPHEROID["GRS_1980",6378137,298.257222101]'),
+            ),
+        )
+        return qs
 
 
 class PathAggregationManager(models.Manager):