本文整理汇总了Python中qgis.core.QgsSpatialIndex.intersects方法的典型用法代码示例。如果您正苦于以下问题:Python QgsSpatialIndex.intersects方法的具体用法?Python QgsSpatialIndex.intersects怎么用?Python QgsSpatialIndex.intersects使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类qgis.core.QgsSpatialIndex
的用法示例。
在下文中一共展示了QgsSpatialIndex.intersects方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: __init__
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
class LayerIndex:
def __init__(self, layer):
self.__layer = layer
self.__index = QgsSpatialIndex()
feats = vector.features(layer)
for ft in feats:
self.__index.insertFeature(ft)
def contains(self, point):
"""Return true if the point intersects the layer"""
intersects = self.__index.intersects(point.boundingBox())
for i in intersects:
request = QgsFeatureRequest().setFilterFid(i)
feat = self.__layer.getFeatures(request).next()
tmpGeom = QgsGeometry(feat.geometry())
if point.intersects(tmpGeom):
return True
return False
def countIntersection(self,bufferGeom,nb):
"""Return true if the buffer intersects enough entities"""
count = 0
intersects = self.__index.intersects(bufferGeom.boundingBox())
for i in intersects:
request = QgsFeatureRequest().setFilterFid(i)
feat = self.__layer.getFeatures(request).next()
tmpGeom = QgsGeometry(feat.geometry())
if bufferGeom.intersects(tmpGeom):
count += 1
if count >= nb:
return True
return False
示例2: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
pointCount = self.parameterAsDouble(parameters, self.POINTS_NUMBER, context)
minDistance = self.parameterAsDouble(parameters, self.MIN_DISTANCE, context)
bbox = source.sourceExtent()
sourceIndex = QgsSpatialIndex(source, feedback)
fields = QgsFields()
fields.append(QgsField('id', QVariant.Int, '', 10, 0))
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, QgsWkbTypes.Point, source.sourceCrs())
nPoints = 0
nIterations = 0
maxIterations = pointCount * 200
total = 100.0 / pointCount if pointCount else 1
index = QgsSpatialIndex()
points = dict()
random.seed()
while nIterations < maxIterations and nPoints < pointCount:
if feedback.isCanceled():
break
rx = bbox.xMinimum() + bbox.width() * random.random()
ry = bbox.yMinimum() + bbox.height() * random.random()
p = QgsPointXY(rx, ry)
geom = QgsGeometry.fromPointXY(p)
ids = sourceIndex.intersects(geom.buffer(5, 5).boundingBox())
if len(ids) > 0 and \
vector.checkMinDistance(p, index, minDistance, points):
request = QgsFeatureRequest().setFilterFids(ids).setSubsetOfAttributes([])
for f in source.getFeatures(request):
if feedback.isCanceled():
break
tmpGeom = f.geometry()
if geom.within(tmpGeom):
f = QgsFeature(nPoints)
f.initAttributes(1)
f.setFields(fields)
f.setAttribute('id', nPoints)
f.setGeometry(geom)
sink.addFeature(f, QgsFeatureSink.FastInsert)
index.insertFeature(f)
points[nPoints] = p
nPoints += 1
feedback.setProgress(int(nPoints * total))
nIterations += 1
if nPoints < pointCount:
feedback.pushInfo(self.tr('Could not generate requested number of random points. '
'Maximum number of attempts exceeded.'))
return {self.OUTPUT: dest_id}
示例3: poly2nb
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def poly2nb(self):
lst = []
index = QgsSpatialIndex()
featsA = self.lyr.getFeatures()
featsB = self.lyr.getFeatures()
for ft in featsA:
index.insertFeature(ft)
featB = QgsFeature()
prv = self.lyr.dataProvider()
while featsB.nextFeature(featB):
geomB = featB.constGeometry()
idb = featB.id()
idxs = index.intersects(geomB.boundingBox())
sor = []
for idx in idxs:
rqst = QgsFeatureRequest().setFilterFid(idx)
featA = prv.getFeatures(rqst).next()
ida = featA.id()
geomA = QgsGeometry(featA.geometry())
if idb!=ida:
if geomB.touches(geomA)==True:
sor.append(ida)
lst.append(sor)
return lst
示例4: __init__
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
class TriangleMesh:
# 0 - 3
# | / |
# 1 - 2
def __init__(self, xmin, ymin, xmax, ymax, x_segments, y_segments):
self.flen = 0
self.quadrangles = []
self.spatial_index = QgsSpatialIndex()
xres = (xmax - xmin) / x_segments
yres = (ymax - ymin) / y_segments
for y in range(y_segments):
for x in range(x_segments):
pt0 = QgsPoint(xmin + x * xres, ymax - y * yres)
pt1 = QgsPoint(xmin + x * xres, ymax - (y + 1) * yres)
pt2 = QgsPoint(xmin + (x + 1) * xres, ymax - (y + 1) * yres)
pt3 = QgsPoint(xmin + (x + 1) * xres, ymax - y * yres)
self._addQuadrangle(pt0, pt1, pt2, pt3)
def _addQuadrangle(self, pt0, pt1, pt2, pt3):
f = QgsFeature(self.flen)
f.setGeometry(QgsGeometry.fromPolygon([[pt0, pt1, pt2, pt3, pt0]]))
self.quadrangles.append(f)
self.spatial_index.insertFeature(f)
self.flen += 1
def intersects(self, geom):
for fid in self.spatial_index.intersects(geom.boundingBox()):
quad = self.quadrangles[fid].geometry()
if quad.intersects(geom):
yield quad
def splitPolygons(self, geom):
for quad in self.intersects(geom):
pts = quad.asPolygon()[0]
tris = [[[pts[0], pts[1], pts[3], pts[0]]], [[pts[3], pts[1], pts[2], pts[3]]]]
if geom.contains(quad):
yield tris[0]
yield tris[1]
else:
for i, tri in enumerate(map(QgsGeometry.fromPolygon, tris)):
if geom.contains(tri):
yield tris[i]
elif geom.intersects(tri):
poly = geom.intersection(tri)
if poly.isMultipart():
for sp in poly.asMultiPolygon():
yield sp
else:
yield poly.asPolygon()
示例5: compute_graph
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def compute_graph(features, feedback, create_id_graph=False, min_distance=0):
""" compute topology from a layer/field """
s = Graph(sort_graph=False)
id_graph = None
if create_id_graph:
id_graph = Graph(sort_graph=True)
# skip features without geometry
features_with_geometry = {f_id: f for (f_id, f) in features.items() if f.hasGeometry()}
total = 70.0 / len(features_with_geometry) if features_with_geometry else 1
index = QgsSpatialIndex()
i = 0
for feature_id, f in features_with_geometry.items():
if feedback.isCanceled():
break
g = f.geometry()
if min_distance > 0:
g = g.buffer(min_distance, 5)
engine = QgsGeometry.createGeometryEngine(g.constGet())
engine.prepareGeometry()
feature_bounds = g.boundingBox()
# grow bounds a little so we get touching features
feature_bounds.grow(feature_bounds.width() * 0.01)
intersections = index.intersects(feature_bounds)
for l2 in intersections:
f2 = features_with_geometry[l2]
if engine.intersects(f2.geometry().constGet()):
s.add_edge(f.id(), f2.id())
s.add_edge(f2.id(), f.id())
if id_graph:
id_graph.add_edge(f.id(), f2.id())
index.insertFeature(f)
i += 1
feedback.setProgress(int(i * total))
for feature_id, f in features_with_geometry.items():
if feedback.isCanceled():
break
if feature_id not in s.node_edge:
s.add_edge(feature_id, None)
return s, id_graph
示例6: LayerIndex
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
class LayerIndex(object):
"""Check an intersection between a QgsGeometry and a QgsVectorLayer."""
def __init__(self, layer):
self.__layer = layer
if QGis.QGIS_VERSION_INT >= 20700:
self.__index = QgsSpatialIndex(layer.getFeatures())
else:
self.__index = QgsSpatialIndex()
for ft in layer.getFeatures():
self.__index.insertFeature(ft)
def contains(self, point):
"""Return true if the point intersects the layer."""
intersects = self.__index.intersects(point.boundingBox())
for i in intersects:
request = QgsFeatureRequest().setFilterFid(i)
feat = self.__layer.getFeatures(request).next()
if point.intersects(QgsGeometry(feat.geometry())):
return True
return False
def count_intersection(self, buffer_geom, nb):
"""Return true if the buffer intersects enough entities."""
count = 0
intersects = self.__index.intersects(buffer_geom.boundingBox())
for i in intersects:
request = QgsFeatureRequest().setFilterFid(i)
feat = self.__layer.getFeatures(request).next()
if buffer_geom.intersects(QgsGeometry(feat.geometry())):
count += 1
if count >= nb:
return True
return False
示例7: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
sourceA = self.parameterAsSource(parameters, self.INPUT, context)
sourceB = self.parameterAsSource(parameters, self.OVERLAY, context)
geomType = QgsWkbTypes.multiType(sourceA.wkbType())
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
sourceA.fields(), geomType, sourceA.sourceCrs())
featB = QgsFeature()
outFeat = QgsFeature()
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs(), context.transformContext())), feedback)
total = 100.0 / (sourceA.featureCount() * sourceB.featureCount()) if sourceA.featureCount() and sourceB.featureCount() else 1
count = 0
for featA in sourceA.getFeatures():
if feedback.isCanceled():
break
if featA.hasGeometry():
geom = featA.geometry()
diffGeom = QgsGeometry(geom)
attrs = featA.attributes()
intersects = indexB.intersects(geom.boundingBox())
request = QgsFeatureRequest().setFilterFids(intersects).setSubsetOfAttributes([])
request.setDestinationCrs(sourceA.sourceCrs(), context.transformContext())
for featB in sourceB.getFeatures(request):
if feedback.isCanceled():
break
tmpGeom = featB.geometry()
if diffGeom.intersects(tmpGeom):
diffGeom = QgsGeometry(diffGeom.difference(tmpGeom))
outFeat.setGeometry(diffGeom)
outFeat.setAttributes(attrs)
sink.addFeature(outFeat, QgsFeatureSink.FastInsert)
else:
sink.addFeature(featA, QgsFeatureSink.FastInsert)
count += 1
feedback.setProgress(int(count * total))
return {self.OUTPUT: dest_id}
示例8: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, progress):
fileName = self.getParameterValue(self.INPUT)
layer = dataobjects.getObjectFromUri(fileName)
fieldName = self.getParameterValue(self.FIELD)
value = self.getParameterValue(self.VALUE)
selected = layer.selectedFeaturesIds()
if len(selected) == 0:
GeoAlgorithmExecutionException(
self.tr('There is no selection in the input layer. '
'Select one feature and try again.'))
ft = layer.selectedFeatures()[0]
geom = ft.geometry()
attrSum = ft[fieldName]
idx = QgsSpatialIndex(layer.getFeatures(QgsFeatureRequest.setSubsetOfAttributes([])))
req = QgsFeatureRequest()
completed = False
while not completed:
intersected = idx.intersects(geom.boundingBox())
if len(intersected) < 0:
progress.setInfo(self.tr('No adjacent features found.'))
break
req = QgsFeatureRequest().setFilterFids(intersected).setSubsetOfAttributes([fieldName], layer.fields())
for ft in layer.getFeatures(req):
tmpGeom = ft.geometry()
if tmpGeom.touches(geom):
geom = tmpGeom.combine(geom)
selected.append(i)
attrSum += ft[fieldName]
if attrSum >= value:
completed = True
break
layer.selectByIds(selected)
self.setOutputValue(self.OUTPUT, fileName)
示例9: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
line_source = self.parameterAsSource(parameters, self.LINES, context)
sameLayer = parameters[self.INPUT] == parameters[self.LINES]
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
source.fields(), QgsWkbTypes.multiType(source.wkbType()), source.sourceCrs())
spatialIndex = QgsSpatialIndex()
splitGeoms = {}
request = QgsFeatureRequest()
request.setSubsetOfAttributes([])
request.setDestinationCrs(source.sourceCrs())
for aSplitFeature in line_source.getFeatures(request):
if feedback.isCanceled():
break
splitGeoms[aSplitFeature.id()] = aSplitFeature.geometry()
spatialIndex.insertFeature(aSplitFeature)
# honor the case that user has selection on split layer and has setting "use selection"
outFeat = QgsFeature()
features = source.getFeatures()
total = 100.0 / source.featureCount() if source.featureCount() else 100
for current, inFeatA in enumerate(features):
if feedback.isCanceled():
break
inGeom = inFeatA.geometry()
attrsA = inFeatA.attributes()
outFeat.setAttributes(attrsA)
if inGeom.isMultipart():
inGeoms = []
for g in inGeom.asGeometryCollection():
inGeoms.append(g)
else:
inGeoms = [inGeom]
lines = spatialIndex.intersects(inGeom.boundingBox())
if len(lines) > 0: # has intersection of bounding boxes
splittingLines = []
engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
engine.prepareGeometry()
for i in lines:
try:
splitGeom = splitGeoms[i]
except:
continue
# check if trying to self-intersect
if sameLayer:
if inFeatA.id() == i:
continue
if engine.intersects(splitGeom.geometry()):
splittingLines.append(splitGeom)
if len(splittingLines) > 0:
for splitGeom in splittingLines:
splitterPList = None
outGeoms = []
split_geom_engine = QgsGeometry.createGeometryEngine(splitGeom.geometry())
split_geom_engine.prepareGeometry()
while len(inGeoms) > 0:
if feedback.isCanceled():
break
inGeom = inGeoms.pop()
if inGeom.isNull(): # this has been encountered and created a run-time error
continue
if split_geom_engine.intersects(inGeom.geometry()):
inPoints = vector.extractPoints(inGeom)
if splitterPList is None:
splitterPList = vector.extractPoints(splitGeom)
try:
result, newGeometries, topoTestPoints = inGeom.splitGeometry(splitterPList, False)
except:
feedback.reportError(self.tr('Geometry exception while splitting'))
result = 1
# splitGeometry: If there are several intersections
# between geometry and splitLine, only the first one is considered.
if result == 0: # split occurred
if inPoints == vector.extractPoints(inGeom):
# bug in splitGeometry: sometimes it returns 0 but
# the geometry is unchanged
#.........这里部分代码省略.........
示例10: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
if source is None:
raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
source.fields(), source.wkbType(), source.sourceCrs())
if sink is None:
raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))
features = source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]))
total = 100.0 / source.featureCount() if source.featureCount() else 0
geoms = dict()
index = QgsSpatialIndex()
for current, f in enumerate(features):
if feedback.isCanceled():
break
geoms[f.id()] = f.geometry()
index.addFeature(f)
feedback.setProgress(int(0.10 * current * total)) # takes about 10% of time
# start by assuming everything is unique, and chop away at this list
unique_features = dict(geoms)
current = 0
for feature_id, geometry in geoms.items():
if feedback.isCanceled():
break
if feature_id not in unique_features:
# feature was already marked as a duplicate
continue
candidates = index.intersects(geometry.boundingBox())
candidates.remove(feature_id)
for candidate_id in candidates:
if candidate_id not in unique_features:
# candidate already marked as a duplicate (not sure if this is possible,
# since it would mean the current feature would also have to be a duplicate!
# but let's be safe!)
continue
if geometry.isGeosEqual(geoms[candidate_id]):
# candidate is a duplicate of feature
del unique_features[candidate_id]
current += 1
feedback.setProgress(int(0.80 * current * total) + 10) # takes about 80% of time
total = 100.0 / len(unique_features) if unique_features else 1
# now, fetch all the feature attributes for the unique features only
# be super-smart and don't re-fetch geometries
request = QgsFeatureRequest().setFilterFids(list(unique_features.keys())).setFlags(QgsFeatureRequest.NoGeometry)
for current, f in enumerate(source.getFeatures(request)):
if feedback.isCanceled():
break
# use already fetched geometry
f.setGeometry(unique_features[f.id()])
sink.addFeature(f, QgsFeatureSink.FastInsert)
feedback.setProgress(int(0.10 * current * total) + 90) # takes about 10% of time
return {self.OUTPUT: dest_id}
示例11: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
poly_source = self.parameterAsSource(parameters, self.POLYGONS, context)
point_source = self.parameterAsSource(parameters, self.POINTS, context)
weight_field = self.parameterAsString(parameters, self.WEIGHT, context)
weight_field_index = -1
if weight_field:
weight_field_index = point_source.fields().lookupField(weight_field)
class_field = self.parameterAsString(parameters, self.CLASSFIELD, context)
class_field_index = -1
if class_field:
class_field_index = point_source.fields().lookupField(class_field)
field_name = self.parameterAsString(parameters, self.FIELD, context)
fields = poly_source.fields()
if fields.lookupField(field_name) < 0:
fields.append(QgsField(field_name, QVariant.Int))
field_index = fields.lookupField(field_name)
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, poly_source.wkbType(), poly_source.sourceCrs())
spatialIndex = QgsSpatialIndex(point_source.getFeatures(
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs(), context.transformContext())), feedback)
point_attribute_indices = []
if weight_field_index >= 0:
point_attribute_indices.append(weight_field_index)
if class_field_index >= 0:
point_attribute_indices.append(class_field_index)
features = poly_source.getFeatures()
total = 100.0 / poly_source.featureCount() if poly_source.featureCount() else 0
for current, polygon_feature in enumerate(features):
if feedback.isCanceled():
break
count = 0
output_feature = QgsFeature()
if polygon_feature.hasGeometry():
geom = polygon_feature.geometry()
engine = QgsGeometry.createGeometryEngine(geom.constGet())
engine.prepareGeometry()
count = 0
classes = set()
points = spatialIndex.intersects(geom.boundingBox())
if len(points) > 0:
request = QgsFeatureRequest().setFilterFids(points).setDestinationCrs(poly_source.sourceCrs(), context.transformContext())
request.setSubsetOfAttributes(point_attribute_indices)
for point_feature in point_source.getFeatures(request):
if feedback.isCanceled():
break
if engine.contains(point_feature.geometry().constGet()):
if weight_field_index >= 0:
weight = point_feature.attributes()[weight_field_index]
try:
count += float(weight)
except:
# Ignore fields with non-numeric values
pass
elif class_field_index >= 0:
point_class = point_feature.attributes()[class_field_index]
if point_class not in classes:
classes.add(point_class)
else:
count += 1
output_feature.setGeometry(geom)
attrs = polygon_feature.attributes()
if class_field_index >= 0:
score = len(classes)
else:
score = count
if field_index == len(attrs):
attrs.append(score)
else:
attrs[field_index] = score
output_feature.setAttributes(attrs)
sink.addFeature(output_feature, QgsFeatureSink.FastInsert)
feedback.setProgress(int(current * total))
return {self.OUTPUT: dest_id}
示例12: processAlgorithm
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def processAlgorithm(self, parameters, context, feedback):
sourceA = self.parameterAsSource(parameters, self.INPUT, context)
sourceB = self.parameterAsSource(parameters, self.INTERSECT, context)
fieldsA = self.parameterAsFields(parameters, self.INPUT_FIELDS, context)
fieldsB = self.parameterAsFields(parameters, self.INTERSECT_FIELDS, context)
fieldListA = QgsFields()
field_indices_a = []
if len(fieldsA) > 0:
for f in fieldsA:
idxA = sourceA.fields().lookupField(f)
if idxA >= 0:
field_indices_a.append(idxA)
fieldListA.append(sourceA.fields()[idxA])
else:
fieldListA = sourceA.fields()
field_indices_a = [i for i in range(0, fieldListA.count())]
fieldListB = QgsFields()
field_indices_b = []
if len(fieldsB) > 0:
for f in fieldsB:
idxB = sourceB.fields().lookupField(f)
if idxB >= 0:
field_indices_b.append(idxB)
fieldListB.append(sourceB.fields()[idxB])
else:
fieldListB = sourceB.fields()
field_indices_b = [i for i in range(0, fieldListB.count())]
fieldListB = vector.testForUniqueness(fieldListA, fieldListB)
for b in fieldListB:
fieldListA.append(b)
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fieldListA, QgsWkbTypes.Point, sourceA.sourceCrs())
spatialIndex = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())), feedback)
outFeat = QgsFeature()
features = sourceA.getFeatures(QgsFeatureRequest().setSubsetOfAttributes(field_indices_a))
total = 100.0 / sourceA.featureCount() if sourceA.featureCount() else 0
for current, inFeatA in enumerate(features):
if feedback.isCanceled():
break
if not inFeatA.hasGeometry():
continue
inGeom = inFeatA.geometry()
has_intersections = False
lines = spatialIndex.intersects(inGeom.boundingBox())
engine = None
if len(lines) > 0:
has_intersections = True
# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
engine.prepareGeometry()
if has_intersections:
request = QgsFeatureRequest().setFilterFids(lines)
request.setDestinationCrs(sourceA.sourceCrs())
request.setSubsetOfAttributes(field_indices_b)
for inFeatB in sourceB.getFeatures(request):
if feedback.isCanceled():
break
tmpGeom = inFeatB.geometry()
points = []
if engine.intersects(tmpGeom.geometry()):
tempGeom = inGeom.intersection(tmpGeom)
out_attributes = [inFeatA.attributes()[i] for i in field_indices_a]
out_attributes.extend([inFeatB.attributes()[i] for i in field_indices_b])
if tempGeom.type() == QgsWkbTypes.PointGeometry:
if tempGeom.isMultipart():
points = tempGeom.asMultiPoint()
else:
points.append(tempGeom.asPoint())
for j in points:
outFeat.setGeometry(tempGeom.fromPoint(j))
outFeat.setAttributes(out_attributes)
sink.addFeature(outFeat, QgsFeatureSink.FastInsert)
feedback.setProgress(int(current * total))
return {self.OUTPUT: dest_id}
示例13: __init__
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
class TriangleMesh:
# 0 - 3
# | / |
# 1 - 2
def __init__(self, xmin, ymin, xmax, ymax, x_segments, y_segments):
self.vbands = []
self.hbands = []
self.vidx = QgsSpatialIndex()
self.hidx = QgsSpatialIndex()
xres = (xmax - xmin) / x_segments
yres = (ymax - ymin) / y_segments
self.xmin, self.ymax, self.xres, self.yres = xmin, ymax, xres, yres
def addVBand(idx, geom):
f = QgsFeature(idx)
f.setGeometry(geom)
self.vbands.append(f)
self.vidx.insertFeature(f)
def addHBand(idx, geom):
f = QgsFeature(idx)
f.setGeometry(geom)
self.hbands.append(f)
self.hidx.insertFeature(f)
for x in range(x_segments):
addVBand(x, QgsGeometry.fromRect(QgsRectangle(xmin + x * xres, ymin, xmin + (x + 1) * xres, ymax)))
for y in range(y_segments):
addHBand(y, QgsGeometry.fromRect(QgsRectangle(xmin, ymax - (y + 1) * yres, xmax, ymax - y * yres)))
def vSplit(self, geom):
"""split polygon vertically"""
for idx in self.vidx.intersects(geom.boundingBox()):
yield idx, geom.intersection(self.vbands[idx].geometry())
def hIntersects(self, geom):
"""indices of horizontal bands that intersect with geom"""
for idx in self.hidx.intersects(geom.boundingBox()):
if geom.intersects(self.hbands[idx].geometry()):
yield idx
def splitPolygons(self, geom):
xmin, ymax, xres, yres = self.xmin, self.ymax, self.xres, self.yres
for x, vi in self.vSplit(geom):
for y in self.hIntersects(vi):
pt0 = QgsPoint(xmin + x * xres, ymax - y * yres)
pt1 = QgsPoint(xmin + x * xres, ymax - (y + 1) * yres)
pt2 = QgsPoint(xmin + (x + 1) * xres, ymax - (y + 1) * yres)
pt3 = QgsPoint(xmin + (x + 1) * xres, ymax - y * yres)
quad = QgsGeometry.fromPolygon([[pt0, pt1, pt2, pt3, pt0]])
tris = [[[pt0, pt1, pt3, pt0]], [[pt3, pt1, pt2, pt3]]]
if geom.contains(quad):
yield tris[0]
yield tris[1]
else:
for i, tri in enumerate(map(QgsGeometry.fromPolygon, tris)):
if geom.contains(tri):
yield tris[i]
elif geom.intersects(tri):
poly = geom.intersection(tri)
if poly.isMultipart():
for sp in poly.asMultiPolygon():
yield sp
else:
yield poly.asPolygon()
示例14: spaced
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
def spaced(bar,buildings_layer_path,receiver_points_layer_path,spaced_pts_distance):
distance_from_facades = 0.1
buildings_layer_name = os.path.splitext(os.path.basename(buildings_layer_path))[0]
buildings_layer = QgsVectorLayer(buildings_layer_path,buildings_layer_name,"ogr")
# cp building layer to delete all fields
buildings_memory_layer = QgsVectorLayer("Polygon?crs=" + str(buildings_layer.crs().authid()), "polygon_memory_layer", "memory")
buildings_memory_layer.dataProvider().addAttributes([])
buildings_feat_all = buildings_layer.dataProvider().getFeatures()
buildings_feat_list = []
for buildings_feat in buildings_feat_all:
buildings_feat_list.append(buildings_feat)
buildings_memory_layer.dataProvider().addFeatures(buildings_feat_list)
buildings_memory_layer.updateExtents()
# this is crazy: I had to addd this line otherwise the first processing doesn't work...
QgsProject.instance().addMapLayers([buildings_memory_layer])
bar.setValue(1)
# this processing alg has as output['OUTPUT'] the layer
output = processing.run("native:buffer", {'INPUT': buildings_memory_layer,
'DISTANCE': distance_from_facades,
'DISSOLVE': False,
'OUTPUT': 'memory:'})
# I can now remove the layer from map...
QgsProject.instance().removeMapLayers( [buildings_memory_layer.id()] )
bar.setValue(25)
# this processing alg has as output['OUTPUT'] the layer
output = processing.run("qgis:polygonstolines", {'INPUT': output['OUTPUT'],
'OUTPUT': 'memory:'})
bar.setValue(50)
# this processing alg has as output['output'] the layer path...
poly_to_lines = output['OUTPUT']
output = processing.run("qgis:pointsalonglines", {'INPUT': poly_to_lines,
'DISTANCE': spaced_pts_distance,
'START_OFFSET': 0,
'END_OFFSET': 0,
'OUTPUT': 'memory:'})
bar.setValue(75)
receiver_points_memory_layer = output['OUTPUT']
del output
## Delete pts in buildings
# creates SpatialIndex
buildings_feat_all = buildings_layer.dataProvider().getFeatures()
buildings_spIndex = QgsSpatialIndex()
buildings_feat_all_dict = {}
for buildings_feat in buildings_feat_all:
buildings_spIndex.insertFeature(buildings_feat)
buildings_feat_all_dict[buildings_feat.id()] = buildings_feat
receiver_points_memory_layer_all = receiver_points_memory_layer.dataProvider().getFeatures()
receiver_points_layer_fields = QgsFields()
receiver_points_layer_fields.append(QgsField("id_pt", QVariant.Int))
receiver_points_layer_fields.append(QgsField("id_bui", QVariant.Int))
receiver_points_layer_writer = QgsVectorFileWriter(receiver_points_layer_path, "System",
receiver_points_layer_fields, QgsWkbTypes.Point,
buildings_layer.crs(), "ESRI Shapefile")
receiver_points_feat_id = 0
receiver_memory_feat_total = receiver_points_memory_layer.dataProvider().featureCount()
receiver_memory_feat_number = 0
for receiver_memory_feat in receiver_points_memory_layer_all:
receiver_memory_feat_number = receiver_memory_feat_number + 1
barValue = receiver_memory_feat_number/float(receiver_memory_feat_total)*25 + 75
bar.setValue(barValue)
rect = QgsRectangle()
rect.setXMinimum(receiver_memory_feat.geometry().asPoint().x() - distance_from_facades)
rect.setXMaximum(receiver_memory_feat.geometry().asPoint().x() + distance_from_facades)
rect.setYMinimum(receiver_memory_feat.geometry().asPoint().y() - distance_from_facades)
rect.setYMaximum(receiver_memory_feat.geometry().asPoint().y() + distance_from_facades)
buildings_selection = buildings_spIndex.intersects(rect)
to_add = True
receiver_geom = receiver_memory_feat.geometry()
building_id_correct = None
for buildings_id in buildings_selection:
#.........这里部分代码省略.........
示例15: run
# 需要导入模块: from qgis.core import QgsSpatialIndex [as 别名]
# 或者: from qgis.core.QgsSpatialIndex import intersects [as 别名]
#.........这里部分代码省略.........
has_hazard_objects = True
if not has_hazard_objects:
message = tr(
'There are no objects in the hazard layer with %s '
'value in %s. Please check your data or use another '
'attribute.') % (
self.hazard_class_attribute,
', '.join(self.hazard_class_mapping[self.wet]))
raise GetDataError(message)
# Filter out just those EXPOSURE features in the analysis extents
transform = QgsCoordinateTransform(
self.requested_extent_crs, self.exposure.layer.crs())
projected_extent = transform.transformBoundingBox(requested_extent)
request = QgsFeatureRequest()
request.setFilterRect(projected_extent)
# We will use this transform to project each exposure feature into
# the CRS of the Hazard.
transform = QgsCoordinateTransform(
self.exposure.crs(), self.hazard.crs())
features = []
for feature in self.exposure.layer.getFeatures(request):
# Make a deep copy as the geometry is passed by reference
# If we don't do this, subsequent operations will affect the
# original feature geometry as well as the copy TS
building_geom = QgsGeometry(feature.geometry())
# Project the building geometry to hazard CRS
building_bounds = transform.transform(building_geom.boundingBox())
affected = False
# get tentative list of intersecting hazard features
# only based on intersection of bounding boxes
ids = hazard_index.intersects(building_bounds)
for fid in ids:
# run (slow) exact intersection test
building_geom.transform(transform)
if hazard_geometries[fid].intersects(building_geom):
affected = True
break
new_feature = QgsFeature()
# We write out the original feature geom, not the projected one
new_feature.setGeometry(feature.geometry())
new_feature.setAttributes(feature.attributes())
new_feature[target_field_index] = 1 if affected else 0
features.append(new_feature)
# every once in a while commit the created features
# to the output layer
if len(features) == 1000:
(_, __) = building_provider.addFeatures(features)
features = []
(_, __) = building_provider.addFeatures(features)
building_layer.updateExtents()
# Generate simple impact report
self.buildings = {}
self.affected_buildings = OrderedDict([
(tr('Flooded'), {})
])
buildings_data = building_layer.getFeatures()
building_type_field_index = building_layer.fieldNameIndex(
self.exposure_class_attribute)
for building in buildings_data:
record = building.attributes()