在地理信息系统（GIS）中，图幅号是用于标识地理区域的一种重要编号，尤其在地图制图和数据管理中具有重要作用。中国的国标 GB/T 13989-2012 《国家基本比例尺地形图分幅和编号》提供了一套完整的图幅编号标准。在 ArcGIS 环境中，我们可以通过编写 Python 脚本来实现自动计算要素的图幅号。本文将详细介绍如何在 ArcGIS 中计算要素涉及的图幅号。
图幅号的计算标准根据 GB/T 13989-2012《国家基本比例尺地形图分幅和编号》规定，涉及到地图制图和数据管理的规范。主要目的是通过规定的规则和比例尺将地理区域划分成标准的图幅，从而实现数据的一致性和管理的便利性。
图幅号的层级分为不同的比例尺，例如：1:1000000、1:500000、1:250000、1:100000 等。每个层级的图幅编号方式不同，随着比例尺的缩小，图幅号的细分程度和编号规则也会相应变化。
大图幅号（基础图幅号）：通常由字母和数字组合表示，例如“B123”。
小图幅号（详细图幅号）：在大图幅号的基础上，进一步细分，通常是字母和数字的组合，例如“B123A456”。
计算基础图幅号：
根据经纬度坐标来确定图幅号的基础编号。例如，经度和纬度分别决定图幅号的列和行。行标识符通常由字母组成（A-Z），列号由数字组成。
分幅编号的细化：
根据比例尺的不同，图幅号进一步细分。例如：
1:1000000 比例尺下，图幅号如“B123”。
1:500000 比例尺下，图幅号如“B123A456”。
各个比例尺的图幅号计算规则不同，需要根据具体的比例尺和图幅号标准进行计算。
在 ArcGIS 中，使用 Python 脚本可以高效地计算要素的图幅号。以下是实现这一功能的基本步骤：
导入所需模块：
# -*- coding: gbk -*-import arcpyimport math
将十进制度转换为度、分、秒格式：
def decimal_to_dms(decimal):degrees = int(decimal)minutes = int((abs(decimal) - abs(degrees)) * 60)seconds = (abs(decimal) - abs(degrees)) * 3600 - minutes * 60return "{}°{}′{:.2f}″".format(degrees, minutes, seconds)
根据经度、纬度和比例尺计算图幅编号：
def calculate_map_number(longitude, latitude, scale):"""根据经度、纬度和比例尺计算图幅编号。"""row_characters = "ABCDEFGHIJKLMNOPQRSTUV"row = row_characters[int(math.floor(latitude / 4))]col = int(math.floor(longitude / 6)) + 31base_number = "{}{}".format(row, col)if scale == 1000000:return base_numberelif scale == 500000:B_number = "B{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) / 2.0) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 3.0) + 1)))[-3:])return "{}{}".format(base_number, B_number)elif scale == 250000:C_number = "C{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) / 1.0) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 1.5) + 1)))[-3:])return "{}{}".format(base_number, C_number)elif scale == 100000:D_number = "D{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 3) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 0.5) + 1)))[-3:])return "{}{}".format(base_number, D_number)elif scale == 50000:E_number = "E{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 6) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 0.25) + 1)))[-3:])return "{}{}".format(base_number, E_number)elif scale == 25000:F_number = "F{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 12) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 0.125) + 1)))[-3:])return "{}{}".format(base_number, F_number)elif scale == 10000:G_number = "G{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 24) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 0.0625) + 1)))[-3:])return "{}{}".format(base_number, G_number)elif scale == 5000:H_number = "H{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 48) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) / 0.03125) + 1)))[-3:])return "{}{}".format(base_number, H_number)elif scale == 2000:I_number = "I{}{}".format(str(1000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 144) + 1)))[-3:],str(1000 + (int(((longitude - int(longitude / 6) * 6) * 96) + 1)))[-3:])return "{}{}".format(base_number, I_number)elif scale == 1000:J_number = "J{}{}".format(str(10000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 288) + 1)))[-4:],str(10000 + (int(((longitude - int(longitude / 6) * 6) * 192) + 1)))[-4:])return "{}{}".format(base_number, J_number)elif scale == 500:K_number = "K{}{}".format(str(10000 + (int(((math.ceil(latitude / 4.0) * 4 - latitude) * 576) + 1)))[-4:],str(10000 + (int(((longitude - int(longitude / 6) * 6) * 384) + 1)))[-4:])return "{}{}".format(base_number, K_number)else:raise ValueError("不支持的比例尺: {}".format(scale))
如字段不存在则添加记录图幅号字段：
def add_map_number_field(fc, field_name):if field_name not in [f.name for f in arcpy.ListFields(fc)]:arcpy.AddField_management(fc, field_name, "TEXT", field_length=255)
投影到2000地理坐标并记录临时路径：
def project_feature_class(input_fc, spatial_reference):temp_fc = arcpy.CreateUniqueName("temp_fc_projected", arcpy.env.scratchGDB) # 使用 scratchGDB 工作空间arcpy.Project_management(input_fc, temp_fc, spatial_reference)return temp_fc
临时要素的值返回到输入要素：
def update_field_from_temp_fc(input_fc, temp_fc, field_name):temp_oid_field = "OBJECTID"input_oid_field = "OBJECTID"# 使用游标遍历临时要素类的每个要素，并将字段值更新到输入要素类中with arcpy.da.UpdateCursor(temp_fc, [temp_oid_field, field_name]) as temp_cursor:temp_data = {row[0]: row[1] for row in temp_cursor}with arcpy.da.UpdateCursor(input_fc, [input_oid_field, field_name]) as input_cursor:for row in input_cursor:oid = row[0]if oid in temp_data:row[1] = temp_data[oid]input_cursor.updateRow(row)
避免线要素有遗漏而生成的点：
def create_sample_points_line(geometry, density=100):points = []length = geometry.length # 获取线的长度step_length = length / density # 根据密度计算步长distance = 0arcpy.AddMessage("线的长度: {}".format(length))while distance <= length:point_geom = geometry.positionAlongLine(distance) # 获取沿线指定距离的点point = point_geom.firstPoint # 从PointGeometry中提取Pointpoints.append(point)distance += step_lengthreturn points
避免面要素有遗漏而生成的点：
def create_sample_points(geometry, density=100):points = []extent = geometry.extent # 获取多边形的范围x_min, x_max, y_min, y_max = extent.XMin, extent.XMax, extent.YMin, extent.YMaxarcpy.AddMessage("范围: x_min={}, x_max={}, y_min={}, y_max={}".format(x_min, x_max, y_min, y_max))if x_max > x_min and y_max > y_min:x_range = x_max - x_miny_range = y_max - y_minstep_x = x_range / density # 根据密度计算x方向的步长step_y = y_range / density # 根据密度计算y方向的步长step = max(step_x, step_y)arcpy.AddMessage("计算的步长: step_x={}, step_y={}, step={}".format(step_x, step_y, step))if step <= 0:step = 1else:step = 1# 在多边形内部生成点x_current = x_minwhile x_current <= x_max:y_current = y_minwhile y_current <= y_max:point = arcpy.Point(x_current, y_current)if geometry.contains(point): # 检查点是否在多边形内points.append(point)y_current += stepx_current += step# 获取多边形的所有边界点for part in geometry:for point in part:if point and point not in points: # 检查点是否为None并避免重复添加points.append(point)return points
主函数调用逻辑：
def main():input_fc = arcpy.GetParameterAsText(0) # 输入要素类scale = int(arcpy.GetParameterAsText(1)) # 比例尺 (例如，1000000表示1:1,000,000)field_name = "MapNumber"# 定义CGCS2000的空间参考 (EPSG:4490)spatial_reference = arcpy.SpatialReference(4490)# 为图幅编号添加字段add_map_number_field(input_fc, field_name)# 如果要素类是投影坐标系，则进行坐标转换if arcpy.Describe(input_fc).spatialReference.type == 'Projected':temp_fc = project_feature_class(input_fc, spatial_reference)is_projected = Trueelse:temp_fc = input_fcis_projected = False# 使用游标遍历要素类的每个要素with arcpy.da.UpdateCursor(temp_fc, ["SHAPE@", field_name]) as cursor:for row in cursor:geometry = row[0]geom_type = geometry.type.lower() # 获取要素的几何类型arcpy.AddMessage("处理几何类型: {}".format(geom_type))# 根据几何类型选择生成采样点的方法if geom_type == "polygon":points = create_sample_points(geometry)elif geom_type == "polyline":points = create_sample_points_line(geometry)elif geom_type == "point":point = geometry.firstPoint # 从PointGeometry中提取Pointpoints = [point]else:arcpy.AddMessage("不支持的几何类型: {}".format(geom_type))continuemap_numbers = set() # 存储生成的图幅编号for point in points:# 获取点的经纬度坐标x, y = point.X, point.Yarcpy.AddMessage("处理点: x={}, y={}".format(x, y))try:map_number = calculate_map_number(x, y, scale)arcpy.AddMessage("计算出的图幅编号: {}".format(map_number))map_numbers.add(map_number)except Exception as e:arcpy.AddMessage("计算图幅编号时出错: {}".format(e))continueif map_numbers:row[1] = "、".join(map_numbers) # 将图幅编号写入字段cursor.updateRow(row)else:arcpy.AddMessage("此要素未计算出图幅编号。")# 清理临时要素类if is_projected:update_field_from_temp_fc(input_fc, temp_fc, field_name)arcpy.Delete_management(temp_fc)if __name__ == "__main__":main()
已经全面支持点、线、面要素的计算，比例尺从1:100万—1:500都可以支持，可以在下拉列表中选择。如果想自己手动用源码建脚本工具可以参考这篇里提到的建脚本工具方法→ArcGIS 批量导出 MXD 地图，高效不加班
通过 Python 脚本计算 ArcGIS 要素的图幅号，不仅能提高工作效率，还能保证编号的准确性和一致性。希望本文提供的代码示例能帮助您在实际工作中实现图幅号的自动化计算。要获取做好的脚本工具在公众号后台回复“图幅号计算”获取。
题外话，我查遍了全网都没有找到比较完整的图幅号计算公开python源码，如果你不懂Arcpy你会满头雾水，只求快快发工具，如果你懂或者你是未来的某天网罗全网的时候中搜到这篇，就知道这篇有多牛了哈哈哈。虽然代码还有些缺陷，不过你们可能暂时发现不了。
