Search code examples
pythongislatitude-longitudegeopandasshapely

Overlapping Polygon in Latitude/Longitude Space

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. enter image description here

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

Solution

  • 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()

    enter image description here