Fixed #11948 -- Added interpolate and project linear referencing methods

Thanks novalis for the report and the initial patch, and Anssi
Kääriäinen and Justin Bronn for the review.

django/contrib/gis/geos/geometry.py

@@ -581,6 +581,20 @@ class GEOSGeometry(GEOSBase, ListMixin):
         "Return the envelope for this geometry (a polygon)."
         return self._topology(capi.geos_envelope(self.ptr))
 
+    def interpolate(self, distance):
+        if not isinstance(self, (LineString, MultiLineString)):
+            raise TypeError('interpolate only works on LineString and MultiLineString geometries')
+        if not hasattr(capi, 'geos_interpolate'):
+            raise NotImplementedError('interpolate requires GEOS 3.2+')
+        return self._topology(capi.geos_interpolate(self.ptr, distance))
+
+    def interpolate_normalized(self, distance):
+        if not isinstance(self, (LineString, MultiLineString)):
+            raise TypeError('interpolate only works on LineString and MultiLineString geometries')
+        if not hasattr(capi, 'geos_interpolate_normalized'):
+            raise NotImplementedError('interpolate_normalized requires GEOS 3.2+')
+        return self._topology(capi.geos_interpolate_normalized(self.ptr, distance))
+
     def intersection(self, other):
         "Returns a Geometry representing the points shared by this Geometry and other."
         return self._topology(capi.geos_intersection(self.ptr, other.ptr))
@@ -590,6 +604,24 @@ class GEOSGeometry(GEOSBase, ListMixin):
         "Computes an interior point of this Geometry."
         return self._topology(capi.geos_pointonsurface(self.ptr))
 
+    def project(self, point):
+        if not isinstance(point, Point):
+            raise TypeError('locate_point argument must be a Point')
+        if not isinstance(self, (LineString, MultiLineString)):
+            raise TypeError('locate_point only works on LineString and MultiLineString geometries')
+        if not hasattr(capi, 'geos_project'):
+            raise NotImplementedError('geos_project requires GEOS 3.2+')
+        return capi.geos_project(self.ptr, point.ptr)
+
+    def project_normalized(self, point):
+        if not isinstance(point, Point):
+            raise TypeError('locate_point argument must be a Point')
+        if not isinstance(self, (LineString, MultiLineString)):
+            raise TypeError('locate_point only works on LineString and MultiLineString geometries')
+        if not hasattr(capi, 'geos_project_normalized'):
+            raise NotImplementedError('project_normalized requires GEOS 3.2+')
+        return capi.geos_project_normalized(self.ptr, point.ptr)
+
     def relate(self, other):
         "Returns the DE-9IM intersection matrix for this Geometry and the other."
         return capi.geos_relate(self.ptr, other.ptr).decode()

django/contrib/gis/geos/prototypes/topology.py

@@ -8,18 +8,18 @@ __all__ = ['geos_boundary', 'geos_buffer', 'geos_centroid', 'geos_convexhull',
            'geos_simplify', 'geos_symdifference', 'geos_union', 'geos_relate']
 
 from ctypes import c_double, c_int
-from django.contrib.gis.geos.libgeos import GEOM_PTR, GEOS_PREPARE
-from django.contrib.gis.geos.prototypes.errcheck import check_geom, check_string
+from django.contrib.gis.geos.libgeos import geos_version_info, GEOM_PTR, GEOS_PREPARE
+from django.contrib.gis.geos.prototypes.errcheck import check_geom, check_minus_one, check_string
 from django.contrib.gis.geos.prototypes.geom import geos_char_p
 from django.contrib.gis.geos.prototypes.threadsafe import GEOSFunc
 
-def topology(func, *args):
+def topology(func, *args, **kwargs):
     "For GEOS unary topology functions."
     argtypes = [GEOM_PTR]
     if args: argtypes += args
     func.argtypes = argtypes
-    func.restype = GEOM_PTR
-    func.errcheck = check_geom
+    func.restype = kwargs.get('restype', GEOM_PTR)
+    func.errcheck = kwargs.get('errcheck', check_geom)
     return func
 
 ### Topology Routines ###
@@ -49,3 +49,16 @@ if GEOS_PREPARE:
     geos_cascaded_union.argtypes = [GEOM_PTR]
     geos_cascaded_union.restype = GEOM_PTR
     __all__.append('geos_cascaded_union')
+
+# Linear referencing routines
+info = geos_version_info()
+if info['version'] >= '3.2.0':
+    geos_project = topology(GEOSFunc('GEOSProject'), GEOM_PTR,
+        restype=c_double, errcheck=check_minus_one)
+    geos_interpolate = topology(GEOSFunc('GEOSInterpolate'), c_double)
+
+    geos_project_normalized = topology(GEOSFunc('GEOSProjectNormalized'),
+        GEOM_PTR, restype=c_double, errcheck=check_minus_one)
+    geos_interpolate_normalized = topology(GEOSFunc('GEOSInterpolateNormalized'), c_double)
+    __all__.extend(['geos_project', 'geos_interpolate',
+        'geos_project_normalized', 'geos_interpolate_normalized'])

django/contrib/gis/geos/tests/test_geos.py

@@ -1023,6 +1023,27 @@ class GEOSTest(unittest.TestCase, TestDataMixin):
 
         print("\nEND - expecting GEOS_NOTICE; safe to ignore.\n")
 
+    @unittest.skipUnless(geos_version_info()['version'] >= '3.2.0', "geos >= 3.2.0 is required")
+    def test_linearref(self):
+        "Testing linear referencing"
+
+        ls = fromstr('LINESTRING(0 0, 0 10, 10 10, 10 0)')
+        mls = fromstr('MULTILINESTRING((0 0, 0 10), (10 0, 10 10))')
+
+        self.assertEqual(ls.project(Point(0, 20)), 10.0)
+        self.assertEqual(ls.project(Point(7, 6)), 24)
+        self.assertEqual(ls.project_normalized(Point(0, 20)), 1.0/3)
+
+        self.assertEqual(ls.interpolate(10), Point(0, 10))
+        self.assertEqual(ls.interpolate(24), Point(10, 6))
+        self.assertEqual(ls.interpolate_normalized(1.0/3), Point(0, 10))
+
+        self.assertEqual(mls.project(Point(0, 20)), 10)
+        self.assertEqual(mls.project(Point(7, 6)), 16)
+
+        self.assertEqual(mls.interpolate(9), Point(0, 9))
+        self.assertEqual(mls.interpolate(17), Point(10, 7))
+
     def test_geos_version(self):
         "Testing the GEOS version regular expression."
         from django.contrib.gis.geos.libgeos import version_regex

docs/ref/contrib/gis/geos.txt

@@ -416,11 +416,36 @@ quarter circle (defaults is 8).
 Returns a :class:`GEOSGeometry` representing the points making up this
 geometry that do not make up other.
 
+.. method:: GEOSGeometry.interpolate(distance)
+.. method:: GEOSGeometry.interpolate_normalized(distance)
+
+.. versionadded:: 1.5
+
+Given a distance (float), returns the point (or closest point) within the
+geometry (:class:`LineString` or :class:`MultiLineString`) at that distance.
+The normalized version takes the distance as a float between 0 (origin) and 1
+(endpoint).
+
+Reverse of :meth:`GEOSGeometry.project`.
+
 .. method:: GEOSGeometry:intersection(other)
 
 Returns a :class:`GEOSGeometry` representing the points shared by this
 geometry and other.
 
+.. method:: GEOSGeometry.project(point)
+.. method:: GEOSGeometry.project_normalized(point)
+
+.. versionadded:: 1.5
+
+Returns the distance (float) from the origin of the geometry
+(:class:`LineString` or :class:`MultiLineString`) to the point projected on the
+geometry (that is to a point of the line the closest to the given point).
+The normalized version returns the distance as a float between 0 (origin) and 1
+(endpoint).
+
+Reverse of :meth:`GEOSGeometry.interpolate`.
+
 .. method:: GEOSGeometry.relate(other)
 
 Returns the DE-9IM intersection matrix (a string) representing the

docs/releases/1.5.txt

@@ -103,10 +103,22 @@ associated with proxy models.
 
 New ``view`` variable in class-based views context
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 In all :doc:`generic class-based views </topics/class-based-views/index>`
 (or any class-based view inheriting from ``ContextMixin``), the context dictionary
 contains a ``view`` variable that points to the ``View`` instance.
 
+GeoDjango
+~~~~~~~~~
+
+* :class:`~django.contrib.gis.geos.LineString` and
+  :class:`~django.contrib.gis.geos.MultiLineString` GEOS objects now support the
+  :meth:`~django.contrib.gis.geos.GEOSGeometry.interpolate()` and
+  :meth:`~django.contrib.gis.geos.GEOSGeometry.project()` methods
+  (so-called linear referencing).
+
+* Support for GDAL < 1.5 has been dropped.
+
 Minor features
 ~~~~~~~~~~~~~~
 
@@ -379,8 +391,6 @@ on the form.
 Miscellaneous
 ~~~~~~~~~~~~~
 
-* GeoDjango dropped support for GDAL < 1.5
-
 * :func:`~django.utils.http.int_to_base36` properly raises a :exc:`TypeError`
   instead of :exc:`ValueError` for non-integer inputs.