transform_points_to_unit_circle.py

import math
from itertools import pairwise

import geopandas as gpd
import numpy as np
from shapely import LineString, Point, Polygon, distance
from shapely.coords import CoordinateSequence


def locate_farthest_intersection_point(
    center_point: Point,
    convex_hull_boundary: LineString,
    angle: float,
    raise_if_multiple: bool = False,
) -> Point | None:
    ray_length = 1e5
    angle_rad = math.radians(angle)
    ray_endpoint = Point(
        center_point.x + ray_length * math.cos(angle_rad),
        center_point.y + ray_length * math.sin(angle_rad),
    )
    ray = LineString([center_point, ray_endpoint])
    intersection = convex_hull_boundary.intersection(ray)

    if intersection.is_empty:
        return None
    elif intersection.geom_type == "Point":
        return intersection
    elif intersection.geom_type in ["MultiPoint", "GeometryCollection"]:
        if raise_if_multiple:
            raise RuntimeError
        points = [geom for geom in intersection.geoms if geom.geom_type == "Point"]
        if not points:
            return None
        closest_point = max(points, key=lambda point: center_point.distance(point))
        return closest_point
    else:
        return None


def get_angle(point1: Point, point2: Point) -> float:
    rads = np.arctan2(point2.y - point1.y, point2.x - point1.x)
    return np.degrees(rads)


def transform_point(
    point: Point, center_point: Point, isochrone_boundary: Polygon
) -> Point:
    angle = get_angle(center_point, point)
    intersection_point = locate_farthest_intersection_point(
        center_point, isochrone_boundary.exterior, angle
    )

    distance_from_isochrone_boundary = distance(center_point, intersection_point)
    distance_from_current_point = distance(center_point, point)
    distance_ratio = min(
        1,
        distance_from_current_point / distance_from_isochrone_boundary
        if distance_from_isochrone_boundary > 0
        else 0,
    )

    length = distance_ratio

    angle_rad = np.radians(angle)

    new_point = Point(length * np.cos(angle_rad), length * np.sin(angle_rad))

    return new_point


def transform_coords(
    coords: CoordinateSequence,
    center_point: Point,
    isochrone_boundary: Polygon,
) -> list[Point]:
    result = []
    coords_list = list(coords)
    for (x1, y1), (x2, y2) in pairwise(coords_list):
        result.append(transform_point(Point(x1, y1), center_point, isochrone_boundary))
        ls = LineString([(x1, y1), (x2, y2)])
        step = 1 / 8
        for idx in np.arange(0, 1, step):
            result.append(
                transform_point(
                    ls.interpolate(idx, normalized=True),
                    center_point,
                    isochrone_boundary,
                )
            )

    result.append(
        transform_point(Point(coords_list[-1]), center_point, isochrone_boundary)
    )
    return result


def transform_geometries(
    gs: gpd.GeoSeries,
    center_point: Point,
    isochrone_boundary: Polygon,
) -> gpd.GeoSeries:
    geoms = []
    for geometry in gs:
        if isinstance(geometry, Polygon):
            transformed_ex = transform_coords(
                geometry.exterior.coords,
                center_point,
                isochrone_boundary,
            )
            transformed_ins = [
                transform_coords(interior.coords, center_point, isochrone_boundary)
                for interior in geometry.interiors
            ]
            geoms.append(Polygon(transformed_ex, transformed_ins))
        elif isinstance(geometry, LineString):
            transformed_coords = transform_coords(
                geometry.coords, center_point, isochrone_boundary
            )
            geoms.append(LineString(transformed_coords))

    return gpd.GeoSeries(geoms)


# Usage: transform_geometries(clipped_buildings)