使用E:\Pro\kk\last\forest.gdb\ds\用地范围,生成图幅，2970*4200m，图幅，重合率10%，结果放在桌面GIS文件夹,tf.shp
import arcpyimport os# 输入参数input_fc = r'E:\Pro\kk\last\forest.gdb\ds\用地范围'output_fc = r'C:\Users\yl\Desktop\GIS\tf.shp'grid_width = 2970 # 米grid_height = 4200 # 米overlap = 0.10 # 10%重合率# 确保输出目录存在output_dir = os.path.dirname(output_fc)if not os.path.exists(output_dir): os.makedirs(output_dir)# 获取输入要素的范围desc = arcpy.Describe(input_fc)extent = desc.extentprint(f"输入范围：{extent}")min_x = extent.XMinmax_x = extent.XMaxmin_y = extent.YMinmax_y = extent.YMaxprint(f"X: [{min_x:.2f}, {max_x:.2f}]")print(f"Y: [{min_y:.2f}, {max_y:.2f}]")# 计算实际步长（减去重叠部分）step_x = grid_width * (1 - overlap)step_y = grid_height * (1 - overlap)print(f"网格大小：{grid_width} x {grid_height} 米")print(f"步长：{step_x:.2f} x {step_y:.2f} 米 (重合率 {overlap*100}%)")# 获取空间参考spatial_ref = desc.spatialReference# 创建输出要素类if arcpy.Exists(output_fc): arcpy.Delete_management(output_fc)arcpy.CreateFeatureclass_management( os.path.dirname(output_fc), os.path.basename(output_fc),'POLYGON', spatial_reference=spatial_ref)# 添加字段arcpy.AddField_management(output_fc, 'ID', 'LONG')arcpy.AddField_management(output_fc, 'X_MIN', 'DOUBLE')arcpy.AddField_management(output_fc, 'X_MAX', 'DOUBLE')arcpy.AddField_management(output_fc, 'Y_MIN', 'DOUBLE')arcpy.AddField_management(output_fc, 'Y_MAX', 'DOUBLE')arcpy.AddField_management(output_fc, 'ROW', 'LONG')arcpy.AddField_management(output_fc, 'COL', 'LONG')# 读取用地范围的所有几何用于相交检查print("正在读取用地范围几何...")boundary_geoms = []with arcpy.da.SearchCursor(input_fc, ['SHAPE@']) as cursor:for row in cursor: boundary_geoms.append(row[0])print(f"用地范围几何数量：{len(boundary_geoms)}")# 生成网格，只保留与用地范围相交的图幅grid_id = 0row_idx = 0y_current = min_ycursor = arcpy.da.InsertCursor(output_fc, ['SHAPE@', 'ID', 'X_MIN', 'X_MAX', 'Y_MIN', 'Y_MAX', 'ROW', 'COL'])total_rows = int((max_y - min_y) / step_y) + 2total_cols = int((max_x - min_x) / step_x) + 2print(f"预计生成约 {total_rows * total_cols} 个图幅...")while y_current < max_y: col_idx = 0x_current = min_xwhile x_current < max_x:# 创建网格多边形array = arcpy.Array([ arcpy.Point(x_current, y_current), arcpy.Point(x_current + grid_width, y_current), arcpy.Point(x_current + grid_width, y_current + grid_height), arcpy.Point(x_current, y_current + grid_height), arcpy.Point(x_current, y_current) ]) polygon = arcpy.Polygon(array, spatial_ref)# 检查是否与任一用地范围几何相交intersects = False for boundary_geom in boundary_geoms:if not polygon.disjoint(boundary_geom): intersects = True break if intersects:# 相交，保存该图幅cursor.insertRow([ polygon, grid_id, x_current, x_current + grid_width, y_current, y_current + grid_height, row_idx, col_idx ]) grid_id += 1col_idx += 1x_current += step_x row_idx += 1y_current += step_yif row_idx % 10 == 0: print(f"已处理 {row_idx} 行...")del cursorprint(f"\n完成！共生成 {grid_id} 个与用地范围相交的完整图幅（保持 2970 x 4200 米）")print(f"输出文件：{output_fc}")
夜雨聆风
