从DSM中提取高程信息

一、ArcMap软件设置

二、坐标提取

1、选取高程点位置

2、更新高程点(X,Y)坐标

3、更新高程点(Z)坐标

  筛选、排序，赋值非-9999高程值到Z。

三、Arcpy生成随机点并获取坐标

  软件版本：
  ArcGis Pro 3.1.5
  PY环境：
  ArcGis Pro 自带 3.9.16
  注意事项：
  1、演示代码为生成随机点并获取坐标过程，随机点是否有效需自行验证。
  2、如已存在坐标点文件，无需执行“arcpy.management.CreateRandomPoints”。

# -*- coding: cp936 -*-
import arcpy
import random
import os

dsm_paths = []
z_field = "Z"

def setValue(dsm_fields, z_field):
    with arcpy.da.UpdateCursor("point", [z_field] + dsm_fields) as cursor:
        for row in cursor:
            for i, dsm_field in enumerate(dsm_fields):
                if row[i + 1] != -9999:  # 跳过第一个字段（Z字段）
                    row[0] = row[i + 1]
                    break
            cursor.updateRow(row)

while True:
    # 输入的矢量面要素和栅格 DSM 数据
    input_path = input("输入矢量面要素路径：")
    while True:
      dsm_path = input("输入 DSM 要素路径 (输入'Q'完成 DSM 要素输入)：")
         
      if dsm_path.upper() == "Q":
          break
      dsm_paths.append(f"{dsm_path} dsm_{len(dsm_paths)}")

    # 输出点要素名称、随机点数量
    default_output_point = "point.shp"
    default_num_field = 5
    output_point = input("输入随机点要素名称，若为SHP数据，要素名称需带扩展名（默认为{}）：".format(default_output_point)) or default_output_point
    num_field = input("输入随机点要素数量（默认为{}）：".format(default_num_field)) or default_num_field
    
    # 设置工作环境
    arcpy.env.workspace = os.path.dirname(input_path)

    # 删除已存在的输出点要素类
    if arcpy.Exists(output_point):
        arcpy.management.Delete(output_point)
    
    # 删除已存在的高程点.dat
    if os.path.exists(os.path.join(os.path.dirname(input_path),"高程点.dat")):
        os.remove(os.path.join(os.path.dirname(input_path),"高程点.dat"))

    # 创建随机点要素, 更新X、Y坐标
    arcpy.management.CreateRandomPoints(os.path.dirname(input_path), output_point, input_path, "", num_field, "0.5 Meters", "POINT")
    arcpy.management.AddField(output_point, "X", "DOUBLE")
    arcpy.management.AddField(output_point, "Y", "DOUBLE")
    arcpy.management.AddField(output_point, "Z", "DOUBLE")
    arcpy.CalculateField_management(output_point, "X", "!SHAPE!.centroid.X", "PYTHON")
    arcpy.CalculateField_management(output_point, "Y", "!SHAPE!.centroid.Y", "PYTHON")

    # 提取随机点的 DSM 高程值到输出点要素类
    #arcpy.sa.ExtractValuesToPoints(output_point, dsm_path, output_point, "NONE", "VALUE_ONLY")
    arcpy.sa.ExtractMultiValuesToPoints(output_point, ";".join(dsm_paths), "NONE")
   
    # 更新点要素类的 Z 字段值
    #arcpy.management.CalculateField(output_point, "Z", "!RASTERVALU!", "PYTHON_9.3")
    dsm_fields = []
    for i in range(len(dsm_paths)):
        dsm_fields.append(f"dsm_{i}")
    setValue(dsm_fields, z_field)

    # 提取随机点的坐标和高程值
    fields = ["Y", "X", "Z"]
    data = []
    with arcpy.da.SearchCursor(output_point, fields) as cursor:
        for row in cursor:
            data.append(row)

    # 写入数据到dat文件
    with open(os.path.join(os.path.dirname(input_path),"高程点.dat"), "w") as file:
        for row in data:
            file.write(f",,{row[1]},{row[0]},{row[2]}\n")

    print("生成结束！")

    # 询问是否继续执行下一个任务
    choice = input("是否继续执行下一个任务？(Y/N) ")
    if choice.upper() != "Y":
        break

print("任务结束。")