Python实战教程一:基于Shapefile格式的不规则多边形等面积缓冲区的建立(内含配套练习数据)

import geopandas as gpdimport numpy as npimport os# 1. 配置路径input_path = r"F:\Desktop\示例数据\示例数据.shp"output_dir = r"F:\Desktop\示例数据"output_path = os.path.join(output_dir, "示例数据等积缓冲区.shp")def get_refined_buffer(geom, target_area, max_iter=5, tolerance=0.001):current_area = geom.areaif current_area <= 0:return geom# 初始估算距离 d (基于二次方程)perimeter = geom.lengthdelta_area_needed = target_area - current_aread = (-perimeter + np.sqrt(perimeter ** 2 + 4 * np.pi * delta_area_needed)) / (2 * np.pi)for i in range(max_iter):buffered_geom = geom.buffer(d)actual_area = buffered_geom.areaerror = target_area - actual_area# 如果误差在允许范围内（例如小于 1 平方米），则停止if abs(error) < tolerance:return buffered_geom# 动态调整距离：根据当前周长近似调整# 导数近似于当前缓冲后的周长current_perimeter = buffered_geom.lengthif current_perimeter == 0: breakadjustment = error / current_perimeterd += adjustmentreturn geom.buffer(d)def process_equal_area_buffer(gdf, factor=2):results = gdf.copy()new_geoms = []print("开始精确迭代计算...")for index, row in gdf.iterrows():geom = row['geometry']target = geom.area * factor# 执行迭代修正refined_geom = get_refined_buffer(geom, target)new_geoms.append(refined_geom)if index % 50 == 0 and index > 0:print(f"已处理 {index} 个城市...")results['geometry'] = new_geoms# 计算新面积 NArea (单位：平方千米)，并保留 6 位小数# 注意：这里计算的是最终几何体的实际面积results['NArea'] = (results['geometry'].area / 1_000_000).round(6)return resultstry:# 读取数据city_gdf = gpd.read_file(input_path)# 执行处理final_gdf = process_equal_area_buffer(city_gdf, 2)# 设置输出 Schema，防止 NArea 写入失败input_schema = gpd.io.file.infer_schema(city_gdf)output_schema = input_schema.copy()output_schema['properties']['NArea'] = 'float:24.6'# 保存结果print("正在写入文件，请稍候...")final_gdf.to_file(output_path,driver='ESRI Shapefile',encoding='utf-8',schema=output_schema,engine='fiona')print(f"✅ 处理成功！结果已保存至: {output_path}")print(f"统计：计算得到的 NArea 均值约为 {final_gdf['NArea'].mean():.4f} km²")except Exception as e:print(f"❌ 运行出错: {e}")