如何计算Shapefile中各区域内的最长距离？

当然可以，你可以用geopandas、shapely和numpy组合完成这个需求，以下是具体步骤和代码：

核心思路 #

1. 投影转换：原始Shapefile的经纬度坐标（WGS84）不适合直接计算距离，需要转成平面投影（比如UTM），确保距离计算的准确性（单位为米）。
2. 凸包简化：凹多边形的最远点对必然在其凸包上，先计算凸包可以大幅减少计算量。
3. 高效计算：用旋转卡壳算法快速找到凸包的最远点对（即区域的最长距离），比暴力遍历点对效率高得多。

完整代码 #

import geopandas as gpd
import numpy as np
from shapely.ops import unary_union

# 加载Shapefile数据
shape = gpd.GeoDataFrame.from_file("India-Districts-2011Census.shp")

# 转换到适合印度的UTM投影（EPSG:32643，覆盖印度北部，若数据跨多带可按需调整）
shape = shape.to_crs(epsg=32643)

def calculate_max_distance(geometry):
    # 处理MultiPolygon类型（如带飞地的区域）
    if geometry.geom_type == "MultiPolygon":
        # 合并所有子多边形的凸包，再取整体凸包
        convex_hulls = [poly.convex_hull for poly in geometry.geoms]
        combined_hull = unary_union(convex_hulls).convex_hull
        geometry = combined_hull
    
    convex_hull = geometry.convex_hull
    
    # 处理特殊几何类型：点或线
    if convex_hull.geom_type == "Point":
        return 0.0
    elif convex_hull.geom_type == "LineString":
        return convex_hull.length
    
    # 提取凸包顶点坐标（去掉最后一个重复的闭合点）
    coords = np.array(convex_hull.exterior.coords)[:-1]
    n = len(coords)
    
    if n == 2:
        return np.linalg.norm(coords[0] - coords[1])
    
    max_dist = 0.0
    j = 1
    # 旋转卡壳算法遍历凸包顶点，找到最远点对
    for i in range(n):
        edge_vector = coords[(i+1) % n] - coords[i]
        # 找到当前边对应的最远顶点
        while True:
            cross_product = np.cross(edge_vector, coords[(j+1) % n] - coords[j])
            if cross_product <= 0:
                break
            j = (j + 1) % n
        # 计算当前顶点对的距离
        current_dist = np.linalg.norm(coords[i] - coords[j])
        if current_dist > max_dist:
            max_dist = current_dist
    
    return max_dist

# 为每个区域计算最长距离（单位：米）
shape["max_distance_m"] = shape["geometry"].apply(calculate_max_distance)

# 查看结果
print(shape[["DISTRICT", "max_distance_m"]].head())

关键说明 #

• 投影选择：这里用了EPSG:32643（UTM 43N），如果你的数据覆盖印度其他区域，可以选择对应UTM带的EPSG码（比如南部用EPSG:32644）。
• MultiPolygon处理：如果区域包含飞地，会先合并所有子区域的凸包，再计算整体最长距离，避免遗漏跨飞地的最远点对。
• 效率对比：旋转卡壳算法的时间复杂度为O(n)（n为凸包顶点数），远优于暴力遍历的O(n²)，适合处理大规模数据。

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