如何无需用户提供内部点，用Python判断点是否在球面多边形内？

嗨，这个问题确实戳中了球面几何和平面几何的核心差异——球面多边形的“内部”和“外部”都是有限区域，常规平面的点-in-多边形方法直接套用会有歧义，但我们完全可以通过定义多边形的顶点环绕方向（顺时针/逆时针）来明确内部区域，这样就不需要用户额外提供参考点了。下面我分享两个经过验证的Python实现方案：

方案一：基于球面环绕数（Winding Number）的方法 #

这个方法的核心是计算目标点相对于多边形的环绕数：如果环绕数为±1，则点在内部；为0则在外部。具体逻辑是通过累加目标点与多边形每条边形成的球面角方向总和，最终根据总和判断位置。

代码实现 #

先实现基础的球面几何工具函数：

import numpy as np

def normalize(v):
    """将向量归一化为单位向量"""
    norm = np.linalg.norm(v)
    return v / norm if norm > 1e-10 else v

def cross_product(u, v):
    """三维向量叉乘"""
    return np.array([
        u[1]*v[2] - u[2]*v[1],
        u[2]*v[0] - u[0]*v[2],
        u[0]*v[1] - u[1]*v[0]
    ])

def dot_product(u, v):
    """三维向量点乘"""
    return np.dot(u, v)

def spherical_angle_sign(p, v1, v2):
    """计算点p相对于边v1->v2的球面角符号（+1逆时针，-1顺时针，0共线）"""
    v1_p = cross_product(v1, p)
    v1_v2 = cross_product(v1, v2)
    
    # 基于右手定则判断方向
    sign = np.sign(dot_product(p, cross_product(v1, v2)))
    if abs(sign) < 1e-10:
        return 0
    
    # 计算角度贡献值
    angle = np.arctan2(np.linalg.norm(cross_product(p, v2)), dot_product(v1, p))
    return sign * angle

然后是核心的包含测试函数：

def point_in_spherical_polygon(p, polygon, tol=1e-8):
    """
    判断点p是否在球面多边形内部
    参数：
        p: 目标点的三维坐标（单位向量）
        polygon: 多边形顶点的三维坐标列表，每个顶点都是单位向量
        tol: 数值容忍度
    返回：
        True（内部）/ False（外部）/ None（共线在边上/顶点上）
    """
    p = normalize(p)
    polygon = [normalize(v) for v in polygon]
    n = len(polygon)
    total_angle = 0.0
    
    for i in range(n):
        v1 = polygon[i]
        v2 = polygon[(i+1)%n]
        
        # 处理点与顶点重合的情况
        if np.linalg.norm(v1 - p) < tol or np.linalg.norm(v2 - p) < tol:
            return None
        
        sign_angle = spherical_angle_sign(p, v1, v2)
        if sign_angle == 0:
            # 检查点是否在边的大圆弧线段上
            angle_v1_p = np.arccos(dot_product(v1, p))
            angle_p_v2 = np.arccos(dot_product(p, v2))
            angle_v1_v2 = np.arccos(dot_product(v1, v2))
            if abs(angle_v1_p + angle_p_v2 - angle_v1_v2) < tol:
                return None
            continue
        
        total_angle += sign_angle
    
    # 环绕数判断：总角度接近±2π则在内部，接近0则在外部
    if abs(total_angle - 2*np.pi) < tol or abs(total_angle + 2*np.pi) < tol:
        return True
    elif abs(total_angle) < tol:
        return False
    else:
        # 数值误差导致的模糊情况
        return None

方案二：基于球面射线法的改进 #

平面射线法通过统计射线穿过边的次数判断内外，在球面上我们可以定义一条从目标点出发的大圆弧射线，统计它与多边形边的交点数：奇数次在内部，偶数次在外部。

代码实现 #

先添加大圆弧交点计算函数：

def great_circle_intersection(a1, a2, b1, b2):
    """计算两条大圆弧a1-a2和b1-b2的交点（返回单位向量，无交点则返回None）"""
    # 计算两个大圆弧所在平面的法向量
    n1 = cross_product(a1, a2)
    n2 = cross_product(b1, b2)
    # 交点为两个法向量的叉乘
    intersect = normalize(cross_product(n1, n2))
    
    # 验证交点是否在两条大圆弧的线段范围内
    def on_segment(p, v1, v2):
        angle_v1_p = np.arccos(dot_product(v1, p))
        angle_p_v2 = np.arccos(dot_product(p, v2))
        angle_v1_v2 = np.arccos(dot_product(v1, v2))
        return abs(angle_v1_p + angle_p_v2 - angle_v1_v2) < 1e-8
    
    if on_segment(intersect, a1, a2) and on_segment(intersect, b1, b2):
        return intersect
    # 检查球面另一侧的交点
    intersect_neg = -intersect
    if on_segment(intersect_neg, a1, a2) and on_segment(intersect_neg, b1, b2):
        return intersect_neg
    return None

然后是射线法的包含测试函数：

def point_in_spherical_polygon_raycast(p, polygon, tol=1e-8):
    """
    用射线法判断点是否在球面多边形内部
    参数：
        p: 目标点的三维坐标（单位向量）
        polygon: 多边形顶点的三维坐标列表，每个顶点都是单位向量
        tol: 数值容忍度
    返回：
        True（内部）/ False（外部）/ None（共线在边上/顶点上）
    """
    p = normalize(p)
    polygon = [normalize(v) for v in polygon]
    n = len(polygon)
    
    # 生成一条与p垂直的射线方向，避免与多边形顶点共线
    if abs(p[0]) > abs(p[1]):
        ray_dir = normalize(np.array([-p[2], 0, p[0]]))
    else:
        ray_dir = normalize(np.array([0, p[2], -p[1]]))
    ray_end = normalize(p + ray_dir)
    
    intersection_count = 0
    
    for i in range(n):
        v1 = polygon[i]
        v2 = polygon[(i+1)%n]
        
        # 处理点与顶点或边重合的情况
        if np.linalg.norm(v1 - p) < tol or np.linalg.norm(v2 - p) < tol:
            return None
        angle_v1_p = np.arccos(dot_product(v1, p))
        angle_p_v2 = np.arccos(dot_product(p, v2))
        angle_v1_v2 = np.arccos(dot_product(v1, v2))
        if abs(angle_v1_p + angle_p_v2 - angle_v1_v2) < tol:
            return None
        
        # 计算射线与边的交点
        intersect = great_circle_intersection(p, ray_end, v1, v2)
        if intersect is not None and np.linalg.norm(intersect - p) > tol:
            intersection_count += 1
    
    # 奇数次交点在内部，偶数次在外部
    return intersection_count % 2 == 1

使用补充说明 #

1. 如果你的输入是经纬度，可通过以下函数转换为单位球面三维坐标：

def latlon_to_cartesian(lat, lon):
    """将纬度（弧度）、经度（弧度）转换为单位球面上的三维坐标"""
    x = np.cos(lat) * np.cos(lon)
    y = np.cos(lat) * np.sin(lon)
    z = np.sin(lat)
    return np.array([x, y, z])

2. 两个方案都依赖多边形顶点的环绕方向来定义内部，无需用户额外提供参考点；
3. 返回None表示点处于多边形边界（顶点或边），可根据业务需求自定义处理逻辑。

内容的提问来源于stack exchange，提问作者Shek