123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- from PyQt5.QtCore import QVariant
- from qgis.core import *
- from qgis.utils import iface
- from qgis.core import QgsExpression, QgsFeatureRequest, QgsVectorLayer
- # 获取渔网图层(假设图层已经加载)
- shp_path = "E:\\projects\\遥感技术部需求\\批量生成元数据\\make\\template\\allfishnet.shp"
- layer = QgsVectorLayer(shp_path, "Fish Layer", "ogr")
- # 获取网格宽度和高度(假设是规则网格)
- extent = layer.extent()
- grid_width = 100 # 替换为实际网格宽度
- grid_height = 100 # 替换为实际网格高度
- # 确保图层有至少一个要素
- if layer.featureCount() > 0:
- # 获取第一个要素
- feature = next(layer.getFeatures())
- geom = feature.geometry()
- # 获取要素的边界矩形(bounding box)
- bounds = geom.boundingBox()
- # 计算宽度和高度
- grid_width = bounds.width()
- grid_height = bounds.height()
- print(f"网格宽度:{grid_width}")
- print(f"网格高度:{grid_height}")
- else:
- print("图层没有任何要素!")
- # 为图层添加邻居字段
- layer.startEditing()
- for field in ["east", "west", "south", "north"]:
- if field not in [f.name() for f in layer.fields()]:
- layer.dataProvider().addAttributes([QgsField(field, QVariant.String)])
- layer.updateFields()
- # 创建空间索引
- spatial_index = QgsSpatialIndex(layer)
- # 设置误差值 过滤掉仅仅是边线重叠的数据
- miss = 0.04
- def cacl(neighbors):
- if len(neighbors) == 0:
- return "自由"
- else:
- return None
- # 遍历图层中的所有要素
- for feature in layer.getFeatures():
- geom = feature.geometry()
- extent = geom.boundingBox()
- # 获取网格宽度和高度
- width = grid_width
- height = grid_height
- # 构造东西南北方向的矩形
- west_extent = QgsRectangle(extent.xMinimum() - width + miss, extent.yMinimum() + miss, extent.xMinimum() - miss,
- extent.yMaximum() - miss)
- east_extent = QgsRectangle(extent.xMaximum() + miss, extent.yMinimum() + miss, extent.xMaximum() + width - miss,
- extent.yMaximum() - miss)
- south_extent = QgsRectangle(extent.xMinimum() + miss, extent.yMinimum() - height + miss, extent.xMaximum() - miss,
- extent.yMinimum() - miss)
- north_extent = QgsRectangle(extent.xMinimum() + miss, extent.yMaximum() + miss, extent.xMaximum() - miss,
- extent.yMaximum() + height - miss)
- # 查询每个方向是否有其他网格
- west_neighbors = spatial_index.intersects(west_extent)
- east_neighbors = spatial_index.intersects(east_extent)
- south_neighbors = spatial_index.intersects(south_extent)
- north_neighbors = spatial_index.intersects(north_extent)
- neighbor_geometries = {
- "north": cacl(north_neighbors),
- "south": cacl(south_neighbors),
- "west": cacl(west_neighbors),
- "east": cacl(east_neighbors)
- }
- for direction, value in neighbor_geometries.items():
- feature[direction] = value
- layer.updateFeature(feature)
- # 保存编辑
- layer.commitChanges()
- # 批量处理每个网格
- # for feature in layer.getFeatures():
- # geom = feature.geometry()
- # rectangle = geom.boundingBox()
- # centroid = geom.centroid().asPoint()
- # # 定义邻居的目标几何
- # neighbor_geometries = {
- # "north": QgsRectangle(
- # rectangle.xMinimum() + 0,
- # rectangle.yMinimum() + grid_height,
- # rectangle.xMaximum() + 0,
- # rectangle.yMaximum() + grid_height
- # ),
- # "south": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() - 1.5 * grid_height,
- # centroid.x() + grid_width / 2, centroid.y() - grid_height / 2),
- # "west": QgsRectangle(centroid.x() - 1.5 * grid_width, centroid.y() - grid_height / 2,
- # centroid.x() - grid_width / 2, centroid.y() + grid_height / 2),
- # "east": QgsRectangle(centroid.x() + grid_width / 2, centroid.y() - grid_height / 2,
- # centroid.x() + 1.5 * grid_width, centroid.y() + grid_height / 2)
- # }
- # # neighbor_geometries = {
- # # "north": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() + grid_height / 2,
- # # centroid.x() + grid_width / 2, centroid.y() + 1.5 * grid_height),
- # # "south": QgsRectangle(centroid.x() - grid_width / 2, centroid.y() - 1.5 * grid_height,
- # # centroid.x() + grid_width / 2, centroid.y() - grid_height / 2),
- # # "west": QgsRectangle(centroid.x() - 1.5 * grid_width, centroid.y() - grid_height / 2,
- # # centroid.x() - grid_width / 2, centroid.y() + grid_height / 2),
- # # "east": QgsRectangle(centroid.x() + grid_width / 2, centroid.y() - grid_height / 2,
- # # centroid.x() + 1.5 * grid_width, centroid.y() + grid_height / 2)
- # # }
- # # neighbor_geometries = {
- # # "north": QgsRectangle(centroid.x(), centroid.y() + grid_height, centroid.x() + grid_width,
- # # centroid.y() + 2 * grid_height),
- # # "south": QgsRectangle(centroid.x(), centroid.y() - 2 * grid_height, centroid.x() + grid_width,
- # # centroid.y() - grid_height),
- # # "west": QgsRectangle(centroid.x() - grid_width, centroid.y() - grid_height, centroid.x(),
- # # centroid.y() + grid_height),
- # # "east": QgsRectangle(centroid.x() + grid_width, centroid.y() - grid_height, centroid.x() + 2 * grid_width,
- # # centroid.y() + grid_height),
- # # }
- # # 查找邻居
- # neighbors = {}
- # for direction, neighbor_geom in neighbor_geometries.items():
- # candidates = spatial_index.intersects(neighbor_geom)
- # neighbors[direction] = int(
- # any(layer.getFeature(fid).geometry().intersects(neighbor_geom) for fid in candidates))
- # # 更新属性
- # for direction, value in neighbors.items():
- # if value == 0:
- # feature[direction] = "自由"
- # else:
- # feature[direction] = None
- # layer.updateFeature(feature)
- # # 保存编辑
- # layer.commitChanges()
- # print("邻居计算完成!")
|