ArcGIS Pro 3.3 Notebook脚本优化：多边形贴合折线转角路径

现有脚本可生成以点为中心、与折线角度对齐的多边形，但在折线转角处，多边形无法贴合路径，需实现多边形随转角适配的效果。当前脚本可运行但未达预期，代码如下：

# Name:        Create Polygons at Points
# Purpose:     Create a polygon centered at each Structure and aligned with the angle of the polyline
# Created:     28/08/2024

import arcpy
import math

# Constants
WORKSPACE = r"PATH TO... Map.gdb"
STRUCTURES_LAYER_NAME = "Structures"
ALIGNMENT_LAYER_NAME = "Route"
OUTPUT_FC_NAME = "Aligned_Polygons"
TOLERANCE = 0.0125  # Updated tolerance

# Set workspace
arcpy.env.workspace = WORKSPACE

def get_map_by_name(aprx, map_name):
    """Retrieve a map by its name."""
    maps = aprx.listMaps(map_name)
    if not maps:
        raise RuntimeError(f"No map named '{map_name}' found in the current project.")
    return maps[0]

def get_layer_by_name(map_obj, layer_name):
    """Retrieve a layer by its name."""
    layers = map_obj.listLayers(layer_name)
    if not layers:
        raise RuntimeError(f"No layer named '{layer_name}' found in the map.")
    return layers[0]

def set_layer_visible_scale(layer):
    """Set the layer to be visible at all scales."""
    layer.minScale = 0
    layer.maxScale = 0
    print(f"Set layer '{layer.name}' to be visible at all scales.")

def create_output_feature_class(output_fc, spatial_ref):
    """Create the output feature class."""
    if arcpy.Exists(output_fc):
        arcpy.management.Delete(output_fc)
        print(f"Deleted existing feature class '{output_fc}'.")
    arcpy.management.CreateFeatureclass(arcpy.env.workspace, output_fc, "POLYGON", spatial_reference=spatial_ref)
    print(f"Output feature class '{output_fc}' created successfully.")

def add_fields_if_not_exist(output_fc, fields):
    """Add fields to the output feature class if they do not exist."""
    existing_fields = {f.name for f in arcpy.ListFields(output_fc)}
    for field_name, field_type in fields.items():
        if field_name not in existing_fields:
            arcpy.management.AddField(output_fc, field_name, field_type)
    print(f"Fields added to '{output_fc}'.")

def calculate_rotated_rectangle(center, width, height, angle):
    """Calculate a rotated rectangle's vertices."""
    angle_rad = math.radians(angle)
    dx = width / 2
    dy = height / 2
    corners = [(-dx, -dy), (dx, -dy), (dx, dy), (-dx, dy)]
    return [(center[0] + x * math.cos(angle_rad) - y * math.sin(angle_rad),
             center[1] + x * math.sin(angle_rad) + y * math.cos(angle_rad)) for x, y in corners]

def is_point_near_vertex(point, alignment_fc, tolerance=TOLERANCE):
    """Check if a point is near any vertex of the alignment features."""
    point_geom = arcpy.PointGeometry(arcpy.Point(*point))
    with arcpy.da.SearchCursor(alignment_fc, ["SHAPE@"]) as cursor:
        for row in cursor:
            geometry = row[0]
            if geometry:
                vertices = [arcpy.PointGeometry(pnt) for pnt in geometry.getPart(0)]
                if any(vertex.distanceTo(point_geom) <= tolerance for vertex in vertices):
                    return True
    return False

def main():
    # Prompt user for map name
    map_name = input("Enter the map name to use: ")
    # Prompt user for rectangle dimensions
    rectangle_length = float(input("Enter the rectangle length (in feet): "))
    rectangle_width = float(input("Enter the rectangle width (in feet): "))
    # Load the ArcGIS Project
    aprx = arcpy.mp.ArcGISProject("CURRENT")
    # Get the map and layers
    map_to_use = get_map_by_name(aprx, map_name)
    structures_layer = get_layer_by_name(map_to_use, STRUCTURES_LAYER_NAME)
    alignment_layer = get_layer_by_name(map_to_use, ALIGNMENT_LAYER_NAME)
    # Set layers to be visible at all scales
    set_layer_visible_scale(structures_layer)
    set_layer_visible_scale(alignment_layer)
    # Get feature class paths
    structures_fc = structures_layer.dataSource
    alignment_fc = alignment_layer.dataSource
    # Extract spatial reference
    spatial_ref = arcpy.Describe(structures_fc).spatialReference
    
    # Create output feature class
    output_fc = OUTPUT_FC_NAME
    create_output_feature_class(output_fc, spatial_ref)
    
    # Define fields to add
    fields_to_add = {
        "Length_ft": "DOUBLE",
        "Width_ft": "DOUBLE",
        "Unique_ID": "TEXT",
        "Type": "TEXT",
        "Feature": "TEXT",
        "Height": "DOUBLE",
        "County": "TEXT"
    }
    add_fields_if_not_exist(output_fc, fields_to_add)
    
    # Collect points data
    points_data = []
    with arcpy.da.SearchCursor(structures_fc, ["OBJECTID", "SHAPE@XY", "Unique_ID", "Type", "Feature", "Height", "County"], sql_clause=(None, 'ORDER BY OBJECTID')) as structures_cursor:
        for structure in structures_cursor:
            current_point = structure[1]
            if is_point_near_vertex(current_point, alignment_fc):
                points_data.append(structure[1:])  # Collect all attributes except OBJECTID
    
    # Create polygons
    with arcpy.da.InsertCursor(output_fc, ["SHAPE@", "Length_ft", "Width_ft", "Unique_ID", "Type", "Feature", "Height", "County"]) as insert_cursor:
        for i, data in enumerate(points_data):
            current_point, unique_id, type_field, feature, height, county = data
            angle = 0  # Default angle for the first point
            
            # Calculate angle for the first point
            if i == 0 and len(points_data) > 1:
                second_point = points_data[1][0]  # Get the second point's coordinates
                dx = second_point[0] - current_point[0]
                dy = second_point[1] - current_point[1]
                angle = math.degrees(math.atan2(dy, dx))
            
            # Calculate angle for all other points
            elif i > 0:
                previous_point = points_data[i - 1][0]  # Get the previous point's coordinates
                dx = current_point[0] - previous_point[0]
                dy = current_point[1] - previous_point[1]
                angle = math.degrees(math.atan2(dy, dx))
            
            # Calculate rectangle points centered on the current point
            rectangle = calculate_rotated_rectangle(current_point, rectangle_length, rectangle_width, angle)
            polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in rectangle]), spatial_ref)
            insert_cursor.insertRow([polygon, rectangle_length, rectangle_width, unique_id, type_field, feature, height, county])
            print(f"Processed {i + 1} features...")
    
    print(f"Polygons created and saved in {output_fc}. Total features processed: {len(points_data)}.")

if __name__ == "__main__":
    main()

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