4

Question: Given a Polygon exhibiting linear-like qualities (long and thin, with clear extremities upon visual inspection), is there a reasonable/defensible order of operations to compute "start" and "end" points from the geometry?

A notebook including a code snippet (somewhat involved, but minimal to generate an example case) is included here and pairs with the current strategy I have proposed below: Gist link. The notebook is more readable but the contents are also included below:

import geopandas as gpd
import numpy as np
import networkx as nx

from shapely.geometry import LineString, MultiPolygon, Point, Polygon
from shapely.ops import triangulate

coords = [ [ -24.43359375, -21.943045533438166 ], [ -24.43359375, -22.39071391683855 ], [ -23.90625, -22.87744046489713 ], [ -22.631835937499996, -22.471954507739213 ], [ -21.665039062499996, -20.797201434306984 ], [ -21.533203125, -18.562947442888298 ], [ -20.830078125, -15.368949896534705 ], [ -19.335937499999996, -12.340001834116316 ], [ -16.34765625, -11.092165893501988 ], [ -13.359375, -10.962764256386809 ], [ -10.810546875, -10.660607953624762 ], [ -7.55859375, -7.231698708367139 ], [ -6.328125, -5.047170736919708 ], [ -6.240234374999999, -2.67968661580376 ], [ -6.240234374999999, -0.3076157096439005 ], [ -7.734374999999999, -0.21972602392080884 ], [ -8.173828125, -2.7235830833483856 ], [ -8.61328125, -5.572249801113899 ], [ -10.107421874999998, -7.231698708367139 ], [ -12.65625, -8.537565350804018 ], [ -17.138671875, -9.188870084473393 ], [ -21.005859375, -11.781325296112277 ], [ -22.060546874999996, -14.902321826141796 ], [ -22.4560546875, -18.771115062337 ], [ -24.0380859375, -21.16648385820657 ], [ -24.43359375, -21.943045533438166 ] ]

# note that the base shape has a clear, linear shape
base_shape = Polygon(coords)

# convert the shape into its composite triangles
tri_cleaned = []
for poly in MultiPolygon(triangulate(base_shape, tolerance=0.0)):
    # and only keep those that are inside of the parent shape
    if poly.centroid.intersects(base_shape):
        tri_cleaned.append(poly)

def great_circle_vec(lat1, lng1, lat2, lng2, earth_radius=6371009):
    phi1 = np.deg2rad(90 - lat1)
    phi2 = np.deg2rad(90 - lat2)

    theta1 = np.deg2rad(lng1)
    theta2 = np.deg2rad(lng2)

    cos = (np.sin(phi1) * np.sin(phi2) * np.cos(theta1 - theta2) + np.cos(phi1) * np.cos(phi2))

    arc = np.arccos(cos)

    # return distance in units of earth_radius
    distance = arc * earth_radius
    return distance

# convert the triangles into a network
G = nx.MultiDiGraph()

node_count = 0
for tri in tri_cleaned:
    coords = list(tri.exterior.coords)
    ids = list(map(lambda x: node_count + x, [1, 2, 3]))

    # add nodes to the network
    for one_id, coord in zip(ids, coords[0:3]):
        G.add_node(one_id, y=coord[0], x=coord[1])

    # add edges to the network
    for a, b, i in zip(coords[:-1], coords[1:], ids):
        d = great_circle_vec(a[1], a[0], b[1], b[0])

        a_id = i
        b_id = i + 1

        if b_id > ids[-1]:
            b_id = ids[0]

        G.add_edge(a_id, b_id, length = d)
        G.add_edge(b_id, a_id, length = d)

    node_count += 4

# add links between the triangles' nodes that are touching
nodes_as_pts = []
node_ids = []

for node_id, xy in G.nodes(data=True):
    nodes_as_pts.append(Point(xy['x'], xy['y']))
    node_ids.append(node_id)

node_gdf = gpd.GeoDataFrame(node_ids, geometry=nodes_as_pts, columns=['node_id'])

for node_id, xy in G.nodes(data=True):
    # get all nodes that are adjacent to the node under review
    node_pt = Point(xy['x'], xy['y'])
    ds = node_gdf.distance(node_pt)
    sub = node_gdf[ds < 0.0001]

    # create the connectors
    for to_id in sub['node_id'].values:
        G.add_edge(node_id, to_id, length = 0)
        G.add_edge(to_id, node_id, length = 0)

# at this point you could potentially calculate the shortest path 
# between all nodes to all nodes and, from that, identify the two nodes
# that are most far apart but generating an entire distance matrix is
# extremely slow/expensive... there must be a better way
nx.shortest_path(G)

Example: enter image description here

In the above image I take a cloud of points from a GPS trace and buffer each. I then union them to create a Polygon shape of the point cloud coverage. This gives me such a "long and skinny" Polygon. This is the type of Polygon I would like to acquire "start" and "end" points from.

Current strategy: Simplify the shape and then triangulate it. From triangulation, convert Polygons into a network graph. With graph, find the two points that are farthest apart in the network.

Shortcomings of this method: This is not easy with the existing API (network is a NetworkX graph) and very slow.

Context: As an aside, I've written more about this process here, if further context is desired.

6
  • 1
    I like your existing approach. An alternative to a triangulation + network graph would be to compute the straight skeleton and then find the pair of leaf nodes with the largest distance. The network of the straight skeleton would be orders of magnitude simpler than one based on triangulation. However it also just feels as though you're missing a temporal dimension with your GPS data that would give you a trip start/end. Commented Aug 25, 2017 at 1:05
  • Thanks, the reason I can't use the temporal attributes is that these point clouds are created by multiple trips of different segments being covered at different times (and some trips potentially even stating in the middle, heading to an extremity, and then finishing in the middle, as mentioned in linked blog). I was interested in using straight skeleton but was hoping to find an implementation of it that could be executed outside of PostGIS. If you know of one I would be very interested.
    – kuanb
    Commented Aug 25, 2017 at 1:19
  • 1
  • Minimum spanning tree of original points and longest path?
    – FelixIP
    Commented Aug 25, 2017 at 3:29
  • 1
    Yes, it is exactly what I am suggesting, this will save time on very heavy operations like buffering and dissolving. At the end of the day this is what you are doing, but hopefully faster
    – FelixIP
    Commented Aug 26, 2017 at 1:49

1 Answer 1

4

You may ned to investigate what is often called a medial axis transform and sometimes a skeletonization. From wiki/Medial_axis:

The medial axis of a simple polygon is a tree whose leaves are the vertices of the polygon, and whose edges are either straight segments or arcs of parabolas. The medial axis together with the associated radius function of the maximally inscribed discs is called the medial axis transform (MAT).

Or this from Umea University, Computing Science may be more helpful:

Input description: A polygon or polyhedron P.

Problem description: What is the set of points within P that have more than one closest point on the boundary of P?

Discussion: The medial-axis transformation is useful in thinning a polygon, or, as is sometimes said, finding its skeleton. The goal is to extract a simple, robust representation of the shape of the polygon.

Stefan Huber's PhD thesis states

The straight skeleton is a geometric structure that is similar to generalized Voronoi diagrams.

...

The primary goal of this thesis is the development of an algorithm that is suitable for implementation and efficient in terms of time and space in order to make the advantages of straight skeletons available to real-world applications.

...

We present an algorithm that constructs planar straight-line graphs whose straight skeleton approximates any given motorcycle graph to an arbitrary precision.

And Yongha Hwang has a Python implementation that generates straight skeletons. You probably want the longest path in that skeleton.

3
  • Thanks much Martin. FelixIP pointed, in the comments under the main question, to such an implementation: bit.ly/2x2YhC9 . Once skeletonized, the next step would be (I assume) to find the longest network in that shape (observed as a network graph). If you'd like to flesh out in pseudocode or prose to confirm that, I'd be happy to accept this answer.
    – kuanb
    Commented Aug 30, 2017 at 1:36
  • @kuanb - I'm not qualified to go into details you ask for but have added notes from a reference that Richard provided and an implementation you say Felix mentions.
    – Martin F
    Commented Aug 31, 2017 at 18:16
  • Thanks - didn't mean to add any excessive work to your answer. I think that's great, just wanted to make sure all the resources from this question were in you answer. :)
    – kuanb
    Commented Aug 31, 2017 at 19:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.