夜雨聆风# -*- coding: utf-8 -*-import arcpyimport osimport math# 输入参数input_fc = r'E:\2026\Pro\data\data\china\2024年数据.gdb\ds\T2024年省级'output_fc = os.path.join(os.environ['USERPROFILE'], 'Desktop', 'GIS', '经纬网.shp')interval = 5 # 经纬度间隔：5度buffer_distance = 80000 # 80公里外扩# 获取输入要素的信息desc = arcpy.Describe(input_fc)spatial_ref = desc.spatialReferenceprint(f"坐标系: {spatial_ref.name}")# 获取要素范围extent = desc.extentprint(f"数据范围: XMin={extent.XMin:.2f}, YMin={extent.YMin:.2f}, XMax={extent.XMax:.2f}, YMax={extent.YMax:.2f}")# 外扩80公里x_min = extent.XMin - buffer_distancey_min = extent.YMin - buffer_distancex_max = extent.XMax + buffer_distancey_max = extent.YMax + buffer_distance# 定义地理坐标系（WGS84或CGCS2000）# 使用CGCS2000地理坐标系（与中国数据匹配）geographic_sr = arcpy.SpatialReference(4490) # CGCS2000地理坐标系# 投影坐标系的范围转换为地理坐标# 创建范围多边形并转换到地理坐标系extent_polygon = arcpy.Polyline(arcpy.Array([ arcpy.Point(x_min, y_min), arcpy.Point(x_max, y_min), arcpy.Point(x_max, y_max), arcpy.Point(x_min, y_max), arcpy.Point(x_min, y_min)]), spatial_ref)# 投影到地理坐标系获取经纬度范围geo_extent_polygon = arcpy.Polygon(arcpy.Array([ arcpy.Point(x_min, y_min), arcpy.Point(x_max, y_min), arcpy.Point(x_max, y_max), arcpy.Point(x_min, y_max)]), spatial_ref)# 获取四个角的经纬度corners = [ (x_min, y_min), (x_min, y_max), (x_max, y_min), (x_max, y_max)]# 转换坐标min_lon, max_lon = float('inf'), float('-inf')min_lat, max_lat = float('inf'), float('-inf')point = arcpy.Point()for x, y in corners: point.X = x point.Y = y point_geom = arcpy.PointGeometry(point, spatial_ref) point_geo = point_geom.projectAs(geographic_sr) lon, lat = point_geo.firstPoint.X, point_geo.firstPoint.Y min_lon = min(min_lon, lon) max_lon = max(max_lon, lon) min_lat = min(min_lat, lat) max_lat = max(max_lat, lat) print(f"角点投影坐标 ({x:.0f}, {y:.0f}) -> 经纬度 ({lon:.2f}, {lat:.2f})")print(f"\n经纬度范围: Lon [{min_lon:.2f}, {max_lon:.2f}], Lat [{min_lat:.2f}, {max_lat:.2f}]")# 计算经纬网起点（取整到5度的倍数）start_lon = math.floor(min_lon / interval) * intervalstart_lat = math.floor(min_lat / interval) * intervalend_lon = math.ceil(max_lon / interval) * intervalend_lat = math.ceil(max_lat / interval) * intervalprint(f"经纬网范围: Lon [{start_lon}, {end_lon}], Lat [{start_lat}, {end_lat}]")print(f"经纬度间隔: {interval}度")# 创建输出要素类arcpy.env.overwriteOutput = Trueoutput_dir = os.path.dirname(output_fc)if not os.path.exists(output_dir): os.makedirs(output_dir)# 创建空的线要素类（使用地理坐标系）arcpy.CreateFeatureclass_management(output_dir, os.path.basename(output_fc), "POLYLINE", spatial_reference=geographic_sr)# 添加字段arcpy.AddField_management(output_fc, "DIRECTION", "TEXT", field_length=10) # 方向：NS（经线）或EW（纬线）arcpy.AddField_management(output_fc, "DEGREE", "DOUBLE") # 经度或纬度值arcpy.AddField_management(output_fc, "LABEL", "TEXT", field_length=20) # 标注文本# 生成经纬网线lines = []# 生成经线（南北向，固定经度）lon_current = start_loncount_meridian = 0while lon_current <= end_lon:# 经线从南到北line = arcpy.Polyline(arcpy.Array([ arcpy.Point(lon_current, start_lat), arcpy.Point(lon_current, end_lat) ]), geographic_sr)# 生成标注if lon_current < 0: label = f"{abs(lon_current):.0f}°W"elif lon_current > 0: label = f"{lon_current:.0f}°E"else: label = "0°"lines.append((line, "NS", lon_current, label)) lon_current += interval count_meridian += 1# 生成纬线（东西向，固定纬度）lat_current = start_latcount_parallel = 0while lat_current <= end_lat:# 纬线从西到东line = arcpy.Polyline(arcpy.Array([ arcpy.Point(start_lon, lat_current), arcpy.Point(end_lon, lat_current) ]), geographic_sr)# 生成标注if lat_current < 0: label = f"{abs(lat_current):.0f}°S"elif lat_current > 0: label = f"{lat_current:.0f}°N"else: label = "0°"lines.append((line, "EW", lat_current, label)) lat_current += interval count_parallel += 1print(f"\n生成经线: {count_meridian}条")print(f"生成纬线: {count_parallel}条")print(f"总计: {len(lines)}条")# 写入要素类with arcpy.da.InsertCursor(output_fc, ["SHAPE@", "DIRECTION", "DEGREE", "LABEL"]) as cursor:for line, direction, degree, label in lines: cursor.insertRow([line, direction, degree, label])print(f"\n经纬网已生成: {output_fc}")print("完成！")
