ArcPy提取栅格高值区域多边形坐标异常至Null Island问题排查

你的代码核心问题是未考虑栅格的实际地理原点（左下角坐标），仅通过行列号和像元大小计算坐标，导致所有多边形从(0,0)位置生成，最终落在Null Island。此外，ArcGIS栅格的行号从左上角开始计数，直接用row计算Y坐标会导致上下颠倒，也需修正。

关键修正点： #

• 获取栅格的extent，提取左下角的X（x_min）、Y（y_min）及栅格总行数（n_rows）
• 修正X坐标：基于左下角X原点计算实际地理坐标
• 修正Y坐标：将左上角起始的行号转换为左下角起始的计数逻辑

修正后的完整代码： #

import arcpy
import numpy as np

# Set the ArcGIS Pro project path
arcpy.env.workspace = r"Set the GDB"

# Path to the input TIFF raster 
input_raster = r"path to tiff file"

# Threshold value to create polygons (change as needed)
threshold_value = 40

output_shapefile = "set a name for the output file"

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

# Convert raster to NumPy array
raster_obj = arcpy.Raster(input_raster)
raster_array = arcpy.RasterToNumPyArray(raster_obj)

# Create a list to store polygon geometries
polygon_geometries = []

# Get raster properties
cell_size = raster_obj.meanCellWidth
raster_spatial_ref = arcpy.Describe(input_raster).spatialReference
raster_extent = raster_obj.extent
x_min = raster_extent.XMin
y_min = raster_extent.YMin
n_rows = raster_obj.height  # 栅格总行数

# Iterate through each cell in the raster
rows, cols = np.where(raster_array > threshold_value)
for row, col in zip(rows, cols):
    # 修正X坐标：基于左下角X原点计算
    x = x_min + col * cell_size + 0.5 * cell_size
    # 修正Y坐标：转换行号为从左下角计数
    y = y_min + (n_rows - row - 1) * cell_size + 0.5 * cell_size
    
    vertex_list = [
        arcpy.Point(x - 0.5 * cell_size, y - 0.5 * cell_size),
        arcpy.Point(x - 0.5 * cell_size, y + 0.5 * cell_size),
        arcpy.Point(x + 0.5 * cell_size, y + 0.5 * cell_size),
        arcpy.Point(x + 0.5 * cell_size, y - 0.5 * cell_size),
        arcpy.Point(x - 0.5 * cell_size, y - 0.5 * cell_size)
    ]
    polygon = arcpy.Polygon(arcpy.Array(vertex_list), raster_spatial_ref)
    polygon_geometries.append(polygon)

# Create the output shapefile
output_shapefile = arcpy.management.CreateFeatureclass(
    arcpy.env.workspace,
    output_shapefile,
    "POLYGON",
    spatial_reference=raster_spatial_ref
)

# Add a field to store the values of each cell (optional)
arcpy.management.AddField(output_shapefile, "CellValue", "DOUBLE")

# Insert the polygons into the feature class
with arcpy.da.InsertCursor(output_shapefile, ["SHAPE@", "CellValue"]) as cursor:
    for idx, polygon in enumerate(polygon_geometries):
        # 若需存储单元格实际值而非阈值，替换为 raster_array[rows[idx], cols[idx]]
        cursor.insertRow([polygon, raster_array[rows[idx], cols[idx]]])

print("Polygon creation completed.")

额外说明： #

1. 提前将栅格对象赋值给raster_obj，避免重复调用arcpy.Raster()提升效率
2. 创建Polygon对象时直接传入空间参考，确保几何对象坐标系与栅格一致
3. 可选择存储单元格实际值，替代固定阈值

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