【无标题】

如何用arcgis制作100米的地球化学采样框

#导入点，量一下距离，然后整除计算，因为奥维上的距离和arcgis上的距离有区别，因此需要重新量一下，然后再计算
import arcpy
import os

设置工作空间

arcpy.env.workspace = r"C:\Users\Desktop\345" # 请确保路径没有中文字符
arcpy.env.overwriteOutput = True

输入点文件路径

point_layer = r"C:\Users\345\TYdian1.shp" # 更改为英文路径

输出正方形文件路径

output_square_layer = r"C:\Users\Desktop\345\squares.shp" # 更改为英文路径

设置正方形的边长（单位：米）

square_length = 127.44737353 # 100米

旋转角度（度）

rotation_angle = 51.5 # 旋转 50 度

检查输入点文件是否存在

if arcpy.Exists(point_layer):
print(“点图层 {} 存在。”.format(point_layer))
else:
print(“点图层 {} 不存在。”.format(point_layer))
exit()

获取点图层的空间参考

point_sr = arcpy.Describe(point_layer).spatialReference

如果点图层没有空间参考，使用 WGS 84 默认空间参考

if point_sr is None:
print(“点图层没有空间参考，使用默认空间参考 WGS 84。”)
point_sr = arcpy.SpatialReference(4490)

如果输出正方形文件已存在，删除它

if os.path.exists(output_square_layer):
print(“文件 {} 已经存在，尝试删除它。”.format(output_square_layer))
os.remove(output_square_layer)

创建新的空图层来存储正方形（使用与点图层相同的空间参考）

try:
arcpy.CreateFeatureclass_management(arcpy.env.workspace, “squares.shp”, “POLYGON”, spatial_reference=point_sr)
print(“Feature class 创建成功！”)
except arcpy.ExecuteError:
print(“创建 Feature class 失败: {}”.format(arcpy.GetMessages()))
exit()

获取点要素

point_cursor = arcpy.da.SearchCursor(point_layer, [“SHAPE@”])

遍历每个点，生成正方形

for point in point_cursor:
point_geom = point[0]

# 获取点的坐标
x, y = point_geom.centroid.X, point_geom.centroid.Y
print("点坐标: ({}, {})".format(x, y))

# 计算正方形的四个角的偏移量（正方形边长为100米，偏移量为50米）
delta = square_length / 2  # 正方形半边长 50米

# 计算四个角的坐标
points = [
    arcpy.Point(x - delta, y - delta),  # 左下
    arcpy.Point(x + delta, y - delta),  # 右下
    arcpy.Point(x + delta, y + delta),  # 右上
    arcpy.Point(x - delta, y + delta)   # 左上
]

 # 创建正方形的多边形
square_polygon = arcpy.Polygon(arcpy.Array(points))

# 旋转正方形 50 度，基于正方形的中心
center = square_polygon.centroid

# 旋转矩阵的计算（绕中心点旋转）
angle_rad = math.radians(rotation_angle)  # 转换为弧度
cos_angle = math.cos(angle_rad)
sin_angle = math.sin(angle_rad)

# 创建仿射变换（包括平移和旋转）
def rotate_point(x, y, cx, cy, angle_rad):
    """将点绕中心点旋转指定角度"""
    dx = x - cx
    dy = y - cy
    new_x = cx + dx * cos_angle - dy * sin_angle
    new_y = cy + dx * sin_angle + dy * cos_angle
    return arcpy.Point(new_x, new_y)

# 旋转正方形的四个角
rotated_points = [rotate_point(p.X, p.Y, center.X, center.Y, angle_rad) for p in points]

# 创建旋转后的正方形多边形
rotated_square = arcpy.Polygon(arcpy.Array(rotated_points))

# 插入正方形到输出图层
try:
    with arcpy.da.InsertCursor(output_square_layer, ["SHAPE@"]) as cursor:
        cursor.insertRow([rotated_square])
    print("正方形生成成功，中心点坐标: ({}, {})".format(x, y))
except arcpy.ExecuteError:
    print("插入正方形失败: {}".format(arcpy.GetMessages()))

print(“所有正方形生成完毕！”)