import time 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(nth, iterable): return (item for i, item in enumerate(iterable, 1) if i % nth != 0) def distance(point_a, point_b): # # TODO: To be the usual euklidean distance formula it lacks a `math.sqrt()` # call! # return sum((a - b)**2 for a, b in zip(point_a, point_b)) def filter_close_points(minimal_distance, points): # # TODO: The reversing is only neccessary if the points have to be looked # at in the order they are given to this function. If going through # them in reverse would be okay we could spare this step. # points = list(reversed(points)) while points: point = points.pop() point_count = len(points) points = [p for p in points if distance(point, p) >= minimal_distance] if len(points) == point_count: 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() points = filter_close_points( self.minimal_distance_input.getValue(), clean_points(iter_points(layer)) ) gml = create_geometry(points).ExportToGML() elapsed_time = time.clock() - start_timestamp self.data_output.setValue(gml) self.time_output.setValue(elapsed_time)