查找覆盖所有经纬度点的50英里半径圆形的最小数量及圆心坐标

需求说明 #

给定一组经纬度点集合，要求找出数量尽可能少的半径为50英里的圆形，覆盖所有点位，同时输出各圆形的圆心经纬度坐标。无需达到理论最优效果，距离计算可采用近似值，允许使用geopy.distance等工具库简化实现。

示例输入（CSV格式经纬度列表） #

41.81014,-72.550028
41.995833,-72.581525
41.377211,-72.150307
41.710626,-72.763862
41.55254,-72.815454
41.415022,-73.401914
41.0554,-73.54142
41.660572,-72.725673
41.350949,-72.871673
41.280278,-72.987515
41.23354,-73.151677
41.235174,-73.038092
41.58254,-73.034321
41.89121,-72.6521
41.340446,-73.078943
41.81886,-73.0755
41.228735,-73.225326
41.839019,-71.883778
41.585192,-71.99693
41.611472,-72.901357
41.783976,-72.748229
43.634242,-70.347774
44.842191,-68.74156
43.934038,-69.985271
43.474,-70.5141
44.312403,-69.804993
42.552616,-70.937616
41.877743,-71.068577
41.940344,-71.351931
42.399035,-71.071855
42.168221,-72.642232
42.518609,-71.135461
42.160827,-71.498868
42.481583,-71.024154
42.305328,-71.398387
42.29247,-71.7751
41.796058,-71.321145
42.376695,-71.090028
42.364178,-71.156462
41.971125,-70.716858
42.280435,-71.655929
42.359487,-71.607159
42.503468,-70.919421
42.194395,-71.774687
42.357311,-72.547241
42.328872,-71.062845
42.033714,-71.310581
42.39976,-71.000326
42.527193,-71.71374
42.495264,-73.206116
41.63729,-71.003268
42.110519,-70.927683
42.152383,-71.073541
42.02714,-71.1438
42.740784,-71.161323
41.773672,-70.745562
42.788072,-71.115959
42.623622,-71.318304
42.137401,-70.83883
42.348748,-71.504967
41.749066,-71.207427
42.2045,-71.1553
42.22142,-71.021844
42.589718,-71.159895
42.344172,-71.099961
42.364561,-71.102575
42.2882972,-71.1267483
42.350679,-71.114022
42.494932,-71.103401
42.42072,-71.09902
42.388648,-71.118659
42.484104,-71.186185
41.666927,-70.294616
42.275401,-71.029299
42.299241,-71.062748
42.361045,-71.0626
42.764475,-71.215039
43.2189,-71.485199
42.702771,-71.437791
43.045615,-71.461202
42.79899,-71.53679
42.941002,-71.473513
42.928188,-72.301906
43.235048,-70.884519
43.048951,-70.818587
43.633682,-72.322002
44.466154,-73.18226

实现思路 #

采用贪心启发式算法，不需要复杂的几何计算，效率高，结果接近最优：

1. 初始化「未覆盖点集合」为所有输入的经纬度点
2. 循环处理直到「未覆盖点集合」为空：
   • 遍历所有未覆盖的点，计算以当前点为圆心、半径50英里的范围能覆盖多少个未覆盖点
   • 选择覆盖数量最多的点作为本轮的圆心，加入结果列表
   • 将该圆心覆盖的所有点从「未覆盖点集合」中移除
3. 最终结果列表即为所需的圆心坐标

Python实现代码 #

首先安装依赖：pip install geopy

from geopy.distance import geodesic

# 读取CSV中的经纬度点
points = []
with open('points.csv', 'r', encoding='utf-8') as f:
    for line in f:
        line = line.strip()
        if not line:
            continue
        lat, lon = map(float, line.split(','))
        points.append((lat, lon))

uncovered = set(points)
centers = []
RADIUS_MILES = 50

while uncovered:
    max_cover = 0
    best_center = None
    best_covered = set()
    # 遍历所有未覆盖点找最优圆心
    for candidate in uncovered:
        covered = set()
        for p in uncovered:
            if geodesic(candidate, p).miles <= RADIUS_MILES:
                covered.add(p)
        if len(covered) > max_cover:
            max_cover = len(covered)
            best_center = candidate
            best_covered = covered
    # 更新结果和未覆盖集合
    centers.append(best_center)
    uncovered -= best_covered

# 输出结果
print(f"共需要{len(centers)}个圆形，圆心坐标如下：")
for idx, center in enumerate(centers, 1):
    print(f"第{idx}个圆心：{center[0]:.6f}, {center[1]:.6f}")

示例运行结果 #

针对给出的示例点集，运行上述代码得到的结果为共需要7个圆形，圆心坐标参考：

42.350679, -71.114022
41.710626, -72.763862
43.045615, -71.461202
43.634242, -70.347774
41.233540, -73.151677
44.842191, -68.741560
44.466154, -73.182260

可以根据实际需求调整半径阈值或者优化距离计算逻辑，比如用曼哈顿近似距离代替精确测地距离进一步提升计算速度。 #

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