Computational geometry + geo-computations: polyline simplification, convex hulls

DP_reg.py

#!/usr/bin/python
"""
Created on Mon Jan 17 2016
Anna M. Kedzierska

Python3 implementation of the Douglas-Peucker algorithm for polyline simplification
"""

from functools import singledispatch
from functools import reduce
import numpy as np
from numpy import dot
import math

import json
import geojson

# define a Point
class PointLoc:
    def __init__(self, p):
        self.x, self.y  = p[0], p[1]
# define a Line (LineString)
class Line:
    def __init__(self, start, end):
        self.start, self.end = start, end

# class with implementations of various distances: point-point, point-segment
class CalculateDistance(object):

    @singledispatch
    def get_dist(shape, point):
        pass

    @get_dist.register(PointLoc)
    def _(arg, point):
        """distance between two points"""
        dx = arg.x - point.x
        dy = arg.y - point.y

        return math.hypot(dx, dy)

    @get_dist.register(Line)
    def _(arg, point):
        """distance between a point and a segment"""
        point = PointLoc(point)
        dvx = arg.end[0] - arg.start[0]
        dvy = arg.end[1] - arg.start[1]
        #distance to the first point
        dwx = point.x - arg.start[0]
        dwy = point.y - arg.start[1]

        v = [dvx, dvy]
        w = [dwx, dwy]
         # c1/2 is an area of the triangle
        c1 = dot_product(w,v)

        # check if the 1st end of the segment is the closest of all (the sign of the dot products checks the side of the points)
        if c1 <= 0:
          return math.hypot(dwx,dwy)

        c2 = dot_product(v,v)
        #distance of the point to the second end of the segment, which is closer then the 1st
        if c2 <= c1:
          dzx = point.x - arg.end[0]
          dzy = point.y - arg.end[1]
          return math.hypot(dzx, dzy)

        # none of the endpoints is closest to he point, so calculate the projection onto the segment
        b = c1 / c2
        proj1 = arg.start[0] + b * v[0]
        proj2 = arg.start[1] + b * v[1]
        out = math.hypot(point.x - proj1, point.y - proj2)

        return out

def check(v1,v2):
    if len(v1)!=len(v2):
        raise ValueError("the length of both vectors/arrays must be the same")
    pass

def cross_product(v1, v2):
    check(v1, v2)
    return np.cross(v1, v2)

def dot_product(v1, v2):
    check(v1, v2)
    return np.dot(v1, v2)

# impementation of Radial Distance alg= remove points 'too close' to each other
# on the line
def estimate_radial_distance(inLine, tolerance):

  def radial_dist_myReduce(resultsList, elem):
     if len(resultsList)==0:
           resultsList = [elem]

     lastElem = resultsList[-1]

     if CalculateDistance.get_dist(PointLoc(lastElem), PointLoc(elem)) > tolerance:

        return resultsList + [elem]

     else:

        return resultsList

  return reduce(radial_dist_myReduce, inLine, [])

#recursive implementation of the DP algorithm for polyline simplification
def get_dpeucker_rec(points, threshold):
    """
    Recursive version of the DP algorithm
    """
    xmax = 0.0 # we want floats
    xindex = 0

    L =len(points)

    TabLines = np.empty(L,dtype=Line)
    TabLines.fill(Line(points[0], points[-1]))
    # alternative recommended for python3: list comprehensions
    x = list(map(CalculateDistance.get_dist, TabLines, points))
    xmax = max(x)

    xindex = np.argmax(x)

    if xmax >= threshold:
        # [:-1] index is set no to add the same element, which appears in the set from the same split
        results = get_dpeucker_rec(points[:xindex+1], threshold)[:-1] + get_dpeucker_rec(points[xindex:], threshold) # adding sets (left and right split)
    else:
        results = [points[0], points[-1]] # collect the segments with no points above the threshold (no more splits)

    return results

def simplify_lineDP(points, tolerance, RadialPass):
    # if Radial Distance == True, it is run before DP
    if RadialPass:
        points = estimate_radial_distance(points, tolerance)
    points = get_dpeucker_rec(points, tolerance)
    return points

# printing simplified points as a geojson
def create_map(points, type_obj, out_file_name):
    """
    Input: list of points as tuples, type_obj= 0: points, else: line
    Returns a GeoJSON file.
    Note: GitHub can automatically render the GeoJSON
    file as a map.
    """

    # Define the GeoJSON type
    geo_map = {"type": "FeatureCollection"}

    # Define empty list to collect each point
    item_list = []
    index = 0
    # Iterate over our data to create GeoJSOn document.
    for point in points:
        index = index + 1
        # beware of 0's
        if point[0] == "0" or point[1] == "0":
            continue

        data = {}
        data['type'] = 'Feature'
        data['id'] = 'Stop: '+ str(index)
        data['geometry'] = {'type': 'Point',
                            'coordinates': (point[0], point[1])}

        item_list.append(data)

    for point in item_list:
        geo_map.setdefault('features', []).append(point)

    with open(out_file_name + '.geojson', 'w') as f:
        f.write(geojson.dumps(geo_map))
        
# the following area calculation will be used for alg improvement + Visvalingam Whyatt implementation        
def get_triangle_area(p0, p1, p2):
    # Note to self: the area of a planar parallelogram or triangle can be expressed
    # by the magnitude of the cross-product of two edge vectors: 0.5|v x w|, v = w,
    vx = p1.x - p0.x
    vy = p1.y - p0.y

    wx = p2.x - p0.x
    wy = p2.y - p0.y


    v = [vx, vy]
    w = [wx, wy]
    area = 0.5 * np.cross(v,w)
    # explicit formulae: vx*wy-wx*vy
    # Important: the sign gives additional info about the placement of the point wrt to a segment/line,
    # which corresponds to the directionality of the polygon created by the points or the associated convex hull
    # e.g. >0  counterclockwise direction (point is to the left of p0p1); <0 is cw, point to the right of p0p1
    return area        
#########

if __name__ == "__main__":

    # read in or process the input data in a geojson format
    data = []
    # the json core:
    with open('route_input.geojson') as json_file:
        json_data = json.load(json_file)
    #for geojson specifically:
    for feature in json_data['features']:
        data.append(feature['geometry']['coordinates'])


    # simplify the polyline using the DP algorithm.
    # call with 'True' if Radial Distance should be a pre-step
    tolerance = 0.003
    data_simplified = simplify_lineDP(data, tolerance, False)
    # print to a geojson under pre-specified name (e.g. simplified_route_out.geojson)
    create_map(data_simplified, 1, 'simplified_route_out')

   
# TODO:
# - Improve implementation to avoid intersections.
#   Step1: compute a convex hull
#   The problem with traditional formulation of the DP algorithm is that
#   the simplified polyline might not be homeomorphic to the oryginal one
#   if self-intersections occur (no topological consistency).
#   Improved algorithm design has at its basis the computation of convex hulls of the polyline(s).
#   1) speeding up the DP algorithm implemented here
#       (https://www.cs.swarthmore.edu/~adanner/cs97/s08/pdf/speedingPeuckerDouglas.pdf) and
#   2) adding conditions that avoid creating self-intersections
#   (http://www.scielo.br/pdf/jbcos/v9n3/06.pdf)
# Best Choice: Melkman's algorithm
# - Visvalingam-Wyatt simplification

process_bus_data.py

#!/usr/bin/python
"""
Created on Mon Jan 17 2016
Anna M. Kedzierska

Process metro data from the LA transport API:
http://developer.metro.net/introduction/realtime-api-overview/realtime-api-returning-json/
used to create: route_input.geojson
"""
import requests
import json
import geojson

def process_bus_data():
    """
    Extract and parse the data
    """
    res_request = requests.get("http://api.metro.net/agencies/%20lametro/routes/704/sequence/")
    res_request.raise_for_status() # check the response object for errors
    data_req = res_request.json()
    i = 0
    loc_data =[]
    while i<len(data_req['items']):
        stops_lon = data_req['items'][i]['longitude']
        stops_lat = data_req['items'][i]['latitude']
        loc_data.append((stops_lon, stops_lat))
        i+=1
    return loc_data

route_input.geojson

simplified_route_out.geojson