Python检查图斑尖锐角

import geopandas as gpd
from shapely.geometry import Point
import math
import os

# 设定尖锐角的阈值为15度
sharp_angle_threshold = 15


def calculate_angle(p1, p2, p3):
    """计算三个点形成的夹角"""
    # 向量 p1p2 和 p2p3
    vector_a = (p1.x - p2.x, p1.y - p2.y)
    vector_b = (p3.x - p2.x, p3.y - p2.y)

    # 计算点积和向量的模
    dot_product = vector_a[0] * vector_b[0] + vector_a[1] * vector_b[1]
    magnitude_a = math.sqrt(vector_a[0] ** 2 + vector_a[1] ** 2)
    magnitude_b = math.sqrt(vector_b[0] ** 2 + vector_b[1] ** 2)

    # 计算角度的余弦值并取反余弦得到角度
    cos_theta = dot_product / (magnitude_a * magnitude_b)
    angle = math.degrees(math.acos(cos_theta))

    return angle


def find_sharp_angles(geom):
    """找到几何图形中的所有尖锐角"""
    sharp_points = []

    if geom.geom_type == 'Polygon':
        polygons = [geom]
    elif geom.geom_type == 'MultiPolygon':
        polygons = list(geom.geoms)
    else:
        return sharp_points  # 处理非多边形几何类型

    for polygon in polygons:
        coords = list(polygon.exterior.coords[:-1])  # 忽略最后一个重复的点

        for i in range(len(coords)):
            p2 = Point(coords[i])
            p1 = Point(coords[i - 1])
            p3 = Point(coords[(i + 1) % len(coords)])

            angle = calculate_angle(p1, p2, p3)
            print(f"顶点坐标: {coords[i]}, 计算角度: {angle}")  # 调试信息

            if angle < sharp_angle_threshold:
                print(f"发现尖锐角: {coords[i]}")  # 调试信息
                sharp_points.append(Point(coords[i]))

    return sharp_points


# 默认文件路径
input_path = r'F:\GIS\PK\New_Shapefile.shp'
output_path = r'F:\GIS\PK\demo\output_sharp_points.shp'

# 加载shapefile
gdf = gpd.read_file(input_path)

# 用于存储所有的尖锐角顶点
sharp_points = []

for _, row in gdf.iterrows():
    geom = row.geometry
    sharp_points.extend(find_sharp_angles(geom))

# 如果存在尖锐角，将它们保存为新的shapefile
if sharp_points:
    sharp_points_gdf = gpd.GeoDataFrame(geometry=sharp_points, crs=gdf.crs)

    # 创建输出目录（如果不存在）
    os.makedirs(os.path.dirname(output_path), exist_ok=True)

    # 保存为新的shapefile
    sharp_points_gdf.to_file(output_path)
    print(f"尖锐角顶点已保存到 {output_path}")
else:
    print("未检测到小于15度的尖锐角。")

规则设置为，遍历一个图斑所有折点，并计算所有三个折点形成的，两条边间的，角的角度（拆开方便理解），并将所有15度的尖锐角的位置记录并导出为点要素。

原文地址：https://blog.csdn.net/qq_40743231/article/details/142463260

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！