在ArcPy中围绕点创建矩形时如何正确定义单位？

我想在WGS 1984投影（单位为十进制度数）的地图上，围绕点创建宽50米、高500米的矩形。但用以下代码运行后，矩形尺寸异常巨大——显然是直接用米作为单位，和十进制度数不匹配导致的。

后来我计算了点东侧50米、南侧500米对应的经纬度差值作为宽高，代码执行后输出的要素类有属性表，但没有要素显示。之前手动设置w=0.1、h=1时能显示，但关闭再打开图层后又看不到了。

inFeatures = "potentialpoint.shp"
outFeatureClass = "rectanglestest.shp"
rectangleWidth= 50
rectangleHeight = 500

# create a new feature class to hold the rectangles
arcpy.CreateFeatureclass_management("D:\OceanWind2023\DepletionExperiments", outFeatureClass, "POLYGON")

# loop through each point and create a rectangle around it
with arcpy.da.SearchCursor(inFeatures, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        x, y = row[0]
        array = arcpy.Array([arcpy.Point(x - rectangleWidth / 2.0, y - rectangleHeight / 2.0),
                             arcpy.Point(x - rectangleWidth / 2.0, y + rectangleHeight / 2.0),
                             arcpy.Point(x + rectangleWidth / 2.0, y + rectangleHeight / 2.0),
                             arcpy.Point(x + rectangleWidth / 2.0, y - rectangleHeight / 2.0)])
        polygon = arcpy.Polygon(array)
        with arcpy.da.InsertCursor(outFeatureClass, ["SHAPE@"]) as insertCursor:
            insertCursor.insertRow([polygon])

1. 核心问题拆解 #

• 单位不匹配：WGS 1984是地理坐标系，用的是经纬度（十进制度数），1度经度在赤道附近≈111公里，直接把50米当度数用，矩形自然大到离谱。
• 空间参考未定义：创建输出要素类时没指定空间参考，导致要素和地图的空间参考不匹配，就算插入成功也可能显示不出来。
• 代码效率问题：把InsertCursor放在循环里，每次循环都创建一次，不仅慢还容易出问题。

2. 最优解决方法：先转平面坐标系计算偏移 #

地理坐标系下经纬度的米数偏移随纬度变化，直接计算误差大。建议先转UTM（通用横轴墨卡托）这类平面坐标系（单位是米），计算完矩形再转回WGS 1984：

import arcpy

inFeatures = "potentialpoint.shp"
outFeatureClass = "rectanglestest.shp"
rectangleWidth = 50  # 单位：米
rectangleHeight = 500  # 单位：米

# 获取输入点的空间参考（WGS 1984）
input_sr = arcpy.Describe(inFeatures).spatialReference

# 创建输出要素类，必须指定空间参考，避免后续显示问题
output_path = r"D:\OceanWind2023\DepletionExperiments"
arcpy.CreateFeatureclass_management(output_path, outFeatureClass, "POLYGON", spatial_reference=input_sr)

# 自动匹配UTM分区（arcpy会根据点的纬度选对应的UTM投影）
utm_sr = input_sr.exportToUTM()

# 把InsertCursor移到循环外，提升效率
with arcpy.da.InsertCursor(outFeatureClass, ["SHAPE@"]) as insert_cursor:
    with arcpy.da.SearchCursor(inFeatures, ["SHAPE@"]) as search_cursor:
        for row in search_cursor:
            point = row[0]
            # 将点从WGS 1984转成UTM坐标系
            point_utm = point.projectAs(utm_sr)
            x, y = point_utm.centroid.X, point_utm.centroid.Y
            
            # 计算矩形四个顶点的UTM坐标
            half_w = rectangleWidth / 2.0
            half_h = rectangleHeight / 2.0
            utm_points = arcpy.Array([
                arcpy.Point(x - half_w, y - half_h),
                arcpy.Point(x - half_w, y + half_h),
                arcpy.Point(x + half_w, y + half_h),
                arcpy.Point(x + half_w, y - half_h)
            ])
            polygon_utm = arcpy.Polygon(utm_points, utm_sr)
            
            # 把矩形转回WGS 1984地理坐标系
            polygon_wgs = polygon_utm.projectAs(input_sr)
            
            # 插入到输出要素类
            insert_cursor.insertRow([polygon_wgs])

3. 若一定要直接用经纬度偏移（精度略低） #

如果不想转投影，可以根据当前纬度计算米和度数的转换系数：

import math

# 假设当前点纬度为y（十进制度数），替换成你的点实际纬度
y = 30.0  
# 1度经度对应的米数：111320 * cos(纬度弧度)
lon_per_m = 1 / (111320 * math.cos(y * math.pi / 180))
# 1度纬度对应的米数≈110574米（全球大致范围）
lat_per_m = 1 / 110574

# 转换50米、500米为度数
width_deg = 50 * lon_per_m
height_deg = 500 * lat_per_m

注意：这种方法在高纬度地区误差会变大，优先用投影转换的方法。

4. 要素不显示的排查步骤 #

1. 检查空间参考：右键输出图层→数据→投影和变换→投影，确保和地图的空间参考一致。
2. 缩放至图层范围：右键图层→缩放至图层，可能要素太小或位置不对，地图没显示到。
3. 检查属性表：打开要素类属性表，看SHAPE字段是否有值，为空则是代码插入失败；有值则大概率是空间参考问题。

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