I have a set of longitudes and latitudes for which I need to create an enclosing Polygon. In some cases, these sets of points overlap and creating a Polygon has been giving me a lot of trouble. I have attempted to use geopandas and shapely to determine the Polygons, find the overlaps, etc but with no luck.
The following image shows what I mean. In it the highlighted areas are what it considers to be the Polygons but it is adding to rings of latitude that I do no have in my points and instead what I need is the the ring that goes behind the Earth in this image and then the two un-highlighted sections (one with a concave end and the other with a convex end). I think this can be accomplished with POLYGON((), (line2)) where this would represent a POLYGON with a hole but I can't seem to create it.
The points are in a csv file, https://drive.google.com/file/d/1Hcj5THoc2jNcY-AEkrmB1bLcCIhCd_FX/view?usp=drive_link, where the the rows are:
longitude, latitude
Some of what I have tried is as follows:
from shapely.geometry import Polygon, Point
import geopandas
p = geopandas.GeoSeries([Point(x,y) for x,y in paired_coords])
p = Polygon(p)
s = geopandas.GeoSeries([Polygon(paired_coords)])
union_all = s.union_all()
gunion_all = geopandas.GeoSeries(union_all)
r = s.remove_repeated_points()
v = r.make_valid()
exploded = v.explode(index_parts=True)
u = s.union(s)
I am 100% sure there will be questions ,I described it the best I could think to but expect there to be something not described well enough
Shapely is not natively suitable for spherical coordiantes, to get around that I will split the input when passing the 180° longitude line.
After that, the geoms 'behave' as one would expect, and a simple intersection between western and eastern part nets the middle section.
from itertools import pairwise
import numpy as np
import pandas as pd
from shapely import Polygon, LineString, LinearRing
def split_ring_at_jump(ring: LinearRing, threshold=90) -> list[LineString]:
"""
Split a LinearRing at points where the step size between consecutive coordinates exceeds a threshold.
"""
coords = np.array(ring.coords)
if len(coords) < 2:
return [ring]
diff_x = np.abs(np.diff(coords[:, 0]))
diff_y = np.abs(np.diff(coords[:, 1]))
step_size = np.sqrt(diff_x ** 2 + diff_y ** 2)
split_indices = np.where(step_size > threshold)[0] + 1
if len(split_indices) == 0:
return [ring]
segments = []
for start, end in pairwise(split_indices.tolist()):
if start + 1 == end:
# single points
continue
segment_coords = coords[start:end]
segments.append(LineString(segment_coords))
start_end = LineString(np.vstack([coords[split_indices[-1]:], coords[0: split_indices[0]]]))
segments.append(start_end)
return segments
df = pd.read_csv(r"C:\Users\victo\Downloads\points.csv", header=None)
ring = LinearRing(df)
splits = split_ring_at_jump(ring)
splits_polygon = [Polygon(geom) for geom in splits]
intersection_polygon = shapely.intersection(Polygon(splits[0]), Polygon(splits[-1]))
polygons_cut = [poly - intersection_polygon for poly in splits_polygon]
For visualization I use shapely.plotting
:
from matplotlib import pyplot as plt
import matplotlib.colors as mcolors
import shapely.plotting
for c, poly in zip(list(mcolors.BASE_COLORS), polygons_cut):
shapely.plotting.plot_polygon(poly, color=c)
shapely.plotting.plot_polygon(intersection_polygon, color='red')
plt.show()