Python QgsSpatialIndex.intersects方法代码示例

# 需要导入模块: 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

# 需要导入模块: 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}

# 需要导入模块: 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

# 需要导入模块: 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()

# 需要导入模块: 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

# 需要导入模块: 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

# 需要导入模块: 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}

# 需要导入模块: 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)

# 需要导入模块: 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
#.........这里部分代码省略.........

# 需要导入模块: 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}

# 需要导入模块: 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}

# 需要导入模块: 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}

# 需要导入模块: 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()

# 需要导入模块: 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:
#.........这里部分代码省略.........

# 需要导入模块: 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()