地块是一些超复杂的多边形，其中一些多边形是简单多边形，有些多边形是包含内部空洞的复杂多边形，而有些多边形是两三个复杂多边形（可能包含内部空洞）的“复合体”。
具有内部空洞的复杂多边形（边界高亮）
两三个复杂多边形（可能包含内部空洞）的“复合体”（边界高亮）
读取地块数据并使用Voronoi图方法顺序计算shapefile文件中每一个复杂多边形最大内切圆的圆心坐标和半径，并添加新字段将圆心坐标和半径写入shapefile文件中
1.导入相关库
# 导入NumPy库，用于数值计算
import numpy as np
# 导入Polygon、MultiPolygon和Point类，用于处理几何对象
from shapely.geometry import Polygon, MultiPolygon, Point
# 导入Voronoi类，用于生成Voronoi图
from scipy.spatial import Voronoi
# 导入geopandas库，用于处理地理空间数据
import geopandas as gpd
# 导入matplotlib.pyplot用于绘图
import matplotlib.pyplot as plt
# 从pyproj模块中导入Transformer类,Transformer类则用于坐标变换
from pyproj import Transformer2.实现基于多边形特性的自适应采样策略
包括对于没有空洞的多边形采用较少的采样点，而对于有空洞的多边形使用更多的采样点，同时根据边界的长度调整采样点的密度
def sample_points_on_edge(geometry, num_points=300):
# 检查几何对象是否有效，如果无效则尝试修复
if not geometry.is_valid:
geometry = geometry.buffer(0)
sampled_points = [] # 初始化采样点列表
# 判断几何对象是否为多边形
if isinstance(geometry, Polygon):
# 获取多边形的外边界和所有内边界（空洞）
exteriors = [geometry.exterior] + list(geometry.interiors)
# 根据是否有内边界（空洞）调整采样点总数
num_points_adjusted = num_points if len(geometry.interiors) > 0 else int(num_points * 0.75)
# 计算所有边界的总长度
total_length = sum([ext.length for ext in exteriors])
# 对每个边界进行处理
for ext in exteriors:
length = ext.length # 计算边界的长度
# 根据边界长度和调整后的采样点数计算此边界的采样点数
# 根据每个边界（外边界和内边界）的长度比例分配采样点数，以实现在边界较长的部分采用更密集的采样，在边界较短的部分采用较稀疏的采样
edge_num_points = int(num_points_adjusted * (length / total_length))
# 在边界上根据计算出的采样间隔插值采样点
distances = np.linspace(0, length, edge_num_points)
sampled_points += [ext.interpolate(distance) for distance in distances]
return sampled_points3.计算多边形的最大内切圆的圆心坐标
def find_largest_inscribed_circle_center(polygon, num_points=300):
# 对多边形的外边界和所有内边界（空洞）上的点进行采样
sampled_points = sample_points_on_edge(polygon, num_points)
# 移除采样点中的重复点，保留唯一点集
points = np.unique(np.array([[point.x, point.y] for point in sampled_points]), axis=0)
# 打印点集，用于调试(看看具体是哪些点)
# print(points)
# 绘制多边形和采样点（可视化出来观察看看采样点是密集还是稀疏）
plt.figure(figsize=(8, 8))
for exterior in [polygon.exterior] + list(polygon.interiors):
plt.plot(*exterior.xy, color='blue') # 绘制多边形边界
plt.scatter(points[:, 0], points[:, 1], color='red') # 绘制采样点
# 确保点集中有足够的唯一点来构建Voronoi图（至少需要3个不同的点）
if len(points) < 3:
raise ValueError("Not enough unique points to construct a Voronoi diagram.")
# 使用点集生成Voronoi图
vor = Voronoi(points)
# 绘制Voronoi图
voronoi_plot_2d(vor, show_vertices=False, line_colors='orange', line_width=2, point_size=2)
plt.title("Polygon and Sampled Points with Voronoi Diagram")
# 过滤掉位于多边形外部和内部空洞中的Voronoi顶点
valid_vertices = [Point(vertex) for vertex in vor.vertices if polygon.contains(Point(vertex))]
# 初始化最大半径和对应的圆心变量
largest_radius = 0
largest_center = None
# 遍历所有有效的Voronoi顶点
for vertex in valid_vertices:
# 计算顶点到多边形边界的最短距离
distance = polygon.boundary.distance(vertex)
# 更新最大半径和对应的圆心
if distance > largest_radius:
largest_radius = distance
largest_center = vertex
return largest_center, largest_radius4.定义坐标转换函数
从CGCS2000_3_Degree_GK_CM_108E转换为WGS84。如果自己不需要可以注释或者删掉
def transform_to_wgs84(x, y):
# 初始化转换器，从CGCS2000（EPSG:4545）到WGS84
transformer = Transformer.from_crs("epsg:4545", "epsg:4326", always_xy=True)
# 执行转换
x2, y2 = transformer.transform(x, y)
return x2, y25.处理shp文件并更新最大内切圆的圆心坐标
def process_shp_file(shp_path, updated_shp_path):
# 使用gpd.read_file(shp_path)读取Shapefile文件并将其转换为GeoDataFrame对象，这使得后续对地理数据的处理更加方便
gdf = gpd.read_file(shp_path)
# 在这里我只是想统计shp文件中没有处理好的地块（即复合地块）的数量以及记录和输出其ID
multi_polygon_count = 0 # 记录MultiPolygon的数量
multi_polygon_ids = [] # 记录MultiPolygon的ID
centers = [] # 初始化圆心坐标列表
radii = [] # 初始化半径列表
# 遍历shp文件中的每个几何对象
for index, geom in enumerate(gdf.geometry):
if isinstance(geom, MultiPolygon):
multi_polygon_count += 1 # 更新MultiPolygon的数量
multi_polygon_ids.append(index) # 记录MultiPolygon的ID
# 初始化临时存储最大圆的变量
max_radius = 0
max_center = None
# 对MultiPolygon中的每个Polygon进行处理
for polygon in geom.geoms:
# 根据是否有空洞选择采样点数量
num_points = 300 if len(polygon.interiors) > 0 else int(300 * 0.75)
center, radius = find_largest_inscribed_circle_center(polygon, num_points)
if radius > max_radius:
max_radius = radius
max_center = center
# 在找到最大内切圆的圆心后，立即转换坐标（可将其写成一个函数调用）
if max_center is not None:
# 转换坐标
center_x_transformed, center_y_transformed = transform_to_wgs84(max_center.x, max_center.y)
# 更新圆心为转换后的坐标
max_center = Point(center_x_transformed, center_y_transformed)
centers.append(max_center)
radii.append(max_radius)
elif isinstance(geom, Polygon):
# 如果多边形有内部空洞，则采样点数量为300，否则为225，以适应多边形的特性
num_points = 300 if len(geom.interiors) > 0 else int(300 * 0.75)
# num_points = 300 if len(geom.interiors) > 0 else 225
center, radius = find_largest_inscribed_circle_center(geom, num_points)
# 在找到最大内切圆的圆心后，立即转换坐标
if center is not None:
# 转换坐标
center_x_transformed, center_y_transformed = transform_to_wgs84(center.x, center.y)
# 更新圆心为转换后的坐标
center = Point(center_x_transformed, center_y_transformed)
centers.append(center)
radii.append(radius)
else:
# 对于非Polygon/MultiPolygon类型，简单地为这些几何对象分配了空的圆心和0半径
centers.append(None)
radii.append(0)
# 将计算得到的圆心坐标和半径添加到GeoDataFrame的新字段中,确保新字段与GeoDataFrame的行数匹配
gdf['center_x'] = [center.x if center else None for center in centers]
gdf['center_y'] = [center.y if center else None for center in centers]
gdf['radius'] = radii
# 使用gdf.to_file(updated_shp_path)将更新后的GeoDataFrame保存为新的Shapefile文件
# 指定编码为UTF-8以正确处理中文（GBK完全不行，Big5有些中文不显示）
gdf.to_file(updated_shp_path, encoding='UTF-8')
# 返回multipolygon的数量和ID
return multi_polygon_count, multi_polygon_ids6.使用matplotlib绘制具有内部空洞的复杂多边形
用于绘制具有外边界和内部空洞的多边形。此函数首先绘制多边形的外边界，然后遍历并绘制每个内部空洞的边界。内部空洞使用虚线样式以便于区分
def draw_polygon_with_holes(geom, ax):
# 通过geom.exterior.xy获取外边界的x和y坐标，使用ax.plot(x, y, 'b')在给定的轴(ax)上绘制外边界
x, y = geom.exterior.xy
# ax.fill(x, y, alpha=0.5, fc='lightblue', ec='black')
ax.plot(x, y, 'b')
# 通过遍历geom.interiors获取每个内部空洞的坐标，并使用ax.plot(x, y, 'b', linestyle='--')以虚线样式绘制内部空洞的边界
for interior in geom.interiors:
x, y = interior.xy
ax.plot(x, y, 'b', linestyle='--') # 使用虚线样式绘制内部空洞的边界7.绘制最大内切圆，并在matplotlib图形界面中标记圆心
使用plt.Circle创建一个圆形对象，其中心位于参数center指定的位置，半径等于参数radius指定的值。然后，使用ax.add_artist(circle)将这个圆形添加到图形界面的指定轴上
def draw_largest_inscribed_circle(ax, center, radius):
circle = plt.Circle((center.x, center.y), radius, color='r', fill=False)
ax.add_artist(circle)
# 在圆心位置绘制一个红色的圆点
ax.plot(center.x, center.y, 'ro') # 红色圆点表示圆心8.处理和绘制GeoDataFrame中指定范围内的几何对象
在这里是想绘制shp文件中指定范围内的一些复杂多边形，将其可视化出来看看结果（最大内切圆和圆心）是否精确
def process_and_draw_polygons(gdf):
# 通过enumerate(gdf.geometry[29:40])遍历GeoDataFrame中的第30至40个几何对象（索引从0开始）
for index, geom in enumerate(gdf.geometry[29:40]): # 仅处理第30至40个几何对象
fig, ax = plt.subplots()
# 对于每个几何对象，根据其类型（MultiPolygon或Polygon），使用draw_polygon_with_holes函数绘制多边形的外边界和内部空洞，
# 然后使用find_largest_inscribed_circle_center函数计算最大内切圆的圆心和半径，并使用draw_largest_inscribed_circle函数进行绘制
if isinstance(geom, MultiPolygon):
for polygon in geom.geoms:
draw_polygon_with_holes(polygon, ax)
center, radius = find_largest_inscribed_circle_center(polygon)
draw_largest_inscribed_circle(ax, center, radius)
elif isinstance(geom, Polygon):
draw_polygon_with_holes(geom, ax)
center, radius = find_largest_inscribed_circle_center(geom)
draw_largest_inscribed_circle(ax, center, radius)
# 为每个几何对象创建的图形设置等比例的视图，并设置标题以表示当前处理的多边形序号和其最大内切圆。最后，使用plt.show()显示绘制的图形
ax.set_aspect('equal')
plt.title(f'Polygon {index + 10} with Largest Inscribed Circle')
plt.show()9.调用示例，查看效果
# 调用示例，路径需要根据实际情况进行替换
shp_path = 'shp/cut.shp' # 输入shp文件的路径
updated_shp_path = 'new_shp/cut.shp' # 输出更新后shp文件的路径
gdf = gpd.read_file(shp_path)
multi_polygon_count, multi_polygon_ids = process_shp_file(shp_path, updated_shp_path)
process_and_draw_polygons(gdf)
# 输出MultiPolygon的数量和对应的ID
print("MultiPolygon count:", multi_polygon_count)
print("MultiPolygon IDs:", multi_polygon_ids)因为MultiPolygon的存在，代码复杂了许多。同时，在某些地方在处理完之后需要可视化出来观察效果，不得已增加了一些代码。整个代码运行起来时没有问题的，下面是运行结果：
（这里中间的多边形的内切圆因为圆心画得太大所有看着不像）
代码有点多，也有点复杂，可能有不足或错误的地方。如果有什么问题可以在评论区留言
扫一扫
837
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
