Python狭长型图斑检测

import os
import sys
import geopandas as gpd
from shapely.geometry import LineString# 设置默认路径
default_path = r"D:\..."
os.chdir(default_path)# 读取矢量图层
vector_layer = gpd.read_file('狭长形图斑.shp')# 定义新的计算最小宽度的函数
def new_min_width(geometry):if geometry.geom_type == 'Polygon':coords = list(geometry.exterior.coords)if len(coords) <= 5:# 如果折点数量少于等于 5 个，直接计算平均宽度bounds = geometry.boundsdx = bounds[2] - bounds[0]dy = bounds[3] - bounds[1]return (dx + dy) / 2else:# 对于不规则多边形，继续使用原来的分割方法min_width_value = float('inf')min_x = float('inf')max_x = float('-inf')min_y = float('inf')max_y = float('-inf')for x, y in coords:if x < min_x:min_x = xif x > max_x:max_x = xif y < min_y:min_y = yif y > max_y:max_y = ydx = max_x - min_xdy = max_y - min_yfor step in range(0, 10):x1 = min_x + step * dx / 10x2 = min_x + (step + 1) * dx / 10y1 = min_y + step * dy / 10y2 = min_y + (step + 1) * dy / 10points_on_line1 = geometry.intersection(gpd.GeoSeries([LineString([(x1, min_y), (x1, max_y)])]))points_on_line2 = geometry.intersection(gpd.GeoSeries([LineString([(x2, min_y), (x2, max_y)])]))for p1 in points_on_line1:for p2 in points_on_line2:if p1 and p2:distance = p1.distance(p2)if distance < min_width_value:min_width_value = distancereturn min_width_valueelif geometry.geom_type == 'MultiPolygon':min_width_value = float('inf')for polygon in geometry.geoms:coords = list(polygon.exterior.coords)if len(coords) <= 5:# 如果折点数量少于等于 5 个，直接计算平均宽度bounds = polygon.boundsdx = bounds[2] - bounds[0]dy = bounds[3] - bounds[1]temp_width = (dx + dy) / 2if temp_width < min_width_value:min_width_value = temp_widthelse:# 对于不规则多边形，继续使用原来的分割方法min_x = float('inf')max_x = float('-inf')min_y = float('inf')max_y = float('-inf')for x, y in coords:if x < min_x:min_x = xif x > max_x:max_x = xif y < min_y:min_y = yif y > max_y:max_y = ydx = max_x - min_xdy = max_y - min_yfor step in range(0, 10):x1 = min_x + step * dx / 10x2 = min_x + (step + 1) * dx / 10y1 = min_y + step * dy / 10y2 = min_y + (step + 1) * dy / 10points_on_line1 = polygon.intersection(gpd.GeoSeries([LineString([(x1, min_y), (x1, max_y)])]))points_on_line2 = polygon.intersection(gpd.GeoSeries([LineString([(x2, min_y), (x2, max_y)])]))for p1 in points_on_line1:for p2 in points_on_line2:if p1 and p2:distance = p1.distance(p2)if distance < min_width_value:min_width_value = distancereturn min_width_valueprint("正在计算图斑最小宽度...")
# 检测每个图斑的最小宽度小于 5m，并筛选出这些图斑
narrow_polygons = vector_layer[vector_layer.geometry.apply(new_min_width) < 5]
print("图斑最小宽度计算完成，正在计算中心点...")# 计算筛选出的图斑的中心点并输出为新的矢量图层
center_points = narrow_polygons.centroid
center_points.to_file('最小宽度检查1.shp')
print("中心点计算完成，结果已保存为'最小宽度检查1.shp'。")

规则设置为，遍历一个图斑的所有折点，数量少于等于5个的直接计算平均宽度是否小于5m，其他的图斑按照分割方法计算近似最小宽度：