如何设置从光栅图像中提取的多边形的正确坐标?

问题描述 投票:0回答:0

我正在尝试使用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 软件。

python gis raster analysis
© www.soinside.com 2019 - 2024. All rights reserved.