|
|
@@ -1,140 +1,57 @@
|
|
|
-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}")
|
|
|
+from qgis._core import QgsWkbTypes
|
|
|
+from qgis.core import QgsVectorLayer, QgsFeature
|
|
|
+
|
|
|
+# 加载Shapefile
|
|
|
+layer = QgsVectorLayer("E:\\projects\\遥感技术部需求\\裁剪栅格\\data\\grid.shp", "grid", "ogr")
|
|
|
+D = 40
|
|
|
+d = 2
|
|
|
+# 检查图层是否加载成功
|
|
|
+if not layer.isValid():
|
|
|
+ print("图层加载失败！")
|
|
|
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("邻居计算完成！")
|
|
|
+ # 遍历每一个Feature
|
|
|
+ for feature in layer.getFeatures():
|
|
|
+ # 获取每个Feature的几何信息
|
|
|
+ geometry = feature.geometry()
|
|
|
+
|
|
|
+ # 判断几何类型，如果是点类型，可以直接获取坐标
|
|
|
+ if geometry.type() == QgsWkbTypes.PointGeometry:
|
|
|
+ point = geometry.asPoint() # 返回点的坐标 (x, y)
|
|
|
+ print(f"点坐标: {point}")
|
|
|
+
|
|
|
+ # 如果是线或面类型，获取所有坐标
|
|
|
+ elif geometry.type() == QgsWkbTypes.LineGeometry:
|
|
|
+ for point in geometry.asPolyline(): # 获取线的所有点
|
|
|
+ print(f"线坐标: {point}")
|
|
|
+
|
|
|
+ elif geometry.type() == QgsWkbTypes.PolygonGeometry:
|
|
|
+ rings = None
|
|
|
+ try:
|
|
|
+ rings = geometry.asPolygon()
|
|
|
+ except Exception as e:
|
|
|
+ rings = geometry.asMultiPolygon()
|
|
|
+ rings = rings[0]
|
|
|
+ finally:
|
|
|
+ print("读取完成")
|
|
|
+ for ring in rings: # 获取面的所有点（包括外环和内环）
|
|
|
+ if len(ring) != 5:
|
|
|
+ print("要素错误！")
|
|
|
+ continue
|
|
|
+ print("开始计算")
|
|
|
+ p1 = ring[0]
|
|
|
+ p2 = ring[1]
|
|
|
+ p3 = ring[2]
|
|
|
+ p4 = ring[3]
|
|
|
+ # ymax = int((max(p1.x(), p2.x(), p3.x(), p4.x()) + D) / d) * d
|
|
|
+ # xmin = int((min(p1.x(), p2.x(), p3.x(), p4.x()) - D) / d) * d
|
|
|
+ # ymin = int((min(p1.y(), p2.y(), p3.y(), p4.y()) - D) / d) * d
|
|
|
+ # xmax = int((max(p1.y(), p2.y(), p3.y(), p4.y()) + D) / d) * d
|
|
|
+
|
|
|
+ ymax = int((max(p1.y(), p2.y(), p3.y(), p4.y()) + D) / d) * d
|
|
|
+ xmin = int((min(p1.x(), p2.x(), p3.x(), p4.x()) - D) / d) * d
|
|
|
+ ymin = int((min(p1.y(), p2.y(), p3.y(), p4.y()) - D) / d) * d
|
|
|
+ xmax = int((max(p1.x(), p2.x(), p3.x(), p4.x()) + D) / d) * d
|
|
|
+
|
|
|
+ print(f"{xmin},{ymin},{xmax},{ymax}")
|
|
|
+ # for point in ring:
|
|
|
+ # print(f"面坐标: {point}")
|
