基于多边形长边边界与内部点计算风向的方法优化咨询

现有一个多边形，需基于其长边边界与内部点计算风向——内部点为风的起始点，风向遵循从长边到短边的长距离方向。

当前实现思路：

• 找出多边形的长边边界；
• 找到边界LineString上距离风起始点最近的点；
• 计算风向。

请问是否存在更优的实现方法？

当前实现代码 #

from shapely import Point, Polygon
from shapely.geometry import LineString
import numpy as np
import matplotlib.pyplot as plt

coords = ((-1, 0.), (0., 1.), (0.5, 1.), (-0.5, 0.), (-1, 0.))
p = Polygon(coords)

p1 = (0, 0.9)
p2 = (-0.5, 0.2)

fig, axs = plt.subplots()
axs.plot(*p.exterior.xy)
axs.scatter(p1[0], p1[1], c='r')
axs.quiver(p1[0], p1[1], -0.5, -0.5, angles='xy', scale_units='xy', scale=1, color='r')
axs.quiver(p2[0], p2[1], 0.5, 0.5, angles='xy', scale_units='xy', scale=1, color='b')

axs.scatter(p2[0], p2[1], c='b')
axs.set_xlim(-2, 2)
axs.set_ylim(-2, 2)

def calculate_wind_direction(p1, p2):
    dx = p2[0] - p1[0]
    dy = p2[1] - p1[1]
    
    # 计算弧度角度
    angle = np.arctan2(dy, dx)
    
    # 转换为角度
    degrees = np.degrees(angle)
    
    # 调整为气象学惯例（顺时针从北开始）
    meteorological_angle = (90 - degrees) % 360
    
    # 转换为"来向"
    from_direction = (meteorological_angle + 180) % 360
    
    return from_direction

# 提取多边形边界为线
l = p.boundary
coords = [c for c in l.coords]
# 生成各边线段
segments = [LineString([a, b]) for a, b in zip(coords, coords[1:])]
# 找到最长边
longest_segment = max(segments, key=lambda x: x.length)

p_start, p_end = [c for c in longest_segment.coords]

# 根据端点距离内部点的远近设置风向
if Point(p_start).distance(Point(p1)) < Point(p_end).distance(Point(p1)):
    wdir = calculate_wind_direction(p_start, p_end)
else:
    wdir = calculate_wind_direction(p_end, p_start)

风向示意图 #

更优实现思路及改进点 #

1. 长边匹配逻辑优化
   当前仅取最长单一边段，但多边形可能存在多条长度接近的"长边"，建议先设定长度阈值（比如取最长边长度的80%为阈值），筛选出所有候选长边；再根据内部点位置，选择距离内部点最远的那条长边（即对应对向短边的长边），避免误选相邻的长侧边。

2. 最近点计算精度提升
   当前仅比较长边两端点到内部点的距离，实际应计算内部点到整条长边线段的垂足点（用Shapely的segment.interpolate(segment.project(Point(inner_point)))可直接获取），这个垂足点才是长边到内部点的真实最近点，以此为参考计算的风向更贴合"从长边指向短边"的逻辑。

3. 风向计算鲁棒性增强

   • 增加内部点的多边形内判断，避免无效计算；
   • 风向计算改为基于"长边垂足点指向内部点"的反方向，再转换为气象学"来向"标准（顺时针从北开始），结果更符合实际风向定义。

4. 代码简化与性能优化

   • 提取边段可直接用list(p.boundary.coords)简化流程；
   • 批量计算时预先缓存候选长边集合，减少重复计算；
   • 用生成器表达式替代部分列表推导，降低内存占用。

改进后的示例代码 #

from shapely import Point, Polygon
from shapely.geometry import LineString
import numpy as np
import matplotlib.pyplot as plt

coords = ((-1, 0.), (0., 1.), (0.5, 1.), (-0.5, 0.), (-1, 0.))
p = Polygon(coords)
p1 = (0, 0.9)
p2 = (-0.5, 0.2)

# 绘制可视化基础图
fig, axs = plt.subplots()
axs.plot(*p.exterior.xy)
axs.scatter(p1[0], p1[1], c='r')
axs.scatter(p2[0], p2[1], c='b')
axs.set_xlim(-2, 2)
axs.set_ylim(-2, 2)

def calculate_wind_direction(from_point, to_point):
    """计算气象学标准的风来向（顺时针从北开始，范围0-360°）"""
    dx = to_point[0] - from_point[0]
    dy = to_point[1] - from_point[1]
    angle_rad = np.arctan2(dy, dx)
    angle_deg = np.degrees(angle_rad)
    # 转换为气象来向：将数学角度转换为顺时针从北，再取反得到来向
    meteorological_from = (90 - angle_deg + 180) % 360
    return meteorological_from

# 提取多边形所有边段
boundary_coords = list(p.boundary.coords)
segments = [LineString([boundary_coords[i], boundary_coords[i+1]]) 
            for i in range(len(boundary_coords)-1)]

# 筛选候选长边（长度不小于最长边的80%）
max_length = max(seg.length for seg in segments)
long_segments = [seg for seg in segments if seg.length >= max_length * 0.8]

def get_optimal_wind_dir(inner_point, long_segments, polygon):
    """计算内部点的最优风向"""
    inner_pt = Point(inner_point)
    if not polygon.contains(inner_pt):
        raise ValueError("内部点不在多边形范围内")
    
    # 选择距离内部点最远的长边（对应对向短边的长边）
    target_segment = max(long_segments, key=lambda seg: seg.distance(inner_pt))
    
    # 获取长边到内部点的垂足点
    proj_distance = target_segment.project(inner_pt)
    foot_point = target_segment.interpolate(proj_distance)
    
    # 计算风向向量（从长边指向内部点，即风的去向）
    wind_vector = (inner_pt.x - foot_point.x, inner_pt.y - foot_point.y)
    # 获取气象学风来向
    wind_dir = calculate_wind_direction(foot_point, inner_pt)
    
    # 绘制风向箭头
    axs.quiver(inner_pt.x, inner_pt.y, wind_vector[0], wind_vector[1], 
               angles='xy', scale_units='xy', scale=1, color='r' if inner_point == p1 else 'b')
    return wind_dir

# 计算两个点的风向
wdir_p1 = get_optimal_wind_dir(p1, long_segments, p)
wdir_p2 = get_optimal_wind_dir(p2, long_segments, p)
print(f"点p1的气象风来向：{wdir_p1:.1f}°")
print(f"点p2的气象风来向：{wdir_p2:.1f}°")

plt.show()

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