test.py 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. from PyQt5.QtCore import QVariant
  2. from qgis.core import *
  3. from qgis.utils import iface
  4. from qgis.core import QgsExpression, QgsFeatureRequest, QgsVectorLayer
  5. # 获取渔网图层(假设图层已经加载)
  6. shp_path = "E:\\projects\\遥感技术部需求\\批量生成元数据\\make\\template\\allfishnet.shp"
  7. layer = QgsVectorLayer(shp_path, "Fish Layer", "ogr")
  8. # 获取网格宽度和高度(假设是规则网格)
  9. extent = layer.extent()
  10. grid_width = 100 # 替换为实际网格宽度
  11. grid_height = 100 # 替换为实际网格高度
  12. # 确保图层有至少一个要素
  13. if layer.featureCount() > 0:
  14. # 获取第一个要素
  15. feature = next(layer.getFeatures())
  16. geom = feature.geometry()
  17. # 获取要素的边界矩形(bounding box)
  18. bounds = geom.boundingBox()
  19. # 计算宽度和高度
  20. grid_width = bounds.width()
  21. grid_height = bounds.height()
  22. print(f"网格宽度:{grid_width}")
  23. print(f"网格高度:{grid_height}")
  24. else:
  25. print("图层没有任何要素!")
  26. # 为图层添加邻居字段
  27. layer.startEditing()
  28. for field in ["east", "west", "south", "north"]:
  29. if field not in [f.name() for f in layer.fields()]:
  30. layer.dataProvider().addAttributes([QgsField(field, QVariant.String)])
  31. layer.updateFields()
  32. # 创建空间索引
  33. spatial_index = QgsSpatialIndex(layer)
  34. # 设置误差值 过滤掉仅仅是边线重叠的数据
  35. miss = 0.04
  36. def cacl(neighbors):
  37. if len(neighbors) == 0:
  38. return "自由"
  39. else:
  40. return None
  41. # 遍历图层中的所有要素
  42. for feature in layer.getFeatures():
  43. geom = feature.geometry()
  44. extent = geom.boundingBox()
  45. # 获取网格宽度和高度
  46. width = grid_width
  47. height = grid_height
  48. # 构造东西南北方向的矩形
  49. west_extent = QgsRectangle(extent.xMinimum() - width + miss, extent.yMinimum() + miss, extent.xMinimum() - miss,
  50. extent.yMaximum() - miss)
  51. east_extent = QgsRectangle(extent.xMaximum() + miss, extent.yMinimum() + miss, extent.xMaximum() + width - miss,
  52. extent.yMaximum() - miss)
  53. south_extent = QgsRectangle(extent.xMinimum() + miss, extent.yMinimum() - height + miss, extent.xMaximum() - miss,
  54. extent.yMinimum() - miss)
  55. north_extent = QgsRectangle(extent.xMinimum() + miss, extent.yMaximum() + miss, extent.xMaximum() - miss,
  56. extent.yMaximum() + height - miss)
  57. # 查询每个方向是否有其他网格
  58. west_neighbors = spatial_index.intersects(west_extent)
  59. east_neighbors = spatial_index.intersects(east_extent)
  60. south_neighbors = spatial_index.intersects(south_extent)
  61. north_neighbors = spatial_index.intersects(north_extent)
  62. neighbor_geometries = {
  63. "north": cacl(north_neighbors),
  64. "south": cacl(south_neighbors),
  65. "west": cacl(west_neighbors),
  66. "east": cacl(east_neighbors)
  67. }
  68. for direction, value in neighbor_geometries.items():
  69. feature[direction] = value
  70. layer.updateFeature(feature)
  71. # 保存编辑
  72. layer.commitChanges()
  73. # 批量处理每个网格
  74. # for feature in layer.getFeatures():
  75. # geom = feature.geometry()
  76. # rectangle = geom.boundingBox()
  77. # centroid = geom.centroid().asPoint()
  78. # # 定义邻居的目标几何
  79. # neighbor_geometries = {
  80. # "north": QgsRectangle(
  81. # rectangle.xMinimum() + 0,
  82. # rectangle.yMinimum() + grid_height,
  83. # rectangle.xMaximum() + 0,
  84. # rectangle.yMaximum() + grid_height
  85. # ),
  86. # "south": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() - 1.5 * grid_height,
  87. # centroid.x() + grid_width / 2, centroid.y() - grid_height / 2),
  88. # "west": QgsRectangle(centroid.x() - 1.5 * grid_width, centroid.y() - grid_height / 2,
  89. # centroid.x() - grid_width / 2, centroid.y() + grid_height / 2),
  90. # "east": QgsRectangle(centroid.x() + grid_width / 2, centroid.y() - grid_height / 2,
  91. # centroid.x() + 1.5 * grid_width, centroid.y() + grid_height / 2)
  92. # }
  93. # # neighbor_geometries = {
  94. # # "north": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() + grid_height / 2,
  95. # # centroid.x() + grid_width / 2, centroid.y() + 1.5 * grid_height),
  96. # # "south": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() - 1.5 * grid_height,
  97. # # centroid.x() + grid_width / 2, centroid.y() - grid_height / 2),
  98. # # "west": QgsRectangle(centroid.x() - 1.5 * grid_width, centroid.y() - grid_height / 2,
  99. # # centroid.x() - grid_width / 2, centroid.y() + grid_height / 2),
  100. # # "east": QgsRectangle(centroid.x() + grid_width / 2, centroid.y() - grid_height / 2,
  101. # # centroid.x() + 1.5 * grid_width, centroid.y() + grid_height / 2)
  102. # # }
  103. # # neighbor_geometries = {
  104. # # "north": QgsRectangle(centroid.x(), centroid.y() + grid_height, centroid.x() + grid_width,
  105. # # centroid.y() + 2 * grid_height),
  106. # # "south": QgsRectangle(centroid.x(), centroid.y() - 2 * grid_height, centroid.x() + grid_width,
  107. # # centroid.y() - grid_height),
  108. # # "west": QgsRectangle(centroid.x() - grid_width, centroid.y() - grid_height, centroid.x(),
  109. # # centroid.y() + grid_height),
  110. # # "east": QgsRectangle(centroid.x() + grid_width, centroid.y() - grid_height, centroid.x() + 2 * grid_width,
  111. # # centroid.y() + grid_height),
  112. # # }
  113. # # 查找邻居
  114. # neighbors = {}
  115. # for direction, neighbor_geom in neighbor_geometries.items():
  116. # candidates = spatial_index.intersects(neighbor_geom)
  117. # neighbors[direction] = int(
  118. # any(layer.getFeature(fid).geometry().intersects(neighbor_geom) for fid in candidates))
  119. # # 更新属性
  120. # for direction, value in neighbors.items():
  121. # if value == 0:
  122. # feature[direction] = "自由"
  123. # else:
  124. # feature[direction] = None
  125. # layer.updateFeature(feature)
  126. # # 保存编辑
  127. # layer.commitChanges()
  128. # print("邻居计算完成!")