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