Programmatically infer start and end points from a Polygon with linear characteristics

Question

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:

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.

Answer 1

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.