计算线段包围一点形成的包围角度

---

前言

通常山脉对降雨的影响比较大，比如秦岭山脉阻挡了北上的水汽，使南北方呈现出差异较大的降水水平。除此之外，其他小的山脉也会对降水有一定的影响，因此，研究降水的影响因子时，山脉是个不可忽略的因素。求山脉所围成的包围角和阻挡方位是量化山脉影响的方法之一，其中，山脉可以用线表示（如等高线），气象站点用点表示。这里提供一种计算山脉所围成包围角度和阻挡方位的算法

---

一、算法原理

对上述实例进行数学化描述：计算线段包围一点形成的包围角度和该线段对该点的阻挡方位。

通常较为简单的一种场景是：

上述线条AB对点O形成的包围角为AOB，只要指导线条起点和终点可以很容易求出角度。

由于显示场景较为复杂，存在多种多样的情形，通过起点和终点方式则不能轻易计算得到，如下面一些情形。


算法完备性：寻找共同规律。

线是由一串点相连构成的，具体为一条线由n个点构成（a, b, c, d, …），由n-1条小线段构成（ab, bc, cd, …）。

我们定义每条小线段与o构成的夹角，即aob,boc,cod,…，可以通过点乘计算；进一步定义角度的正负，叉乘向里为负，向外为正。并对所有小线段角度累加，以下是不同案例的统计情况：

可以发现，对角度累积求和，极大值和极小值所在点即是A和B；包围角度为极大值减去极小值。

二、代码实现（python）

NORTH_VEC = (0, 100) # 定义正北向量
EAST_VEC = (100, 0)  # 定义正东向量

def distance(x1, y1, x2, y2):   # (x1, y1),(x2, y2) 两点距离
    return math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))

def vector(point1, point2):   #  两点形成的向量
    return (point2[0] - point1[0], point2[1] - point1[1])

def angle(vec1, vec2):  #  向量夹角--带方向
    vec1dotvec2 = vec1[0] * vec2[0] + vec1[1] * vec2[1]
    absVec1 = math.sqrt(vec1[0] * vec1[0] + vec1[1] * vec1[1])
    absVec2 = math.sqrt(vec2[0] * vec2[0] + vec2[1] * vec2[1])
    ang = math.acos(vec1dotvec2 / (absVec1 * absVec2))
    if vec1[0] * vec2[1] - vec1[1] * vec2[0] < 0:
        ang = -ang
    return ang

def myMethod(point, line):
	"""
	point:点
	line:线
	"""
    middlePoint = line[int(len(line)/2)]
    endPoint = line[len(line) - 1]

    vecEndPoint = vector(point, endPoint)
    
    sumAngleList = []
    sumAngle = 0
    DirecList = []
    sumAngleList.append(sumAngle)
    for i in range(len(line) - 1):
        point1 = line[i]
        point2 = line[i + 1]
        vec1 = vector(point, point1)
        vec2 = vector(point, point2)
        ang = angle(vec1, vec2)

        Direc = math.fabs(angle(vec1, NORTH_VEC))   #  与正北方向形成的夹角，有两种情况，在2/3象限或在1/4象限。如果在2/3象限需要再运算：2*pi-夹角
        if math.fabs(angle(vec1, EAST_VEC)) > math.pi / 2:  #  如果同时与正东方向夹角大于pi/2，可以确定在2/3象限
            Direc = 2 * math.pi - Direc
        DirecList.append(Direc)

        sumAngle += ang
        sumAngleList.append(sumAngle)

    Direc = math.fabs(angle(vecEndPoint, NORTH_VEC))
    if math.fabs(angle(vecEndPoint, EAST_VEC)) > math.pi / 2:
        Direc = 2 * math.pi - Direc
    DirecList.append(Direc)

    sumAngleMax = -7
    sumAngleMin = 7
    startDirec = 0
    endDirec = 0
    for k in range(len(sumAngleList)):  #  找累加最大值和最小值,最大值-最小值为包围角度，对应的方向为起始包围方向和终止包围方向。
        if sumAngleList[k] > sumAngleMax:
            startDirec = DirecList[k]
            sumAngleMax = sumAngleList[k]
        if sumAngleList[k] < sumAngleMin:
            endDirec = DirecList[k]
            sumAngleMin = sumAngleList[k]
    bwAngle = sumAngleMax-sumAngleMin