geoalg / geoalg /

#!/usr/bin/env python
# Author:  mozman
# Purpose: polylines
# Created: 31.03.2010

from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import

class Polyline(list):
    def __init__(self, points=[], closed=False):
        super(Polyline, self).__init__(points)
        self.closed = closed

    def nsegments(self):
        if self.closed:
            return len(self)
            return len(self)-1

    def length(self):
        plen = 0.
        for index in range(self.nsegments):
            plen += self.segment_length(index)
        return plen

    def point_at(self, distance):
        """Get point with <distance> from polyline start-point."""
        plen = 0.
        index = 0
        while plen < distance:
                plen += self.segment_length(index)
            except IndexError:
                return None
            index += 1
        if index == len(self):
            start = self[-1]
            end = self[0]
            start = self[index-1]
            end = self[index]
        return _interpolate(start, end, plen-distance)

    def segment_length(self, index):
        seg = self.segment(index)
        return _dist(seg[0], seg[1])

    def segment(self, index):
            return (self[index], self[index+1])
        except IndexError:
            # for closed polylines exists an extra segment from tail to head
            # with the index == len(self)
            if self.closed:
                return (self[index], self[0])
                raise # reraise IndexError

    def itersegments(self):
        """Iterate polyline segments as (startpoint, endpoint) tuples."""
        for index in range(self.nsegments):
            yield self.segment(index)

def _dist(p1, p2):
    return sum( [(c1-c2)**2 for (c1, c2) in zip(p1, p2)] )**0.5

def _interpolate(start, end, dist):
    """Interpolate point with distance <dist> from start."""
    def ip_2d():
        dx = end[0] - start[0]
        dy = end[1] - start[1]
        factor = dist / (dx**2 + dy**2)**0.5
        return (start[0]+dx*factor, start[1]+dy*factor)

    def ip_3d():
        dx = end[0] - start[0]
        dy = end[1] - start[1]
        dz = end[2] - start[2]
        factor = dist / (dx**2 + dy**2 + dz**2)**0.5
        return (start[0]+dx*factor, start[1]+dy*factor, start[2]+dz*factor)

    return (ip_2d() if len(start) == 2 else ip_3d())