Hi Sean Gilles,
I created a patch against Shapely-1.0.14 that exposes the GEOS linear
referencing methods GEOSProject and GEOSInterpolate. I also created a
basic test case. There is no error-checking though. Please tell me if
you need more information. Thanks for creating this package.
Index: tests/LinearReferencing.txt
===================================================================
--- tests/LinearReferencing.txt (revision 0)
+++ tests/LinearReferencing.txt (revision 0)
@@ -0,0 +1,33 @@
+Test linear projection and interpolation
+========================================
+
+ >>> from shapely.geometry import Point, LineString, MultiLineString
+ >>> from numpy import array
+ >>> point = Point(1, 1)
+ >>> line1 = LineString(([0, 0], [2, 0]))
+ >>> line2 = LineString(([3, 0], [3, 6]))
+ >>> multiLine = MultiLineString((array(line1), array(line2)))
+ >>> line1.project(point)
+ 1.0
+ >>> line1.projectNormalized(point)
+ 0.5
+ >>> line2.project(point)
+ 1.0
+ >>> line2.projectNormalized(point)
+ 0.16666666666666666
+ >>> multiLine.project(point)
+ 1.0
+ >>> multiLine.projectNormalized(point)
+ 0.125
+ >>> line1.interpolate(0.5).wkt
+ 'POINT (0.5000000000000000 0.0000000000000000)'
+ >>> line1.interpolateNormalized(0.5).wkt
+ 'POINT (1.0000000000000000 0.0000000000000000)'
+ >>> line2.interpolate(0.5).wkt
+ 'POINT (3.0000000000000000 0.5000000000000000)'
+ >>> line2.interpolateNormalized(0.5).wkt
+ 'POINT (3.0000000000000000 3.0000000000000000)'
+ >>> multiLine.interpolate(0.5).wkt
+ 'POINT (0.5000000000000000 0.0000000000000000)'
+ >>> multiLine.interpolateNormalized(0.5).wkt
+ 'POINT (3.0000000000000000 2.0000000000000000)'
Index: shapely/ctypes_declarations.py
===================================================================
--- shapely/ctypes_declarations.py (revision 1524)
+++ shapely/ctypes_declarations.py (working copy)
@@ -210,3 +210,14 @@
lgeos.GEOSDistance.restype = ctypes.c_int
lgeos.GEOSDistance.argtypes = [ctypes.c_void_p, ctypes.c_void_p,
ctypes.c_void_p]
+ lgeos.GEOSProject.restype = ctypes.c_double
+ lgeos.GEOSProject.argtypes = [ctypes.c_void_p, ctypes.c_void_p]
+
+ lgeos.GEOSProjectNormalized.restype = ctypes.c_double
+ lgeos.GEOSProjectNormalized.argtypes = [ctypes.c_void_p,
ctypes.c_void_p]
+
+ lgeos.GEOSInterpolate.restype = ctypes.c_void_p
+ lgeos.GEOSInterpolate.argtypes = [ctypes.c_void_p, ctypes.c_double]
+
+ lgeos.GEOSInterpolateNormalized.restype = ctypes.c_void_p
+ lgeos.GEOSInterpolateNormalized.argtypes = [ctypes.c_void_p,
ctypes.c_double]
Index: shapely/geometry/base.py
===================================================================
--- shapely/geometry/base.py (revision 1524)
+++ shapely/geometry/base.py (working copy)
@@ -396,6 +396,20 @@
retval = lgeos.GEOSDistance(self._geom, other._geom, byref(d))
return d.value
+ # Linear referencing
+
+ def project(self, other):
+ return lgeos.GEOSProject(self._geom, other._geom)
+
+ def projectNormalized(self, other):
+ return lgeos.GEOSProjectNormalized(self._geom, other._geom)
+
+ def interpolate(self, distance):
+ return geom_factory(lgeos.GEOSInterpolate(self._geom, distance))
+
+ def interpolateNormalized(self, distance):
+ return geom_factory(lgeos.GEOSInterpolateNormalized(self._geom,
distance))
+
# Topology operations
#
# These use descriptors to reduce the amount of boilerplate.
On 01/22/2010 11:16 AM, Roy Hyunjin Han wrote:
When I use novalis's patch from
http://code.djangoproject.com/attachment/ticket/11948/contrib-gis-linearref.patch
on Django 1.1.1, projection and interpolation seem to work properly,
including on MultiLineStrings for multiple line segments. I will try
to figure out why the same code isn't working in shapely.geometry.base.
from django.contrib.gis.geos import geometry
point = geometry.Point(1, 1)
line = geometry.LineString([[0, 0], [2, 0]])
assert line.project(point) == 1
print line.interpolate(line.project(point))
multiLine =
geometry.MultiLineString(geometry.LineString([[0,0],[2,0]]),
geometry.LineString([[3,0],[3,3]]))
assert multiLine.project(point) == 1
print multiLine.interpolate(multiLine.project(point))
On 01/22/2010 09:37 AM, Roy Hyunjin Han wrote:
I have been trying to use the GEOS linear referencing functions
through the base geometry class as Sean suggested in an earlier
email, and although interpolation works, projection does not
(lgeos.GEOSProject gives a return code of 0 instead of 1). Do you
have any idea why this is happening?
from shapely.geometry import base, LineString
from shapely.geos import lgeos
line = LineString([[0, 0], [2, 0]])
print '==== Interpolation ===='
point = base.geom_factory(lgeos.GEOSInterpolateNormalized(line._geom,
base.c_double(0.5)))
print point.wkt
print '==== Projection ===='
scalar = base.c_double()
returnCode = lgeos.GEOSProject(point._geom, line._geom,
base.byref(scalar))
if returnCode == 1:
print 'Projection of %s onto %s is %s' % (point.wkt, line.wkt,
scalar.value)
else:
print 'Projection failed'
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
Index: tests/LinearReferencing.txt
===================================================================
--- tests/LinearReferencing.txt (revision 0)
+++ tests/LinearReferencing.txt (revision 0)
@@ -0,0 +1,33 @@
+Test linear projection and interpolation
+========================================
+
+ >>> from shapely.geometry import Point, LineString, MultiLineString
+ >>> from numpy import array
+ >>> point = Point(1, 1)
+ >>> line1 = LineString(([0, 0], [2, 0]))
+ >>> line2 = LineString(([3, 0], [3, 6]))
+ >>> multiLine = MultiLineString((array(line1), array(line2)))
+ >>> line1.project(point)
+ 1.0
+ >>> line1.projectNormalized(point)
+ 0.5
+ >>> line2.project(point)
+ 1.0
+ >>> line2.projectNormalized(point)
+ 0.16666666666666666
+ >>> multiLine.project(point)
+ 1.0
+ >>> multiLine.projectNormalized(point)
+ 0.125
+ >>> line1.interpolate(0.5).wkt
+ 'POINT (0.5000000000000000 0.0000000000000000)'
+ >>> line1.interpolateNormalized(0.5).wkt
+ 'POINT (1.0000000000000000 0.0000000000000000)'
+ >>> line2.interpolate(0.5).wkt
+ 'POINT (3.0000000000000000 0.5000000000000000)'
+ >>> line2.interpolateNormalized(0.5).wkt
+ 'POINT (3.0000000000000000 3.0000000000000000)'
+ >>> multiLine.interpolate(0.5).wkt
+ 'POINT (0.5000000000000000 0.0000000000000000)'
+ >>> multiLine.interpolateNormalized(0.5).wkt
+ 'POINT (3.0000000000000000 2.0000000000000000)'
Index: shapely/ctypes_declarations.py
===================================================================
--- shapely/ctypes_declarations.py (revision 1524)
+++ shapely/ctypes_declarations.py (working copy)
@@ -210,3 +210,14 @@
lgeos.GEOSDistance.restype = ctypes.c_int
lgeos.GEOSDistance.argtypes = [ctypes.c_void_p, ctypes.c_void_p,
ctypes.c_void_p]
+ lgeos.GEOSProject.restype = ctypes.c_double
+ lgeos.GEOSProject.argtypes = [ctypes.c_void_p, ctypes.c_void_p]
+
+ lgeos.GEOSProjectNormalized.restype = ctypes.c_double
+ lgeos.GEOSProjectNormalized.argtypes = [ctypes.c_void_p, ctypes.c_void_p]
+
+ lgeos.GEOSInterpolate.restype = ctypes.c_void_p
+ lgeos.GEOSInterpolate.argtypes = [ctypes.c_void_p, ctypes.c_double]
+
+ lgeos.GEOSInterpolateNormalized.restype = ctypes.c_void_p
+ lgeos.GEOSInterpolateNormalized.argtypes = [ctypes.c_void_p,
ctypes.c_double]
Index: shapely/geometry/base.py
===================================================================
--- shapely/geometry/base.py (revision 1524)
+++ shapely/geometry/base.py (working copy)
@@ -396,6 +396,20 @@
retval = lgeos.GEOSDistance(self._geom, other._geom, byref(d))
return d.value
+ # Linear referencing
+
+ def project(self, other):
+ return lgeos.GEOSProject(self._geom, other._geom)
+
+ def projectNormalized(self, other):
+ return lgeos.GEOSProjectNormalized(self._geom, other._geom)
+
+ def interpolate(self, distance):
+ return geom_factory(lgeos.GEOSInterpolate(self._geom, distance))
+
+ def interpolateNormalized(self, distance):
+ return geom_factory(lgeos.GEOSInterpolateNormalized(self._geom,
distance))
+
# Topology operations
#
# These use descriptors to reduce the amount of boilerplate.
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
