基于3点构建圆/弧

# coding:utf-8
import numpy as np
import pandas as pd
import geopandas as gpd
from scipy import interpolate
from functools import partial
from shapely.geometry import Point, LineString
from shapely.affinity import rotate, scale
from shapely.ops import split, linemerge, snap
import logging
logging.basicConfig(level=logging.WARNING,
                    format='%(asctime)s-%(filename)s[line:%(lineno)d]-%(levelname)s:%(message)s',
                    datefmt='%Y-%m-%d %H:%M:%S')


class arcBase3pnts():
    # 基于3点构建圆/弧
    
    def __init__(self, pnts):
        p1, p2, p3 = pnts
        self.p1 = p1
        self.p2 = p2
        self.p3 = p3
    
    def getCenterAndRadius(self):
        # pnts = [[x1,y1],[x2,y2],[x3,y3]]
        p1, p2, p3 = self.p1, self.p2, self.p3
        
        ls1, ls2 = LineString([p2, p3]), LineString([p1, p3])
        fx = partial(scale, xfact=100, yfact=100)
        ls1, ls2 = fx(ls1), fx(ls2)
        vs1, vs2 = rotate(ls1, 90), rotate(ls2, 90)
        
        cp = vs1.intersection(vs2)
        r = cp.distance(Point(p1))
        return cp, r 
            
    
    def getCircle(self):
        # pnts = [[x1,y1],[x2,y2],[x3,y3]]
        cp, r = self.getCenterAndRadius()
        
        return cp.buffer(r)
    
    def getArc(self):
        # 获取过三点的弧长
        # 要求弧长过点的顺序按照点的顺序，即中间点在弧长的中间部分
        p1, p2, p3 = self.p1, self.p2, self.p3
        circle = self.getCircle()
        ls = LineString([p1, p3])
        parts = split(circle, ls)
        logging.warning(f"[{p1},{p2},{p3}]")
        logging.warning(f"nparts:{len(parts)}")
        
        lg = [part.distance(Point(p2)) for part in parts]
        idx = np.argmin(lg)
        part = parts[idx]
        
        res = circle.boundary.intersection(part)
        tmp = []
        for i in range(len(res)):  # 处理可能出行的多线情况
            gtype = res[i].geom_type
            # logging.warning(f"gtype:{gtype}")
            if gtype == "LineString":
                tmp.append(res[i])
        res2 = linemerge(tmp)        
        return res2
    
    def getLastArc(self):
        # 获取过后两个点之间的弧
        p2 = self.p2
        arc = self.getArc()
        pp2 = Point(p2)
        toler = np.ceil(arc.distance(pp2))
        arc = snap(arc, pp2, tolerance=toler)
        parcs = split(arc, pp2)
        index = 0 if parcs[0].distance(pp2)<parcs[0].distance(pp2) else 1
        return parcs[index]


if __name__ == "__main__":
    
    pnts = [(13490633.7293, 4161607.44897),(13497014.3574, 4149089.65039),(13499473.5673, 4141429.48048)]
    oarc = arcBase3pnts(pnts=pnts)
    arc = oarc.getArc()
    prc = oarc.getLastArc()