我正在尝试使用Python创建一个多边形特征,它表示光栅图像中值超过一定量的所有单元格。问题是我的输出特征最终出现在所谓的“空岛”(位于大西洋),这是坐标系出现问题的迹象。这是我最后使用的代码示例:
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
# Path to the output shapefile (replace 'path\to\your\output.gdb\output_polygons' with the actual path)
output_shapefile = "set a name for the output file"
# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")
# Convert raster to NumPy array
raster_array = arcpy.RasterToNumPyArray(input_raster)
# Create a list to store polygon geometries
polygon_geometries = []
# Get the cell size of the raster
cell_size = arcpy.Raster(input_raster).meanCellWidth
# Get the spatial reference of the raster
raster_spatial_ref = arcpy.Describe(input_raster).spatialReference
# Iterate through each cell in the raster
rows, cols = np.where(raster_array > threshold_value)
for row, col in zip(rows, cols):
x = col * cell_size + 0.5 * cell_size
y = row * 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))
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 polygon in polygon_geometries:
cursor.insertRow([polygon, threshold_value])
print("Polygon creation completed.")
到目前为止,我检查了每个图层的属性:地理数据库、数据集、输入要素、输出要素(代码运行后)以及我正在处理的地图选项卡,它们的所有坐标系都是相同的。为什么输出层无法获得正确的坐标对我来说没有意义。我使用 ArcGIS Pro 作为我的主要 GIS 软件。