Python中求解覆盖多边形指定比例的圆形标记尺寸/半径

核心逻辑是先建立markersize与实际地理半径的映射关系，再通过迭代法计算圆形与多边形的重叠面积，逐步调整markersize直至达到目标覆盖比例。

关键步骤 #

• 计算多边形地理面积：借助地理空间库（如geopandas）精准获取目标多边形的实际面积。
• 建立尺寸转换规则：将地理半径通过绘图轴的比例尺、屏幕DPI，转换为matplotlib的markersize（注意markersize是标记面积，单位为点的平方）。
• 迭代调整+面积校验：用二分法不断调整markersize，计算圆形与多边形的重叠面积占比，直到符合设定的目标比例。

完整实现代码 #

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

def markersize_to_radius(markersize, ax, dpi=96):
    """将matplotlib的markersize转换为经纬度坐标系下的地理半径（单位：度）"""
    # markersize是标记面积（点²），1点=1/72英寸
    inch_area = markersize / (72 ** 2)
    pixel_area = inch_area * (dpi ** 2)
    radius_pixel = np.sqrt(pixel_area / np.pi)
    
    # 获取轴的坐标范围，计算像素与地理坐标的转换比例
    x_min, x_max = ax.get_xlim()
    y_min, y_max = ax.get_ylim()
    ax_width_pixel = ax.get_window_extent().width
    ax_height_pixel = ax.get_window_extent().height
    
    lon_per_pixel = (x_max - x_min) / ax_width_pixel
    lat_per_pixel = (y_max - y_min) / ax_height_pixel
    
    # 取经纬度方向的平均比例作为近似转换系数
    radius_deg = radius_pixel * np.mean([lon_per_pixel, lat_per_pixel])
    return radius_deg

def radius_to_markersize(radius_deg, ax, dpi=96):
    """将经纬度坐标系下的地理半径（度）转换为matplotlib的markersize"""
    x_min, x_max = ax.get_xlim()
    y_min, y_max = ax.get_ylim()
    ax_width_pixel = ax.get_window_extent().width
    ax_height_pixel = ax.get_window_extent().height
    
    lon_per_pixel = (x_max - x_min) / ax_width_pixel
    lat_per_pixel = (y_max - y_min) / ax_height_pixel
    pixel_per_deg = 1 / np.mean([lon_per_pixel, lat_per_pixel])
    
    radius_pixel = radius_deg * pixel_per_deg
    pixel_area = np.pi * (radius_pixel ** 2)
    inch_area = pixel_area / (dpi ** 2)
    markersize = inch_area * (72 ** 2)
    return markersize

def find_target_markersize(epicenter, polygon, target_ratio=0.2, ax=None, tol=1e-3, max_iter=100):
    """找到使圆形覆盖多边形目标比例的markersize"""
    if ax is None:
        # 临时绘图初始化轴范围
        fig, ax = plt.subplots()
        gpd.GeoSeries([polygon]).plot(ax=ax)
        ax.scatter(epicenter.x, epicenter.y, s=1)
        plt.close(fig)
    
    # 计算多边形总面积
    poly_area = polygon.area
    
    # 初始化二分法搜索范围
    low_ms = 1
    high_ms = 10000
    best_ms = None
    best_ratio = 0
    
    for _ in range(max_iter):
        current_ms = (low_ms + high_ms) / 2
        radius_deg = markersize_to_radius(current_ms, ax)
        # 创建震中对应的圆形几何对象
        circle = epicenter.buffer(radius_deg)
        # 计算圆形与多边形的重叠面积
        intersection = polygon.intersection(circle)
        overlap_area = intersection.area
        current_ratio = overlap_area / poly_area
        
        # 调整搜索区间
        if current_ratio < target_ratio:
            low_ms = current_ms
        else:
            high_ms = current_ms
        
        # 检查是否达到精度要求
        if abs(current_ratio - target_ratio) < tol:
            best_ms = current_ms
            best_ratio = current_ratio
            break
    
    return best_ms, best_ratio

# 示例使用
if __name__ == "__main__":
    # 模拟多边形（经纬度坐标）
    polygon_coords = [(100.0, 30.0), (100.5, 30.0), (100.5, 30.5), (100.0, 30.5)]
    polygon = Polygon(polygon_coords)
    
    # 模拟震中坐标
    epicenter = Point(100.25, 30.25)
    
    # 创建绘图轴
    fig, ax = plt.subplots(figsize=(8, 8))
    gpd.GeoSeries([polygon]).plot(ax=ax, color='blue', alpha=0.5)
    
    # 寻找目标markersize
    target_ms, actual_ratio = find_target_markersize(epicenter, polygon, target_ratio=0.2, ax=ax)
    
    # 绘制结果
    ax.scatter(epicenter.x, epicenter.y, color='green', s=target_ms, alpha=0.5)
    ax.set_title(f"实际覆盖比例: {actual_ratio:.2%} | Markersize: {target_ms:.0f}")
    plt.show()

注意事项 #

• 坐标适配：如果使用投影坐标系（如UTM），可直接用x方向的像素-地理比例转换，无需经纬度平均，精度更高。
• 重叠处理：通过shapely的intersection方法直接计算重叠面积，天然适配多边形的复杂形状，无需额外处理重叠逻辑。
• 精度控制：调整tol参数可控制覆盖比例的精度，max_iter参数限制迭代次数，避免无限循环。

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