如何在地理空间数据上应用带约束的K-Means聚类？

需求概述 #

我需要在Python中编写针对地理空间数据的带约束K-Means聚类程序，具体要求如下：

• 输入数据：包含经纬度坐标对、quantity约束列的GeoDataFrame
• 聚类规则：
  • 每个聚类内的点地理位置相近
  • 各聚类对应的多边形互不相交
  • 每个聚类的quantity列总和必须在15000至18000之间

数据相关参考：

• QGIS中点数据呈现自然地理分布特征
• GeoDataFrame前38条数据包含唯一ID、quantity数值、几何坐标信息
• 目标聚类效果为无重叠的多边形区域，每个区域内的点归属唯一聚类

现有问题 #

我自行编写了K-Means聚类调整脚本，尝试解决聚类总和低于下限（15000）的问题，但出现了总和超出上限（18000）的情况。

现有实现代码 #

n_territories = number_of_cluster(geodataframe=gdf, col_name=col_name, threshold_value=threshold_value)

if num_of_cluster is not None: 
    if num_of_cluster > n_territories:n_territories= int(num_of_cluster)
    #n_territories = int(num_of_cluster)
    # if n_territories<4: n_territories= 4
# else: n_territories = int(gdf[col_name].sum()/threshold_value)-2
# if n_territories > (int(gdf[col_name].sum()/threshold_value)-2): n_territories = int(gdf[col_name].sum()/threshold_value)-2  
# Generate 2000 unique random hex color codes with high contrast
# n_territories = int(gdf[col_name].sum()/threshold_value)


threshold_value1= threshold_value*1.2
#unique_high_contrast_color_codes = [random_high_contrast_color(used_color_codes) for _ in range(n_territories)]
unique_high_contrast_color_codes = generate_unique_complementary_colors(n_territories)
coordinates = np.array(gdf.geometry.apply(lambda geom: (geom.x, geom.y)).tolist())
gdf_array = gdf[['ttpl_id',col_name,'geometry']].values
points = gdf_array[:, 2]
first_point = min(points, key=lambda point: (point.x, point.y))
gdf_array = sorted(gdf_array, key=lambda row:  first_point.distance(row[2]))
gdf_array = np.vstack(gdf_array)
#Initializing the centroid 
centroids = kmeans_plusplus_init(coordinates,n_territories, see=see)
# Find the assignments of data points to centroids
assignments = assign_to_nearest_centroid(gdf_array, centroids)
gdf_array = np.column_stack((gdf_array,np.full(len(gdf_array), '0')))
# Create an array of values corresponding to the assignments
values = np.concatenate([np.full(len(assignment), i) for i, assignment in enumerate(assignments[0])])
# Assign the values to the appropriate positions in gdf_array
gdf_array[np.concatenate(assignments[0]), 3] = values
categories = np.unique(gdf_array[:, 3]) ###finding the unique cluster id (3rd Column) from each points 
sum_per_category = [np.sum(gdf_array[gdf_array[:, 3] == cat][:, 1].astype(int)) for cat in categories]  # Calculate the sum of the cluster id (unique value 3rd Column) for each category
result = np.column_stack((categories, sum_per_category)) # Combine the category labels and the sums into a new array
dist_val, sum_constraint = calculate_sse(gdf_array, centroids, threshold_value)
dist_old= 1
sum_constraint_old= 1
dist_tol= abs(dist_val-dist_old)/dist_old
const_tol= abs(sum_constraint-sum_constraint_old)/sum_constraint_old
result1= result[result[:,-1]<threshold_value]
result2= result[result[:,-1]>threshold_value]
it_= 0
while (dist_tol> 0.0001 or const_tol>0.0001) and it_<20:
    arr1 = np.empty_like(gdf_array[:1,:])
    for item in result2:
        id_ = item[0]
        gdf_array = gdf_array[np.unique(gdf_array[:, 0], return_index=True)[1]]
        gdf_array = sorted(gdf_array, key=lambda row:  first_point.distance(row[2]))
        gdf_array = np.vstack(gdf_array)
        gdf1_array1= gdf_array[gdf_array[:,3]==id_]
        points = gdf1_array1[:, 2]
        median_point = Point(np.median([point.x for point in points]), np.median([point.y for point in points]))
        distances = np.array([round(haversine_distance(row[2], centroids[id_]),2) for row in gdf1_array1])
        gdf1_array1 = np.column_stack((gdf1_array1, distances))
        gdf1_array1 = gdf1_array1[np.argsort(gdf1_array1[:, -1])]
        cumulative_sum = np.cumsum(gdf1_array1[:, 1])
        gdf1_array2 = gdf1_array1[:np.argmax(cumulative_sum > threshold_value1),:]
        gdf3_array = gdf1_array1[np.isin(gdf1_array1[:, 0], gdf1_array2[:,0], invert=True)]
        gdf1_array1 = gdf1_array1[np.isin(gdf1_array1[:, 0], gdf1_array2[:,0], invert=False)]
        if np.sum(gdf1_array2[:, 1])>=threshold_value1:
            centroids[id_] = MultiPoint(gdf1_array1[:,2]).centroid
    for item in result1:
        try:
            id_ = item[0]
            cluster_sum = item[-1]
            gdf_array = gdf_array[np.unique(gdf_array[:, 0], return_index=True)[1]]
            gdf_array1= gdf_array
            distances = np.array([round(haversine_distance(row[2], centroids[id_]),2) for row in gdf_array1])
            gdf_array1 = np.column_stack((gdf_array1, distances))
            gdf_array1 = gdf_array1[np.argsort(gdf_array1[:, -1])]
            #gdf_array1 = np.vstack(sorted(gdf_array1,key=lambda row:  centroids[id_].distance(row[2]) ))
            cumulative_sum = np.cumsum(gdf_array1[:, 1])
            gdf_array2 = gdf_array1[:np.argmax(cumulative_sum > (threshold_value1)),:]
            gdf_array2[:, 3] = id_
            gdf_array2= np.concatenate((gdf_array2[:,:-1], gdf_array[gdf_array[:,3]==id_]), axis=0)
            gdf_array = gdf_array[np.isin(gdf_array[:, 0], gdf_array2[:,0], invert=True)]
            arr1 = np.concatenate((arr1, gdf_array2), axis=0)
        except:pass
    gdf_array = np.concatenate((arr1[1:, :], gdf_array), axis=0)
    # Extract the first column
    gdf_array = gdf_array[np.unique(gdf_array[:, 0], return_index=True)[1]]
    centroids = np.empty((n_territories,), dtype=object)
    for j in range(0, n_territories):
        centroids[j] = MultiPoint(gdf_array[gdf_array[:,3]==j][:,2]).centroid
    ############condition checking:
    sum_per_category = [np.sum(gdf_array[gdf_array[:, 3] == cat][:, 1].astype(int)) for cat in categories]
    # Combine the category labels and the sums into a new array
    result = np.column_stack((result, sum_per_category))
    result1= result[result[:,-1]<threshold_value]
    result2= result[result[:,-1]>threshold_value]
    dist_old= dist_val
    sum_constraint_old= sum_constraint
    dist_val, sum_constraint = calculate_sse(gdf_array, centroids, threshold_value)
    dist_tol= abs(dist_val-dist_old)/dist_old
    const_tol= abs(sum_constraint-sum_constraint_old)/sum_constraint_old
    it_= it_+1

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