找到覆盖所有点的 50 英里半径的最小圆圈数
Find Minimum number of Circles with 50-Mile Radius that Covers All Points
给定一个包含 lat/lon 个点的列表，我们如何找到最小数量的 50 英里半径圆（及其 lat/lon 个点），使这些圆覆盖列表中的所有点？
解不需要最优，radiuses/distances的计算可以近似计算，简单。或者使用像 geopy.distance.
这样的辅助库
例如，这里是 lat/lon 个点的 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
根据评论更新的答案：
你有很多选择。
这里有 3 种不同的方法：
1.随着 scipy.CKDTree:
优点：
- 这会很快
缺点：
- 不太准确，因为计算的距离是欧氏距离
- 并且半径将与您输入的相同，所以这里以度为单位
我会使用 cKDTree 和半径查询来查找半径内的所有点，从列表中删除这些点，然后继续剩余的点。这不是最优的，但可以作为一个很好的基础。
from scipy.spatial import cKDTree
points = [(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)]
# Radius of circle. Note that the unit is the same as in your list (here, degrees.)
radius = 1
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = points[0]
# Create a KDTree
tree = cKDTree(points)
# Find indexes of all points in radius
indexes_of_points_in_radius = tree.query_ball_point(start_point, radius)
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
2。随着 sklearn.neighbors.BallTree:
优点：
- 这样会更准确，因为这里计算的距离是大圆距离（Haversine公式）。
缺点：
- 与 cKDTree 一样，半径将与您的输入相同，因此这里以度为单位。
- 比scipy.cKDTree慢一点（我测试时慢了2倍）
另请注意，我发现有些人建议将您的输入转换为弧度，因为这是 Haversine 公式 (https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html) 所必需的。但根据我对 scikit-learn 1.0.1 的测试，这不是必需的。但以防万一，您可以这样做：
from math import radians
points = [tuple(map(radians, point)) for point in points]
start_point = (radians(start_point[0]), radians(start_point[1]))
radius = radians(radius)
使用 BallTree 的代码：
from sklearn.neighbors import BallTree
import numpy as np
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = np.array([points[0]])
# Create a BallTree, and chose the Haversine formula (great circle distance)
tree = BallTree(points, metric="haversine")
# Find indexes of all points in radius
indexes_of_points_in_radius = tree.query_radius(start_point, r=radius)[0]
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
3。使用sklearn.neighbors.BallTree，使用用户定义的距离函数：
优点：
- 我们将能够在这里使用非常精确的距离
- 我们将能够以英里（或米）为单位指定此距离
缺点：
- 比 cKDTree 慢很多（我测试时是 10 倍）
from pyproj import Geod
from sklearn.neighbors import BallTree
import numpy as np
# Create a WGS84 ellipsoid
geod = Geod(ellps='WGS84')
# User-defined function for BallTree
# We use the "inv" method of pyproj in order to get the distance in meters between 2 points
# It computes the geodesic distance using the wonderful Karney's algorithm
def geodedsic_distance(point_01, point_02):
lat1,lon1 = point_01
lat2,lon2 = point_02
az12,az21,distance_in_meters = geod.inv(lon1,lat1,lon2,lat2)
return distance_in_meters
def miles_to_meters(miles):
return miles * 1609.344
# Radius in miles
radius_in_miles = 50
radius_in_meters = miles_to_meters(50)
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = np.array([points[0]])
# Create a BallTree, and chose our custom function
tree = BallTree(points, metric=geodedsic_distance)
# Find indexes of all points in radius, specified in meters
indexes_of_points_in_radius = tree.query_radius(start_point, r=radius_in_meters)[0]
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
如果您想了解更多关于英里到度的转换（以及为什么，事实上，我们不能）和计算地球上的距离：
给定一个包含 lat/lon 个点的列表，我们如何找到最小数量的 50 英里半径圆（及其 lat/lon 个点），使这些圆覆盖列表中的所有点？
解不需要最优，radiuses/distances的计算可以近似计算，简单。或者使用像 geopy.distance.
例如，这里是 lat/lon 个点的 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
根据评论更新的答案：
你有很多选择。
这里有 3 种不同的方法：
1.随着 scipy.CKDTree:
优点：
- 这会很快
缺点：
- 不太准确，因为计算的距离是欧氏距离
- 并且半径将与您输入的相同，所以这里以度为单位
我会使用 cKDTree 和半径查询来查找半径内的所有点，从列表中删除这些点，然后继续剩余的点。这不是最优的，但可以作为一个很好的基础。
from scipy.spatial import cKDTree
points = [(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)]
# Radius of circle. Note that the unit is the same as in your list (here, degrees.)
radius = 1
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = points[0]
# Create a KDTree
tree = cKDTree(points)
# Find indexes of all points in radius
indexes_of_points_in_radius = tree.query_ball_point(start_point, radius)
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
2。随着 sklearn.neighbors.BallTree:
优点：
- 这样会更准确，因为这里计算的距离是大圆距离（Haversine公式）。
缺点：
- 与 cKDTree 一样，半径将与您的输入相同，因此这里以度为单位。
- 比scipy.cKDTree慢一点（我测试时慢了2倍）
另请注意，我发现有些人建议将您的输入转换为弧度，因为这是 Haversine 公式 (https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html) 所必需的。但根据我对 scikit-learn 1.0.1 的测试，这不是必需的。但以防万一，您可以这样做：
from math import radians
points = [tuple(map(radians, point)) for point in points]
start_point = (radians(start_point[0]), radians(start_point[1]))
radius = radians(radius)
使用 BallTree 的代码：
from sklearn.neighbors import BallTree
import numpy as np
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = np.array([points[0]])
# Create a BallTree, and chose the Haversine formula (great circle distance)
tree = BallTree(points, metric="haversine")
# Find indexes of all points in radius
indexes_of_points_in_radius = tree.query_radius(start_point, r=radius)[0]
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
3。使用sklearn.neighbors.BallTree，使用用户定义的距离函数：
优点：
- 我们将能够在这里使用非常精确的距离
- 我们将能够以英里（或米）为单位指定此距离
缺点：
- 比 cKDTree 慢很多（我测试时是 10 倍）
from pyproj import Geod
from sklearn.neighbors import BallTree
import numpy as np
# Create a WGS84 ellipsoid
geod = Geod(ellps='WGS84')
# User-defined function for BallTree
# We use the "inv" method of pyproj in order to get the distance in meters between 2 points
# It computes the geodesic distance using the wonderful Karney's algorithm
def geodedsic_distance(point_01, point_02):
lat1,lon1 = point_01
lat2,lon2 = point_02
az12,az21,distance_in_meters = geod.inv(lon1,lat1,lon2,lat2)
return distance_in_meters
def miles_to_meters(miles):
return miles * 1609.344
# Radius in miles
radius_in_miles = 50
radius_in_meters = miles_to_meters(50)
num_circles = 0
list_is_no_empty = True
while(list_is_no_empty):
# Take the first point in order to find all points within distance radius
start_point = np.array([points[0]])
# Create a BallTree, and chose our custom function
tree = BallTree(points, metric=geodedsic_distance)
# Find indexes of all points in radius, specified in meters
indexes_of_points_in_radius = tree.query_radius(start_point, r=radius_in_meters)[0]
# Create the list of points to remove (points that were found within distance radius)
points_to_remove = [points[i] for i in indexes_of_points_in_radius]
# Remove these points
points = list(set(points) - set(points_to_remove))
# Increment the number of circles
num_circles += 1
# If no points remain, exit loop
if points == []:
list_is_no_empty = False
print("Number of circles:", num_circles)
如果您想了解更多关于英里到度的转换（以及为什么，事实上，我们不能）和计算地球上的距离：