xmet/lib/xmet/geo.py

import re
import enum
import math
import shapely

from typing import Self

class PointDirection(enum.Enum):
    EQUAL = 0
    LEFT  = enum.auto()
    RIGHT = enum.auto()

class PointSequence(list):
    linestring: shapely.LineString
    polygon:    shapely.Polygon

    @staticmethod
    def from_file(path: str) -> Self:
        ret = PointSequence()

        with open(path, 'r') as fh:
            data = fh.read()

            for line in data.split('\n'):
                if line == '':
                    continue

                lat, lon = re.split(r'\s*,\s*', line)

                ret.add(float(lon), float(lat))

        return ret

    def __init__(self, points: list=None):
        super().__init__()

        self.linestring = None
        self.polygon    = None

        if points is not None:
            for point in points:
                typeof = type(point)

                if typeof is tuple:
                    self.add(*point)
                elif typeof is shapely.Point:
                    self.append(point)

    def add(self, lon: float, lat: float):
        self.append(shapely.Point(lon, lat))

    def is_closed(self) -> bool:
        return self[-1] == self[0]

    def close(self):
        if not self.is_closed():
            self.append(self[0])

        self.linestring = shapely.LineString(self)
        self.polygon    = shapely.Polygon(self)

    def nearest_index(self, point: shapely.Point) -> int:
        indices = list()

        for i in range(0, len(self)):
            indices.append((i, self[i].distance(point)))

        indices.sort(key=lambda i: i[1])

        return indices[0][0]

    def index_distance(self, i1: int, i2: int) -> int:
        """
        Returns the index distance of i1 relative to i2, and whether i1 is
        considered left of, equal to, or right of i2.
        """
        count = len(self)
        value = count - ((i1 - i2) % count)

        if value == 0:
            direction = PointDirection.EQUAL
        elif value > count / 2:
            direction = PointDirection.RIGHT
        else:
            direction = PointDirection.LEFT

        return value, direction

class PolygonBuilder():
    point_first: shapely.Point
    point_last:  shapely.Point

    def __init__(self, sequence: PointSequence, bounds: PointSequence):
        self.sequence    = sequence
        self.bounds      = bounds
        self.point_first = None
        self.point_last  = None
        self.total       = 0

    def yield_point(self, point: shapely.Point):
        """
        Yield the single point to the caller, while maintaining state of
        whether this point is the first or last point seen by the polygon
        builder.
        """
        if self.point_first is None:
            self.point_first = point

        self.point_last = point

        yield point

        self.total += 1

    def each_intermediate_point(self, p1: shapely.Point, p2: shapely.Point):
        count = len(self.bounds)

        i1 = self.bounds.nearest_index(p1)
        i2 = self.bounds.nearest_index(p2)

        dist, direction = self.bounds.index_distance(i1, i2)

        if direction is not PointDirection.LEFT:
            return

        for i in range(i1, i2+1):
            yield from self.yield_point(self.bounds[i % count])

    def each_point_within(self):
        last = None

        for point in self.sequence:
            if last is None:
                last = point
                continue

            last_within  = self.bounds.polygon.contains(last)
            point_within = self.bounds.polygon.contains(point)

            #
            # If the first point in the current line exists within the other
            # geometry, then yield it.
            #
            if last_within:
                yield from self.yield_point(last)

            #
            # Check for intersections with the line in the other geometry.
            #
            inter = self.bounds.linestring.intersection(shapely.LineString([last, point]))

            if inter.geom_type == 'Point':
                #
                # If the intersection is a single point, then yield that
                # point.
                #
                yield from self.yield_point(inter)
            elif inter.geom_type == 'MultiPoint':
                #
                # If the intersection is multiple points, then yield those
                # points, as well as all between on the other geometry, if and
                # only if the intersection does not constitute the first and
                # last point.
                #
                last_geom = None

                for geom in inter.geoms:
                    if last_geom is None:
                        last_geom = geom
                        continue

                    yield from self.yield_point(last_geom)
                    yield from self.each_intermediate_point(last_geom, geom)
                    yield from self.yield_point(geom)

                    last_geom = geom

            #
            # If the second point in the current line exists within the other
            # geometry, then yield that.
            #
            if point_within:
                yield from self.yield_point(point)

            #
            # Yield all intermediate points if the first point is to the right
            # of the last point.
            #
            if self.point_first is not None and self.point_last is not None:
                yield from self.each_intermediate_point(self.point_last,
                                                        self.point_first)

            last = point

    def process(self) -> shapely.Polygon:
        pass

def heading(p1: shapely.Point, p2: shapely.Point) -> float:
    dx = p2.x - p1.x
    dy = p2.y - p1.y

    return math.atan2(dy, dx)