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, x=coord) # add edges to the network for a, b, i in zip(coords[:-1], coords[1:], ids): d = great_circle_vec(a, a, b, b) a_id = i b_id = i + 1 if b_id > ids[-1]: b_id = ids 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)
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.