最近由于项目需求，业主需要圆弧形的期望线图，这样让图片看起来更有活力，且对方向属性的表征也更加明确，但是大概看了一圈常用的建模软件，好像没有提供弧线标注的选项，我也没继续找下去，不如自己造一个吧。正常来说，期望线图层都是两个区域形心点之间的连线，是一条只包括首尾两个点的线段，我们需要基于这两个点去生成一段圆弧。

一开始准备做的时候觉得蛮简单，后面发现角度问题很麻烦，想不清楚的话，结果总是错的。下面记录一下整个过程。

整个过程总结为：基于一对起始坐标点，找到一个半径为$R$的圆，让这两个点落在这个圆上，然后得到这两点之间的圆弧线(短圆弧)线型坐标。

计算圆心位置

首先我们确定$R$，我们不希望圆弧对应的圆心角$\theta$太大，依据三角关系：

$\theta = 2* arcsin(0.5*d/r)$

当然$\theta$是一个锐角，才最符合我们的需求，可以推导出$r/d∈[\sqrt2/2,+\infty]$,这个比例可以依据需求选择，经过多次尝试，我发现$r/d=1.5$，圆弧的效果最舒服。

然后我们利用几何关系计算圆心，这里一看就知道会有两个圆心，一个在路段前进方向的左侧、一个在右侧，我们规定取在前进方向左侧的圆心，利用解析几何的简单知识，可以对圆心进行求解：

def get_cen_loc(o_loc, d_loc, r):
    """
    基于给定的起终点坐标以及半径值计算圆心坐标
    """
    if d_loc[0] == o_loc[0]:
        y0 = y1 = (o_loc[1] + d_loc[1]) / 2
        deltay = (y0 - o_loc[1]) ** 2
        deltax = math.sqrt(r ** 2 - deltay)
        x0 = d_loc[0] - deltax
        x1 = d_loc[0] + deltax
        if np.cross(np.array([d_loc[0] - o_loc[0], d_loc[1] - o_loc[1]]), np.array([o_loc[0] - x0, o_loc[1] - y0])) < 0:
            return [x0, y0]
        else:
            return [x1, y1]
    else:
        C1 = (d_loc[0] ** 2 + d_loc[1] ** 2 - o_loc[0] ** 2 - o_loc[1] ** 2) / 2 / (d_loc[0] - o_loc[0])
        C2 = (d_loc[1] - o_loc[1]) / (d_loc[0] - o_loc[0])
        A = 1 + C2 ** 2
        B = 2 * (o_loc[0] - C1) * C2 - 2 * o_loc[1]
        C = (o_loc[0] - C1) ** 2 + o_loc[1] ** 2 - r ** 2
        y0 = (-B + math.sqrt(B * B - 4 * A * C)) / 2 / A
        y1 = (-B - math.sqrt(B * B - 4 * A * C)) / 2 / A
        x0 = C1 - C2 * y0
        x1 = C1 - C2 * y1
        if np.cross(np.array([d_loc[0] - o_loc[0], d_loc[1] - o_loc[1]]), np.array([o_loc[0] - x0, o_loc[1] - y0])) < 0:
            return [x0, y0]
        else:
            return [x1, y1]

利用参数方程求解起始点所对的弧度

我们的思路是基于起终点的弧度，在两个弧度之间去生成均匀的中间弧度点，然后利用参数方程求解圆弧上的$xy$坐标。本来打算用笛卡尔坐标，但是笛卡尔坐标想要生成均匀的中间点，还要考虑曲率的变化，很麻烦，遂放弃。

然后分别计算两个点($O$和$D$点)在圆中对应的弧度$\alpha，\beta$,利用第一个参数方程来求解，也就是要涉及到$arccos$函数，每次遇到这个函数，都很头疼，涉及到象限和$2\pi$的问题，一不小心就出错。

我用的是numpy的np.arccos函数，他的输出值域是$[0,\pi]$之间，但是我们会做处理，我们在计算OD两个点所对应的弧度$\theta$时，会先计算圆心到这个点的向量，先判断这个向量在第几象限：

如果在1、2象限，直接返回arccos的值，如果在3、4象限，返回$2\pi - arccos$的值。

这样我们就可以计算到起始弧度，正常情况下，起点$O$对应的弧度比终点$D$对应的弧度小，且相差的是锐角(基于我们之前的前提，$r/d$的比例关系)，但是当起点在4象限、终点在1象限时，会有角度的突变：起点的弧度比终点的弧度大，所以我们需要做处理：

得到起终点弧度的对应值后，我们只需要在区间内均匀取样弧度值，然后带入到参数方程中就可以得到弧上的均匀坐标点位。

完整代码

import math
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import LineString


# 主函数
def get_arrow(o_loc=None, d_loc=None, ratio=1.0):
    d = calc_dis(o_loc=o_loc, d_loc=d_loc)
    print(d)
    r = ratio * d
    # 计算圆心, 前进方向右侧
    cen_loc = get_cen_loc(o_loc=o_loc, d_loc=d_loc, r=r)

    print('circle:')
    print(cen_loc)
    start_theta = get_rad_loc(x=o_loc[0], y=o_loc[1], x_r=cen_loc[0], y_r=cen_loc[1], r=r)

    end_theta = get_rad_loc(x=d_loc[0], y=d_loc[1], x_r=cen_loc[0], y_r=cen_loc[1], r=r)

    if 180 * start_theta / np.pi - 180 * end_theta / np.pi >= 180:
        start_theta = - (2 * np.pi - start_theta)

    print(180 * start_theta / np.pi, 180 * end_theta / np.pi)

    # 采样30个点
    delta_theta = np.linspace(start_theta, end_theta, num=30)
    print(delta_theta)
    x_array = [cen_loc[0] + r * np.cos(delta) for delta in delta_theta]
    y_array = [cen_loc[1] + r * np.sin(delta) for delta in delta_theta]

    arc_line = LineString([(x, y) for x, y in zip(x_array, y_array)])

    return arc_line


def get_cen_loc(o_loc, d_loc, r):
    """
    基于给定的起终点坐标以及半径值计算圆心坐标
    """
    if d_loc[0] == o_loc[0]:
        y0 = y1 = (o_loc[1] + d_loc[1]) / 2
        deltay = (y0 - o_loc[1]) ** 2
        deltax = math.sqrt(r ** 2 - deltay)
        x0 = d_loc[0] - deltax
        x1 = d_loc[0] + deltax
        if np.cross(np.array([d_loc[0] - o_loc[0], d_loc[1] - o_loc[1]]), np.array([o_loc[0] - x0, o_loc[1] - y0])) < 0:
            return [x0, y0]
        else:
            return [x1, y1]
    else:
        C1 = (d_loc[0] ** 2 + d_loc[1] ** 2 - o_loc[0] ** 2 - o_loc[1] ** 2) / 2 / (d_loc[0] - o_loc[0])
        C2 = (d_loc[1] - o_loc[1]) / (d_loc[0] - o_loc[0])
        A = 1 + C2 ** 2
        B = 2 * (o_loc[0] - C1) * C2 - 2 * o_loc[1]
        C = (o_loc[0] - C1) ** 2 + o_loc[1] ** 2 - r ** 2
        y0 = (-B + math.sqrt(B * B - 4 * A * C)) / 2 / A
        y1 = (-B - math.sqrt(B * B - 4 * A * C)) / 2 / A
        x0 = C1 - C2 * y0
        x1 = C1 - C2 * y1
        if np.cross(np.array([d_loc[0] - o_loc[0], d_loc[1] - o_loc[1]]), np.array([o_loc[0] - x0, o_loc[1] - y0])) < 0:
            return [x0, y0]
        else:
            return [x1, y1]


def get_rad_loc(x=1.1, y=2.1, x_r=0.1, y_r=0.1, r=1.2):
    vec = np.array([x - x_r, y - y_r])
    # print(vec)
    if vec[0] >= 0 and vec[1] >= 0:
        theta = np.arccos((x - x_r) / r)
        return theta
    elif vec[0] <= 0 and vec[1] >= 0:
        theta = np.arccos((x - x_r) / r)
        return theta

    elif vec[0] <= 0 and vec[1] <= 0:
        theta = np.arccos((x - x_r) / r)
        return 2 * np.pi - theta

    else:
        theta = np.arccos((x - x_r) / r)
        return 2 * np.pi - theta


# 这个距离计算, 可以传入经纬度坐标直接计算
# 如果两个点之间的经度纬度跨越很大, 如果传入的是平面投影坐标, 计算完毕后再转为地理坐标, 弧线会变形
def calc_dis(o_loc=None, d_loc=None):
    d = ((d_loc[0] - o_loc[0]) ** 2 + (d_loc[1] - o_loc[1]) ** 2) ** 0.5
    if d <= 0.00000000001:
        raise  ValueError(r'起终点重合!')
    else:
        return d


if __name__ == '__main__':
    link_gdf = gpd.read_file(r'./line.geojson')
    print(link_gdf)

    link_gdf.plot()
    plt.show()


    link_gdf['coord_list'] = link_gdf['geometry'].apply(lambda x: list((x.coords)))
    link_gdf['o_loc'] = link_gdf['coord_list'].apply(lambda x: x[0])
    link_gdf['d_loc'] = link_gdf['coord_list'].apply(lambda x: x[-1])
    link_gdf['geometry'] = link_gdf[['o_loc', 'd_loc']].apply(lambda x: get_arrow(o_loc=x[0], d_loc=x[1], ratio=1.5), axis=1)

    link_gdf.drop(columns=['coord_list', 'o_loc', 'd_loc'], axis=1, inplace=True)

    link_gdf.plot()
    plt.show()

    link_gdf.to_file(r'./arc_line.geojson', encoding='gbk', driver='GeoJSON')