点和线的扇形缓冲区生成

# coding:utf-8
# 扇形缓冲区
from shapely.ops import substring
from shapely.geometry import LineString, Polygon
from math import pi
import numpy as np
import logging
logging.basicConfig(level=logging.WARNING,
                    format='%(asctime)s-%(filename)s[line:%(lineno)d]-%(levelname)s:%(message)s',
                    datefmt='%H:%M:%S')


def pnt2polar(p):
    # 平面点转极坐标
    x, y = p
    length = np.sqrt(x * x + y * y)
    x = round(x, 6)
    y = round(y, 6)
    alpha = np.arctan(abs(y / x)) if x != 0. else 0.
    if x > 0:
        if y > 0:
            angle = alpha
        elif y == 0:
            angle = 0
        else:
            angle = pi * 2 - alpha
    elif x == 0:
        if y > 0:
            angle = pi / 2
        elif y == 0:
            angle = 0
        else:
            angle = 3 * pi / 2
    else:
        if y > 0:
            angle = pi - alpha
        elif y == 0:
            angle = pi
        else:
            angle = pi + alpha

    return round(length, 3), round(angle, 3)


def polar2pnt(p):
    # 极坐标转平面坐标
    length, angle = p
    eps = 1e-8

    if length == 0:
        x, y = 0
    else:
        if angle < eps:
            x = length
            y = 0
        elif angle - 3 * pi / 2 > eps:
            x = length * np.cos(2 * pi - angle)
            y = -length * np.sin(2 * pi - angle)
        elif angle - 3 * pi / 2 == eps:
            x = 0
            y = length
        elif angle - pi > eps:
            x = -length * np.cos(angle - pi)
            y = -length * np.sin(angle - pi)
        elif angle - pi == eps:
            x = -length
            y = 0
        elif angle - pi / 2 > eps:
            x = -length * np.cos(pi - angle)
            y = length * np.sin(pi - angle)
        elif angle - pi / 2 == eps:
            x = 0
            y = -length
        else:
            x = length * np.cos(angle)
            y = length * np.sin(angle)

    return round(x, 3), round(y, 3)


def arcPointBuffer(p, r, fangle, tangle):
    # p 是目标点
    # r 是缓冲半径
    # fangle, tangle 是扇形角度
    # 0 -> 2*pi x轴正方向，逆时针

    if tangle > fangle:
        angles = np.linspace(fangle, tangle, num=64)
    else:
        angles = np.linspace(fangle, 2*pi, num=32)
        angles_ = np.linspace(0, tangle, num=32)
        angles = np.r_[angles, angles_]
    polars = np.c_[[r] * 64, angles]
    vector = [polar2pnt(_) for _ in polars]

    x, y = np.array(p).tolist()
    vector = [[x+i, y+j] for i, j in vector]

    ply = [[x, y]] + vector + [[x, y]]

    return Polygon(ply)


def arcLineBuffer(geom, r, arc):
    # 在线的终点进行arcBuffer
    # geom
    # r 缓冲半径
    # arc 扇形圆心角

    length = geom.length
    if length < r:
        fl = geom
    else:
        fl = substring(geom, length-r, length)

    ps, pe = fl.coords[0], fl.coords[-1]
    vec = pe[0]-ps[0], pe[1]-ps[1]
    logging.warning(f"vec:{vec}")
    lgth, angle = pnt2polar(vec)

    fangle = angle - arc/2
    tangle = angle + arc/2
    if fangle < 0:
        fangle = 2*pi - (arc/2 - angle)
    if tangle > 2*pi:
        tangle = tangle - pi*2
    logging.warning(f"angle:{angle:.3f}")
    logging.warning(f"ftangle:{fangle:.3f},{tangle:.3f}")
    res = arcPointBuffer(pe, lgth, fangle, tangle)
    return res


if __name__ == '__main__':

    ls = LineString([[0, 0], [10, 0]])
    arc = pi/4
    res = arcLineBuffer(ls, 2, arc)
    res = arcPointBuffer((10, 0), 1, 0, arc)