Order a set of coordinates to create the maximum straight lines

Question

I have a list with geographical coordinates that form a polygon. My goal is to create a polygon with the maximum straight lines as possible, I will then plot the polygon into a map. I have already checked the code described here but it doesn't seem to retrieve the output I'm looking for.

Example:

def order_points(pp):
    lat, lon=zip(*pp)
    coords=np.array(pp)
    centroid=(mean(lat), mean(lon))
    
    # compute centroid
    cent=(sum([p[0] for p in pp])/len(pp),sum([p[1] for p in pp])/len(pp))
    # sort by polar angle
    pp.sort(key=lambda p: math.atan2(p[1]-cent[1],p[0]-cent[0]))
    
    return pp

polygon_coordinates=[(40.5332064016964, -8.78427738454133),
 (41.25358575070879, -8.662562329821816),
 (41.34395074064536, -8.78162348005526),
 (41.34366188895376, -8.662097742704521),
 (41.3427118223502, -8.423053362078234),
 (41.25263867182198, -8.423846579910334),
 (41.339322521192834, -7.945013920157079),
 (40.52870776612322, -7.957832320269227)]

m = folium.Map(location=(41.25387369401857, -8.781923733405996), zoom_start=12)
polygon_coordinates_reorder=order_points(polygon_coordinates)
polygon = folium.Polygon(locations=polygon_coordinates_reorder, color='red', fill_color='red')
m.save("17_05_map.html")
m

The thing is, the output of this is:

But I was looking for something like this (the outline in black):

"polygon with the maximum straight lines as possible" What do you mean by this? — Mushegh, May 21 '23 at 16:05
@Mushegh I meant to say that I want the inner angle between consecutive edges to be aprox 90 degrees. — Beatriz Santos, May 21 '23 at 16:08

score 0 · Answer 1 · answered May 22 '23 at 09:48

your order_points function sorts points by their polar coordinates, NOT by a condition that determines the resulting angle of the lines...

not sure if there's a general (non-ambiguous) way to sort your points so that they form straight lines... (after all there are might many straight lines possible with a set of 8 points...)

if it's just this one polygon, I'd recommend the "manual" way :-) e.g. like this: (done with EOmaps)

from eomaps import Maps
from shapely.geometry import Polygon
import geopandas as gpd
import numpy as np

polygon_coordinates=[(40.5332064016964, -8.78427738454133),
                     (41.25358575070879, -8.662562329821816),
                     (41.34395074064536, -8.78162348005526),
                     (41.34366188895376, -8.662097742704521),
                     (41.3427118223502, -8.423053362078234),
                     (41.25263867182198, -8.423846579910334),
                     (41.339322521192834, -7.945013920157079),
                     (40.52870776612322, -7.957832320269227)]
# reorder coordinates to lon/lat 
polygon_coordinates = np.array(polygon_coordinates)[:,::-1]


m = Maps(Maps.CRS.GOOGLE_MERCATOR)
m.add_wms.GOOGLE.add_layer.roadmap_standard()

# plot data points so we can add nice annotations
m.set_data(None, *polygon_coordinates.T)
m.set_shape.scatter_points()
m.plot_map(fc="k", zorder=2)
m.add_annotation(ID=range(len(polygon_coordinates)), 
                 text=lambda ind, **kwargs: str(ind))

# re-order polygon points and add the resulting polygon to the map
gdf = gpd.GeoDataFrame(
    geometry=[Polygon(polygon_coordinates[[2,3,1,5,4,6,7,0,2]])], 
    crs=4326)
m.add_gdf(gdf, fc="r", ec="k", alpha=0.5)

unfortunately there are several more polygons with different vertices, so manually isn't an option. — Beatriz Santos, May 23 '23 at 14:12

score 0 · Answer 2 · answered Jul 11 '23 at 04:03

Fairly sure this ends up as something of an optimization problem (possible running a restricted version of the traveling salesman problem (TSP) could solve this where the score is the sum of how far the angles are off from [angle % 90, 90 - angle % 90].min rather than the distance.

Start by building the convex hull, call this the candidate polygon
Window through the candidate polygon (including the first and last point) two points at a time. For each point that remains, check if it could be a right angle (or straight line) by adding the point before or after that pair of points. Record the results in a bi-directional hash. IE: hash[(0,7)] = Set(6) and hash[6]= [(0,7)]
Once you have tried for all points in the hull, look for any single vertices with only 1. For example, 0, 7, 1, 5 (though 0 and 7 are part of the hull). Add those to the canidate, removing the paths from the hashes they already add. Repeat step 2, but only for the new points.
At this point you can solve it as a Hamiltonian path problem (given a graph, is it possible to find a path that connects all points). Though this is probably less ideal than the TSP approach, but in theory you may have good odds of being able to solve it quickly and easily.

The gist is trying to build paths from the triples. IE: (2,3,1) -> (1,5,4) -> (4,6,7) -> (7,0,2). Though that only works flawlessly if there is a "perfect" solution that allows for all points to be 180 or 90 degrees. If there is even one that isn't 90/180, then the TSP approach is really needed unless you want to re-run it with a greater allowance. IE: +-10 degrees from 90.

Order a set of coordinates to create the maximum straight lines

2 Answers2