如何用Python围绕点簇绘制非重叠多边形并聚合至GeoDataFrame？

首先修正你初始代码中的两处问题：一是gpd.points_from_xy调用时用了不存在的Long/Lat列（实际是小写的lon/lat）；二是直接为经纬度坐标设置投影坐标系ESRI:102003，正确流程应先指定原始坐标为WGS84（EPSG:4326）再投影。

修正后的基础代码 #

import pandas as pd
import numpy as np
import geopandas as gpd
from sklearn.cluster import DBSCAN
from shapely.geometry import Polygon
from shapely.ops import unary_union

# 初始化数据
df = pd.DataFrame({
    'yr': [2018, 2017, 2018, 2016],
    'id': [0, 1, 2, 3],
    'v': [10, 12, 8, 10],
    'lat': [32.7418248, 32.8340583, 32.8340583, 32.7471895],
    'lon': [-97.524066, -97.0805484, -97.0805484, -96.9400779]
})

# 创建GeoDataFrame并完成坐标投影
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df['lon'], df['lat']),
    crs="EPSG:4326"  # 原始经纬度坐标的CRS
)
gdf = gdf.to_crs("ESRI:102003")  # 转换为适合距离计算的投影CRS

步骤1：识别空间点簇 #

用DBSCAN算法将空间上邻近的点归为同一簇，参数eps（聚类半径）和min_samples（簇内最小点数）可根据你的数据调整：

# 提取点坐标作为聚类输入
coords = np.array([(geom.x, geom.y) for geom in gdf.geometry])

# 执行DBSCAN聚类
dbscan = DBSCAN(eps=10000, min_samples=1)
gdf['cluster_id'] = dbscan.fit_predict(coords)

步骤2：生成自定义多边形（两种可选方案） #

方案1：Alpha形状（贴合点簇的不规则多边形） #

Alpha形状能生成贴合点集分布的不规则多边形，比凸包更灵活，且确保顶点数≥5：

def alpha_shape(points, alpha=0.0001):
    """基于点集生成Alpha形状多边形"""
    from shapely.ops import triangulate
    from shapely.geometry import MultiPoint
    
    if len(points) < 3:
        return None
    
    multi_point = MultiPoint(points)
    # 生成Delaunay三角网
    triangles = triangulate(multi_point)
    alpha_sq = alpha ** 2
    edges = set()
    
    # 筛选符合Alpha阈值的三角形边
    for tri in triangles:
        vertices = list(tri.exterior.coords)[:-1]
        dists = [
            (vertices[i][0]-vertices[j][0])**2 + (vertices[i][1]-vertices[j][1])**2
            for i, j in [(0,1), (1,2), (2,0)]
        ]
        if all(d <= alpha_sq for d in dists):
            for i in range(3):
                edge = tuple(sorted([vertices[i], vertices[(i+1)%3]]))
                edges.add(edge)
    
    # 构建多边形环
    if not edges:
        return multi_point.convex_hull
    
    edge_map = {}
    for p1, p2 in edges:
        edge_map.setdefault(p1, []).append(p2)
        edge_map.setdefault(p2, []).append(p1)
    
    start = next(iter(edge_map.keys()))
    current, next_p = start, edge_map[start][0]
    poly_coords = [start]
    
    while next_p != start:
        poly_coords.append(next_p)
        neighbors = edge_map[next_p]
        next_p = neighbors[0] if neighbors[1] == current else neighbors[1]
        current = poly_coords[-1]
    
    poly_coords.append(start)
    
    # 若顶点数不足5，用凸包替代
    return Polygon(poly_coords) if len(poly_coords) >=5 else multi_point.convex_hull

# 为每个簇生成Alpha形状多边形
cluster_polys = {}
for cluster_id in gdf['cluster_id'].unique():
    cluster_points = gdf[gdf['cluster_id'] == cluster_id].geometry
    point_coords = [(p.x, p.y) for p in cluster_points]
    cluster_polys[cluster_id] = alpha_shape(point_coords)

# 将多边形匹配到每个点的新列
gdf['cluster_polygon'] = gdf['cluster_id'].map(cluster_polys)

方案2：正N边形（规则非圆形多边形） #

如果需要规则形状的多边形，可生成边数≥5的正多边形，围绕簇中心扩展：

def regular_polygon(center, radius, num_sides=5):
    """生成围绕中心点的正N边形"""
    from math import pi, cos, sin
    
    angles = [2 * pi * i / num_sides for i in range(num_sides)]
    coords = [
        (center.x + radius * cos(angle), center.y + radius * sin(angle))
        for angle in angles
    ]
    coords.append(coords[0])  # 闭合多边形
    return Polygon(coords)

# 为每个簇生成正5边形
cluster_polys = {}
for cluster_id in gdf['cluster_id'].unique():
    cluster_gdf = gdf[gdf['cluster_id'] == cluster_id]
    # 计算簇中心
    cluster_center = cluster_gdf.geometry.unary_union.centroid
    # 计算簇内点到中心的最大距离，作为基础半径并添加缓冲
    max_dist = max(cluster_gdf.geometry.distance(cluster_center))
    poly = regular_polygon(cluster_center, max_dist + 2000, num_sides=5)
    cluster_polys[cluster_id] = poly

# 匹配多边形到数据列
gdf['cluster_polygon'] = gdf['cluster_id'].map(cluster_polys)

步骤3：处理多边形重叠 #

若生成的多边形存在重叠，可通过以下方式处理：

# 合并所有多边形并拆分非重叠部分
all_polys = gpd.GeoSeries(cluster_polys.values())
union = unary_union(all_polys)
# 若合并后为MultiPolygon，可拆分后重新匹配到簇（需根据业务逻辑调整）
if union.geom_type == 'MultiPolygon':
    split_polys = list(union.geoms)
    # 此处需添加逻辑将拆分后的多边形与原簇匹配

最终结果 #

处理后的gdf新增了cluster_id（簇标识）和cluster_polygon（对应多边形）列，可直接用于分析或导出：

# 查看结果
print(gdf[['id', 'cluster_id', 'cluster_polygon']])

# 可选：转换回WGS84坐标系便于查看
gdf_wgs84 = gdf.to_crs("EPSG:4326")

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