tin_reducer_gml

import time
from functools import partial
from itertools import chain, imap
from osgeo import ogr
from pywps.Process.Process import WPSProcess


def iter_points(layer):
    for feature in layer:
        exterior_ring = feature.GetGeometryRef().GetGeometryRef(0)
        for i in xrange(exterior_ring.GetPointCount()):
            yield exterior_ring.GetPoint(i)


def clean_points(points):
    #
    # Delete the last point of each triangle (we only need three point for our
    # calculations)
    #
    # TODO: The code below has the effect of the original code but it doesn't
    #       match the description above and it repeats points which will get
    #       filtered out in a later processing step anyway.  I guess the
    #       `drop_every_nth()` function below is actually wanted here instead.
    #
    points = iter(points)
    try:
        while True:
            yield next(points)
            next(points)
            point = next(points)
            yield point
            yield point
            next(points)
    except StopIteration:
        pass  # Ignored intentionally.


def drop_every_nth(iterable, nth):
    return (item for i, item in enumerate(iterable) if i % nth != 0)


def calculate_mid_value(a, b):
    return min(a, b) + abs(b - a) / 2.0


def partition_points(points, bounding_box):
    bbox_xmin, bbox_xmax, bbox_ymin, bbox_ymax = bounding_box
    vertical_middle = calculate_mid_value(bbox_xmin, bbox_xmax)
    horizontal_middle = calculate_mid_value(bbox_ymin, bbox_ymax)
    result = [[list(), list()], [list(), list()]]
    for point in points:
        x, y, _ = point
        #
        # TODO: Points *on* the lines won't get filtered out here.  If that's
        #       a bug we need an extra test here.
        #
        result[y < horizontal_middle][x < vertical_middle].append(point)

    for row in result:
        for cell in row:
            yield cell


def distance(point_a, point_b):
    return sum((a - b)**2 for a, b in zip(point_a, point_b))


def filter_close_points(minimal_distance, points):
    seen = set()
    for point in points:
        if (
            point not in seen
            and all(distance(point, p) >= minimal_distance for p in points)
        ):
            seen.add(point)
            yield point


def create_geometry(points):
    geometry = ogr.Geometry(ogr.wkbMultiPoint)
    wkb_point = ogr.Geometry(ogr.wkbPoint)
    for point in points:
        wkb_point.AddPoint(*point)
        geometry.AddGeometry(wkb_point)
    return geometry


class Process(WPSProcess):
    def __init__(self):
        WPSProcess.__init__(
            self,
            identifier=__name__,
            title='Pollmueller Tin-Reducer',
            version = '0.1',
            storeSupported='true',
            statusSupported='true',
            abstract=(
                'reduces the points of a given triangulated irregular network'
                ' by a variable factor'
            ),
            grassLocation=False
        )
        self.data_output = self.addComplexOutput(
            identifier='output',
            title='Output vector data',
            formats=[{'mimeType': 'text/xml'}]
        )
        self.time_output = self.addLiteralOutput(
            identifier='output2',
            title='computation time (to compare different approaches)'
        )
        self.data_filename_input = self.addComplexInput(
            identifier='data',
            title='Input vector data',
            formats=[{'mimeType': 'text/xml'}]
        )
        self.minimal_distance_input = self.addLiteralInput(
            identifier='input2', title='minimal distance between points'
        )

    def execute(self):
        #
        # TODO: Different operating systems have different semantics and
        #       resolutions for `time()` and `clock()`.  Make sure `clock()` is
        #       the correct one here.
        #
        start_timestamp = time.clock()

        shape_file = ogr.Open(self.data_filename_input.getValue())
        if shape_file is None:
            raise Exception('failed to open shape file.')
        #
        # Get the first layer.
        #
        layer = shape_file.GetLayer(0)
        layer.ResetReading()
        #
        # Read points.
        #
        cleaned_points = clean_points(iter_points(layer))
        #
        # Create raster cells with the corresponding points (5% overlap added).
        #
        # TODO: Where is the 5% overlap actually realized?
        #
        cells = partition_points(cleaned_points, layer.GetExtent())
        minimal_distance = self.minimal_distance_input.getValue()
        result_points = chain.from_iterable(
            imap(partial(filter_close_points, minimal_distance), cells)
        )
        gml = create_geometry(result_points).ExportToGML()
        elapsed_time = time.clock() - start_timestamp
        self.data_output.setValue(gml)
        self.time_output.setValue(elapsed_time)