夜雨聆风import arcpyimport os# 设置工作空间desktop = os.path.join(os.path.expanduser("~"), "Desktop", "GIS")os.makedirs(desktop, exist_ok=True)arcpy.env.workspace = desktoparcpy.env.overwriteOutput = True# 输出路径output_path = os.path.join(arcpy.env.workspace, "world_diagonal_lines_5deg.shp")# 第一条线：从 (-180, -90) 到 (180, 90)# 斜率 = (90 - (-90)) / (180 - (-180)) = 180 / 360 = 0.5# 方程：lat = 0.5 * lon + b# 第一条线通过 (-180, -90)：-90 = 0.5 * (-180) + b => b = 0# 所以第一条线：lat = 0.5 * lon# 平行线间距：5 度（纬度方向偏移）lat_offset = 5# 删除已存在的文件if arcpy.Exists(output_path): arcpy.Delete_management(output_path) print("已删除旧文件")# 创建线要素类arcpy.CreateFeatureclass_management( out_path=arcpy.env.workspace, out_name="world_diagonal_lines_5deg.shp", geometry_type="POLYLINE", spatial_reference=arcpy.SpatialReference(4326) # WGS84)# 添加字段arcpy.AddField_management(output_path, "jd", "DOUBLE")arcpy.AddField_management(output_path, "wd", "DOUBLE")arcpy.AddField_management(output_path, "offset", "DOUBLE")arcpy.AddField_management(output_path, "direction", "TEXT", field_length=10)print("要素类已创建，开始添加对角线...")# 计算线与地球边界的交点def calc_line_points(slope, b):"""计算斜线与地球边界的交点"""points = []# 左边界 (lon = -180): lat = slope * (-180) + blat_left = slope * (-180) + bif -90 <= lat_left <= 90: points.append((-180, lat_left))# 右边界 (lon = 180): lat = slope * 180 + blat_right = slope * 180 + bif -90 <= lat_right <= 90: points.append((180, lat_right))# 下边界 (lat = -90): lon = (-90 - b) / slopeif slope != 0: lon_bottom = (-90 - b) / slopeif -180 <= lon_bottom <= 180: points.append((lon_bottom, -90))# 上边界 (lat = 90): lon = (90 - b) / slopeif slope != 0: lon_top = (90 - b) / slopeif -180 <= lon_top <= 180: points.append((lon_top, 90))# 去重并排序points = list(set(points))if len(points) >= 2: points.sort(key=lambda p: p[0])return points# 使用插入游标添加要素with arcpy.da.InsertCursor(output_path, ["SHAPE@", "jd", "wd", "offset", "direction"]) as cursor: line_count = 0# ========== 第一组：从 (-180, -90) 到 (180, 90)，斜率 0.5 ==========print("正在生成第一组对角线（左下→右上）...") slope1 = 0.5for b in range(-180, 181, lat_offset): points = calc_line_points(slope1, b)if len(points) >= 2: array = arcpy.Array([arcpy.Point(lon, lat) for lon, lat in points]) polyline = arcpy.Polyline(array, arcpy.SpatialReference(4326)) mid_lon = (points[0][0] + points[-1][0]) / 2mid_lat = (points[0][1] + points[-1][1]) / 2cursor.insertRow([polyline, mid_lon, mid_lat, b, "NE-SW"]) line_count += 1# ========== 第二组：从 (-180, 90) 到 (180, -90)，斜率 -0.5 ==========print("正在生成第二组对角线（左上→右下）...") slope2 = -0.5# 第二条线通过 (-180, 90)：90 = -0.5 * (-180) + b => b = 0for b in range(-180, 181, lat_offset): points = calc_line_points(slope2, b)if len(points) >= 2: array = arcpy.Array([arcpy.Point(lon, lat) for lon, lat in points]) polyline = arcpy.Polyline(array, arcpy.SpatialReference(4326)) mid_lon = (points[0][0] + points[-1][0]) / 2mid_lat = (points[0][1] + points[-1][1]) / 2cursor.insertRow([polyline, mid_lon, mid_lat, b, "NW-SE"]) line_count += 1# 统计结果result_count = int(arcpy.GetCount_management(output_path)[0])# 分别统计两组线layer1 = arcpy.SelectLayerByAttribute_management(output_path, "NEW_SELECTION", "direction = 'NE-SW'")count1 = int(arcpy.GetCount_management(layer1)[0])layer2 = arcpy.SelectLayerByAttribute_management(output_path, "NEW_SELECTION", "direction = 'NW-SE'")count2 = int(arcpy.GetCount_management(layer2)[0])arcpy.SelectLayerByAttribute_management(output_path, "CLEAR_SELECTION")print(f"\n[OK] 文件已生成：{output_path}")print(f"[OK] 对角线总数：{result_count} 条")print(f"[OK] 第一组（左下→右上）：{count1} 条，斜率 0.5，第一条线：(-180, -90) 到 (180, 90)")print(f"[OK] 第二组（左上→右下）：{count2} 条，斜率 -0.5，第一条线：(-180, 90) 到 (180, -90)")print(f"[OK] 平行线间距：{lat_offset} 度（纬度方向）")print(f"[OK] 字段：jd (经度), wd (纬度), offset (截距), direction (方向)")print(f"[OK] 坐标系统：EPSG:4326 (WGS84)")
