Python QgsGeometry.intersection方法代码示例

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT))
        vlayerB = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT2))
        vproviderA = vlayerA.dataProvider()

        fields = vector.combineVectorFields(vlayerA, vlayerB)
        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
            fields, vproviderA.geometryType(), vproviderA.crs()
        )
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        index = vector.spatialindex(vlayerB)
        nElement = 0
        selectionA = vector.features(vlayerA)
        nFeat = len(selectionA)
        for inFeatA in selectionA:
            nElement += 1
            progress.setPercentage(nElement / float(nFeat) * 100)
            geom = QgsGeometry(inFeatA.geometry())
            atMapA = inFeatA.attributes()
            intersects = index.intersects(geom.boundingBox())
            for i in intersects:
                request = QgsFeatureRequest().setFilterFid(i)
                inFeatB = vlayerB.getFeatures(request).next()
                tmpGeom = QgsGeometry(inFeatB.geometry())
                try:
                    if geom.intersects(tmpGeom):
                        atMapB = inFeatB.attributes()
                        int_geom = QgsGeometry(geom.intersection(tmpGeom))
                        if (
                            int_geom.wkbType() == QGis.WKBUnknown
                            or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection
                        ):
                            int_com = geom.combine(tmpGeom)
                            int_sym = geom.symDifference(tmpGeom)
                            int_geom = QgsGeometry(int_com.difference(int_sym))
                        try:
                            if int_geom.wkbType() in wkbTypeGroups[wkbTypeGroups[int_geom.wkbType()]]:
                                outFeat.setGeometry(int_geom)
                                attrs = []
                                attrs.extend(atMapA)
                                attrs.extend(atMapB)
                                outFeat.setAttributes(attrs)
                                writer.addFeature(outFeat)
                        except:
                            ProcessingLog.addToLog(
                                ProcessingLog.LOG_INFO,
                                self.tr(
                                    "Feature geometry error: One or more output features ignored due to invalid geometry."
                                ),
                            )
                            continue
                except:
                    break

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = dataobjects.getObjectFromUri(
            self.getParameterValue(self.INPUT))
        vlayerB = dataobjects.getObjectFromUri(
            self.getParameterValue(self.INPUT2))
        vproviderA = vlayerA.dataProvider()

        geomType = vproviderA.geometryType()
        if geomType in GEOM_25D:
            raise GeoAlgorithmExecutionException(
                self.tr('Input layer has unsupported geometry type {}').format(geomType))

        fields = vector.combineVectorFields(vlayerA, vlayerB)
        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fields,
                                                                     geomType, vproviderA.crs())
        outFeat = QgsFeature()
        index = vector.spatialindex(vlayerB)
        selectionA = vector.features(vlayerA)
        total = 100.0 / len(selectionA)
        for current, inFeatA in enumerate(selectionA):
            progress.setPercentage(int(current * total))
            geom = QgsGeometry(inFeatA.geometry())
            atMapA = inFeatA.attributes()
            intersects = index.intersects(geom.boundingBox())
            for i in intersects:
                request = QgsFeatureRequest().setFilterFid(i)
                inFeatB = vlayerB.getFeatures(request).next()
                tmpGeom = QgsGeometry(inFeatB.geometry())
                if geom.intersects(tmpGeom):
                    atMapB = inFeatB.attributes()
                    int_geom = QgsGeometry(geom.intersection(tmpGeom))
                    if int_geom.wkbType() == QGis.WKBUnknown or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                        int_com = geom.combine(tmpGeom)
                        int_sym = geom.symDifference(tmpGeom)
                        int_geom = QgsGeometry(int_com.difference(int_sym))
                    if int_geom.isGeosEmpty() or not int_geom.isGeosValid():
                        ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                                               self.tr('GEOS geoprocessing error: One or '
                                                       'more input features have invalid '
                                                       'geometry.'))
                        break
                    try:
                        if int_geom.wkbType() in wkbTypeGroups[wkbTypeGroups[int_geom.wkbType()]]:
                            outFeat.setGeometry(int_geom)
                            attrs = []
                            attrs.extend(atMapA)
                            attrs.extend(atMapB)
                            outFeat.setAttributes(attrs)
                            writer.addFeature(outFeat)
                    except:
                        ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                               self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
                        continue

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
 def compute(self, inPoly, inLns, inField, outPath, progressBar):
     polyLayer = ftools_utils.getVectorLayerByName(inPoly)
     lineLayer = ftools_utils.getVectorLayerByName(inLns)
     polyProvider = polyLayer.dataProvider()
     lineProvider = lineLayer.dataProvider()
     if polyProvider.crs() != lineProvider.crs():
         QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
     fieldList = ftools_utils.getFieldList(polyLayer)
     index = polyProvider.fieldNameIndex(unicode(inField))
     if index == -1:
         index = polyProvider.fields().count()
         fieldList.append(QgsField(unicode(inField), QVariant.Double, "real", 24, 15, self.tr("length field")))
     sRs = polyProvider.crs()
     inFeat = QgsFeature()
     inFeatB = QgsFeature()
     outFeat = QgsFeature()
     inGeom = QgsGeometry()
     outGeom = QgsGeometry()
     distArea = QgsDistanceArea()
     start = 0.00
     add = 100.00 / polyProvider.featureCount()
     check = QFile(self.shapefileName)
     if check.exists():
         if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
             return
     writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, polyProvider.geometryType(), sRs)
     spatialIndex = ftools_utils.createIndex(lineProvider)
     polyFit = polyProvider.getFeatures()
     while polyFit.nextFeature(inFeat):
         inGeom = QgsGeometry(inFeat.geometry())
         atMap = inFeat.attributes()
         lineList = []
         length = 0
         lineList = spatialIndex.intersects(inGeom.boundingBox())
         if len(lineList) > 0:
             check = 0
         else:
             check = 1
         if check == 0:
             for i in lineList:
                 lineProvider.getFeatures(QgsFeatureRequest().setFilterFid(int(i))).nextFeature(inFeatB)
                 tmpGeom = QgsGeometry(inFeatB.geometry())
                 if inGeom.intersects(tmpGeom):
                     outGeom = inGeom.intersection(tmpGeom)
                     length = length + distArea.measure(outGeom)
         outFeat.setGeometry(inGeom)
         atMap.append(length)
         outFeat.setAttributes(atMap)
         writer.addFeature(outFeat)
         start = start + 1
         progressBar.setValue(start * (add))
     del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
 def testIntersection(self):
     myLine = QgsGeometry.fromPolyline([QgsPoint(0, 0),QgsPoint(1, 1),QgsPoint(2, 2)])
     myPoint = QgsGeometry.fromPoint(QgsPoint(1, 1))
     intersectionGeom = QgsGeometry.intersection(myLine, myPoint)
     myMessage = ('Expected:\n%s\nGot:\n%s\n' %
                   (QGis.Point, intersectionGeom.type()))
     assert intersectionGeom.wkbType() == QGis.WKBPoint, myMessage
     
     layer = QgsVectorLayer("Point", "intersection", "memory")
     assert layer.isValid(), "Failed to create valid point memory layer"
     
     provider = layer.dataProvider()
             
     ft = QgsFeature()
     ft.setGeometry(intersectionGeom)
     provider.addFeatures([ft])
     
     myMessage = ('Expected:\n%s\nGot:\n%s\n' %
                   (1, layer.featureCount()))
     assert layer.featureCount() == 1, myMessage

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        lineLayer = dataobjects.getObjectFromUri(self.getParameterValue(self.LINES))
        polyLayer = dataobjects.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        lengthFieldName = self.getParameterValue(self.LEN_FIELD)
        countFieldName = self.getParameterValue(self.COUNT_FIELD)

        (idxLength, fieldList) = vector.findOrCreateField(polyLayer,
                                                          polyLayer.fields(), lengthFieldName)
        (idxCount, fieldList) = vector.findOrCreateField(polyLayer, fieldList,
                                                         countFieldName)

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
            fieldList.toList(), polyLayer.wkbType(), polyLayer.crs())

        spatialIndex = vector.spatialindex(lineLayer)

        ftLine = QgsFeature()
        ftPoly = QgsFeature()
        outFeat = QgsFeature()
        inGeom = QgsGeometry()
        outGeom = QgsGeometry()
        distArea = QgsDistanceArea()

        features = vector.features(polyLayer)
        total = 100.0 / len(features)
        hasIntersections = False
        for current, ftPoly in enumerate(features):
            inGeom = ftPoly.geometry()
            attrs = ftPoly.attributes()
            count = 0
            length = 0
            hasIntersections = False
            lines = spatialIndex.intersects(inGeom.boundingBox())
            engine = None
            if len(lines) > 0:
                hasIntersections = True
                # use prepared geometries for faster intersection tests
                engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
                engine.prepareGeometry()

            if hasIntersections:
                request = QgsFeatureRequest().setFilterFids(lines).setSubsetOfAttributes([])
                for ftLine in lineLayer.getFeatures(request):
                    tmpGeom = ftLine.geometry()
                    if engine.intersects(tmpGeom.geometry()):
                        outGeom = inGeom.intersection(tmpGeom)
                        length += distArea.measureLength(outGeom)
                        count += 1

            outFeat.setGeometry(inGeom)
            if idxLength == len(attrs):
                attrs.append(length)
            else:
                attrs[idxLength] = length
            if idxCount == len(attrs):
                attrs.append(count)
            else:
                attrs[idxCount] = count
            outFeat.setAttributes(attrs)
            writer.addFeature(outFeat)

            progress.setPercentage(int(current * total))

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        lineLayer = dataobjects.getObjectFromUri(self.getParameterValue(self.LINES))
        polyLayer = dataobjects.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        lengthFieldName = self.getParameterValue(self.LEN_FIELD)
        countFieldName = self.getParameterValue(self.COUNT_FIELD)

        polyProvider = polyLayer.dataProvider()

        (idxLength, fieldList) = vector.findOrCreateField(polyLayer,
                polyLayer.pendingFields(), lengthFieldName)
        (idxCount, fieldList) = vector.findOrCreateField(polyLayer, fieldList,
                countFieldName)

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
            fieldList.toList(), polyProvider.geometryType(), polyProvider.crs())

        spatialIndex = vector.spatialindex(lineLayer)

        ftLine = QgsFeature()
        ftPoly = QgsFeature()
        outFeat = QgsFeature()
        inGeom = QgsGeometry()
        outGeom = QgsGeometry()
        distArea = QgsDistanceArea()

        current = 0
        features = vector.features(polyLayer)
        total = 100.0 / float(len(features))
        hasIntersections = False
        for ftPoly in features:
            inGeom = QgsGeometry(ftPoly.geometry())
            attrs = ftPoly.attributes()
            count = 0
            length = 0
            hasIntersections = False
            lines = spatialIndex.intersects(inGeom.boundingBox())
            if len(lines) > 0:
                hasIntersections = True

            if hasIntersections:
                for i in lines:
                    request = QgsFeatureRequest().setFilterFid(i)
                    ftLine = lineLayer.getFeatures(request).next()
                    tmpGeom = QgsGeometry(ftLine.geometry())
                    if inGeom.intersects(tmpGeom):
                        outGeom = inGeom.intersection(tmpGeom)
                        length += distArea.measure(outGeom)
                        count += 1

            outFeat.setGeometry(inGeom)
            if idxLength == len(attrs):
                attrs.append(length)
            else:
                attrs[idxLength] = length
            if idxCount == len(attrs):
                attrs.append(count)
            else:
                attrs[idxCount] = count
            outFeat.setAttributes(attrs)
            writer.addFeature(outFeat)

            current += 1
            progress.setPercentage(int(current * total))

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        layerA = dataobjects.getObjectFromUri(
            self.getParameterValue(Clip.INPUT))
        layerB = dataobjects.getObjectFromUri(
            self.getParameterValue(Clip.OVERLAY))

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
            layerA.pendingFields(),
            layerA.dataProvider().geometryType(),
            layerA.dataProvider().crs())

        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()

        index = vector.spatialindex(layerB)

        selectionA = vector.features(layerA)

        current = 0
        total = 100.0 / float(len(selectionA))

        for inFeatA in selectionA:
            geom = QgsGeometry(inFeatA.geometry())
            attrs = inFeatA.attributes()
            intersects = index.intersects(geom.boundingBox())
            first = True
            found = False
            if len(intersects) > 0:
                for i in intersects:
                    layerB.getFeatures(
                        QgsFeatureRequest().setFilterFid(i)).nextFeature(
                            inFeatB)
                    tmpGeom = QgsGeometry(inFeatB.geometry())
                    if tmpGeom.intersects(geom):
                        found = True
                        if first:
                            outFeat.setGeometry(QgsGeometry(tmpGeom))
                            first = False
                        else:
                            try:
                                cur_geom = QgsGeometry(outFeat.geometry())
                                new_geom = QgsGeometry(
                                    cur_geom.combine(tmpGeom))
                                outFeat.setGeometry(QgsGeometry(new_geom))
                            except:
                                ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                                                       self.tr('GEOS geoprocessing error: One or '
                                                               'more input features have invalid '
                                                               'geometry.'))
                                break
                if found:
                    try:
                        cur_geom = QgsGeometry(outFeat.geometry())
                        new_geom = QgsGeometry(geom.intersection(cur_geom))
                        if new_geom.wkbType() == QGis.WKBUnknown or QgsWKBTypes.flatType(new_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                            int_com = QgsGeometry(geom.combine(cur_geom))
                            int_sym = QgsGeometry(geom.symDifference(cur_geom))
                            new_geom = QgsGeometry(int_com.difference(int_sym))
                        try:
                            outFeat.setGeometry(new_geom)
                            outFeat.setAttributes(attrs)
                            writer.addFeature(outFeat)
                        except:
                            ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                                                   self.tr('Feature geometry error: One or more '
                                                           'output features ignored due to '
                                                           'invalid geometry.'))
                            continue
                    except:
                        ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                                               self.tr('GEOS geoprocessing error: One or more '
                                                       'input features have invalid geometry.'))
                        continue

            current += 1
            progress.setPercentage(int(current * total))

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT))
        vlayerB = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT2))

        geomType = vlayerA.wkbType()
        fields = vector.combineVectorFields(vlayerA, vlayerB)
        writer = self.getOutputFromName(Union.OUTPUT).getVectorWriter(fields,
                                                                      geomType,
                                                                      vlayerA.crs())
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        indexA = vector.spatialindex(vlayerB)
        indexB = vector.spatialindex(vlayerA)

        count = 0
        nElement = 0
        featuresA = vector.features(vlayerA)
        nFeat = len(featuresA) if len(featuresA) > 0 else 1
        for inFeatA in featuresA:
            progress.setPercentage(nElement / float(nFeat) * 50)
            nElement += 1
            lstIntersectingB = []
            geom = QgsGeometry(inFeatA.geometry())
            atMapA = inFeatA.attributes()
            intersects = indexA.intersects(geom.boundingBox())
            if len(intersects) < 1:
                try:
                    outFeat.setGeometry(geom)
                    outFeat.setAttributes(atMapA)
                    writer.addFeature(outFeat)
                except:
                    # This really shouldn't happen, as we haven't
                    # edited the input geom at all
                    ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                           self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
            else:
                for id in intersects:
                    count += 1
                    request = QgsFeatureRequest().setFilterFid(id)
                    inFeatB = vlayerB.getFeatures(request).next()
                    atMapB = inFeatB.attributes()
                    tmpGeom = QgsGeometry(inFeatB.geometry())

                    if geom.intersects(tmpGeom):
                        int_geom = geom.intersection(tmpGeom)
                        lstIntersectingB.append(tmpGeom)

                        if int_geom is None:
                            # There was a problem creating the intersection
                            ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                   self.tr('GEOS geoprocessing error: One or more input features have invalid geometry.'))
                            int_geom = QgsGeometry()
                        else:
                            int_geom = QgsGeometry(int_geom)

                        # TODO: the result may have a different dimension (e.g. intersection of two polygons may result in a single point)
                        # or the result may be a collection of geometries (e.g. intersection of two polygons results in three polygons and one linestring).
                        # We need to filter out all acceptable geometries into a single (possibly multi-part) geometry - and we need
                        # to do it consistently also in the code further below

                        if int_geom.wkbType() == QGis.WKBUnknown or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                            # Intersection produced different geomety types
                            temp_list = int_geom.asGeometryCollection()
                            for i in temp_list:
                                if i.type() == geom.type():
                                    int_geom = QgsGeometry(i)
                                    try:
                                        outFeat.setGeometry(int_geom)
                                        outFeat.setAttributes(atMapA + atMapB)
                                        writer.addFeature(outFeat)
                                    except:
                                        ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                               self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
                        else:
                            # Geometry list: prevents writing error
                            # in geometries of different types
                            # produced by the intersection
                            # fix #3549
                            if int_geom.wkbType() in wkbTypeGroups[wkbTypeGroups[int_geom.wkbType()]]:
                                try:
                                    outFeat.setGeometry(int_geom)
                                    outFeat.setAttributes(atMapA + atMapB)
                                    writer.addFeature(outFeat)
                                except:
                                    ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                           self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))

                # the remaining bit of inFeatA's geometry
                # if there is nothing left, this will just silently fail and we're good
                diff_geom = QgsGeometry(geom)
                if len(lstIntersectingB) != 0:
                    intB = QgsGeometry.unaryUnion(lstIntersectingB)
                    diff_geom = diff_geom.difference(intB)
                    if diff_geom is None:
                        ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                               self.tr('GEOS geoprocessing error: One or more input features have invalid geometry.'))
                        diff_geom = QgsGeometry()
                    if diff_geom.isGeosEmpty() or not diff_geom.isGeosValid():
                        ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]

#.........这里部分代码省略.........
            start = 20.00
            add = 80.00 / len(featToEliminate)
        else:
            start = 100

        progress.setPercentage(start)
        madeProgress = True

        # We go through the list and see if we find any polygons we can
        # merge the selected with. If we have no success with some we
        # merge and then restart the whole story.
        while madeProgress:  # Check if we made any progress
            madeProgress = False
            featNotEliminated = []

            # Iterate over the polygons to eliminate
            for i in range(len(featToEliminate)):
                feat = featToEliminate.pop()
                geom2Eliminate = QgsGeometry(feat.geometry())
                bbox = geom2Eliminate.boundingBox()
                fit = processLayer.getFeatures(
                    QgsFeatureRequest().setFilterRect(bbox))
                mergeWithFid = None
                mergeWithGeom = None
                max = 0
                min = -1
                selFeat = QgsFeature()

                while fit.nextFeature(selFeat):
                    selGeom = QgsGeometry(selFeat.geometry())

                    if geom2Eliminate.intersects(selGeom):
                        # We have a candidate
                        iGeom = geom2Eliminate.intersection(selGeom)

                        if iGeom is None:
                            continue

                        if boundary:
                            selValue = iGeom.length()
                        else:
                            # area. We need a common boundary in
                            # order to merge
                            if 0 < iGeom.length():
                                selValue = selGeom.area()
                            else:
                                selValue = -1

                        if -1 != selValue:
                            useThis = True

                            if smallestArea:
                                if -1 == min:
                                    min = selValue
                                else:
                                    if selValue < min:
                                        min = selValue
                                    else:
                                        useThis = False
                            else:
                                if selValue > max:
                                    max = selValue
                                else:
                                    useThis = False

                            if useThis:

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        layerA = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT_A))
        layerB = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT_B))
        fieldA = self.getParameterValue(self.FIELD_A)
        fieldB = self.getParameterValue(self.FIELD_B)

        idxA = layerA.fieldNameIndex(fieldA)
        idxB = layerB.fieldNameIndex(fieldB)

        fieldList = [layerA.pendingFields()[idxA],
                     layerB.pendingFields()[idxB]]

        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fieldList,
                QGis.WKBPoint, layerA.dataProvider().crs())

        spatialIndex = vector.spatialindex(layerB)

        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        inGeom = QgsGeometry()
        tmpGeom = QgsGeometry()

        features = vector.features(layerA)

        current = 0
        total = 100.0 / float(len(features))
        hasIntersections = False

        for inFeatA in features:
            inGeom = inFeatA.geometry()
            hasIntersections = False
            lines = spatialIndex.intersects(inGeom.boundingBox())

            if len(lines) > 0:
                hasIntersections = True

            if hasIntersections:
                for i in lines:
                    request = QgsFeatureRequest().setFilterFid(i)
                    inFeatB = layerB.getFeatures(request).next()
                    tmpGeom = QgsGeometry(inFeatB.geometry())

                    points = []
                    attrsA = inFeatA.attributes()
                    attrsB = inFeatB.attributes()

                    if inGeom.intersects(tmpGeom):
                        tempGeom = inGeom.intersection(tmpGeom)
                        if tempGeom.type() == QGis.Point:
                            if tempGeom.isMultipart():
                                points = tempGeom.asMultiPoint()
                            else:
                                points.append(tempGeom.asPoint())

                            for j in points:
                                outFeat.setGeometry(tempGeom.fromPoint(j))
                                outFeat.setAttributes([attrsA[idxA],
                                        attrsB[idxB]])
                                writer.addFeature(outFeat)

            current += 1
            progress.setPercentage(int(current * total))

        del writer

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def intersect(featureId, feature1, feature2, mousePoint):
        """
        To check if there is an intersection between 2 features close to a given point
        :param featureId: if we want to snap on a given feature
        :param feature1: the first feature
        :param feature2: the second feature
        :param mousePoint: the given point
        :return: the intersection as QgsPoint or none
        """

        if featureId is None or feature1.id() == featureId or feature2.id() == featureId:
            geometry1 = feature1.geometry()
            geometry2 = feature2.geometry()
            if geometry1.type() == 0:
                return geometry1.asPoint()
            if geometry2.type() == 0:
                return geometry2.asPoint()
            if geometry1.type() == 2:
                polygon = geometry1.geometry()
                newG = polygon.boundary()
                geometry1 = QgsGeometry(newG)
            if geometry2.type() == 2:
                polygon = geometry2.geometry()
                newG = polygon.boundary()
                geometry2 = QgsGeometry(newG)

            intersection = geometry1.intersection(geometry2)
            if intersection.type() == 0:
                intersectionP = intersection.asPoint()
                intersectionMP = intersection.asMultiPoint()
                if intersectionMP is not None and len(intersectionMP) > 0:
                    for point in intersectionMP:
                        if intersectionP is None:
                            intersectionP = point
                        elif mousePoint.sqrDist(point) < mousePoint.sqrDist(intersectionP):
                            intersectionP = QgsPoint(point.x(), point.y())
            elif intersection.type() == 1:
                intersectionPL = intersection.asPolyline()
                intersectionMPL = intersection.asMultiPolyline()
                intersectionP = None
                if intersectionMPL is not None and len(intersectionMPL) > 0:
                    for line in intersectionMPL:
                        for point in line:
                            if intersectionP is None:
                                intersectionP = point
                            elif mousePoint.sqrDist(point) < mousePoint.sqrDist(intersectionP):
                                intersectionP = QgsPoint(point.x(), point.y())
                else:
                    for point in intersectionPL:
                        if intersectionP is None:
                            intersectionP = point
                        elif mousePoint.sqrDist(point) < mousePoint.sqrDist(intersectionP):
                            intersectionP = QgsPoint(point.x(), point.y())
            elif intersection.type() == 2:
                intersectionMPL = intersection.asMultiPolyline()
                intersectionP = None
                for line in intersectionMPL:
                    for point in line:
                        if intersectionP is None:
                            intersectionP = point
                        elif mousePoint.sqrDist(point) < mousePoint.sqrDist(intersectionP):
                            intersectionP = QgsPoint(point.x(), point.y())
            else:
                return None

            if intersectionP and intersectionP != QgsPoint(0, 0):
                return intersectionP

        return None

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT))
        vlayerB = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT2))
        GEOS_EXCEPT = True
        FEATURE_EXCEPT = True
        vproviderA = vlayerA.dataProvider()

        fields = vector.combineVectorFields(vlayerA, vlayerB)
        names = [field.name() for field in fields]
        ProcessingLog.addToLog(ProcessingLog.LOG_INFO, unicode(names))
        writer = self.getOutputFromName(Union.OUTPUT).getVectorWriter(fields,
                                                                      vproviderA.geometryType(), vproviderA.crs())
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        indexA = vector.spatialindex(vlayerB)
        indexB = vector.spatialindex(vlayerA)

        count = 0
        nElement = 0
        featuresA = vector.features(vlayerA)
        nFeat = len(featuresA)
        for inFeatA in featuresA:
            progress.setPercentage(nElement / float(nFeat) * 50)
            nElement += 1
            lstIntersectingB = []
            geom = QgsGeometry(inFeatA.geometry())
            atMapA = inFeatA.attributes()
            intersects = indexA.intersects(geom.boundingBox())
            if len(intersects) < 1:
                try:
                    outFeat.setGeometry(geom)
                    outFeat.setAttributes(atMapA)
                    writer.addFeature(outFeat)
                except:
                    # This really shouldn't happen, as we haven't
                    # edited the input geom at all
                    raise GeoAlgorithmExecutionException(
                        self.tr('Feature exception while computing union'))
            else:
                for id in intersects:
                    count += 1
                    request = QgsFeatureRequest().setFilterFid(id)
                    inFeatB = vlayerB.getFeatures(request).next()
                    atMapB = inFeatB.attributes()
                    tmpGeom = QgsGeometry(inFeatB.geometry())

                    if geom.intersects(tmpGeom):
                        int_geom = geom.intersection(tmpGeom)
                        lstIntersectingB.append(tmpGeom)

                        if int_geom is None:
                           # There was a problem creating the intersection
                            raise GeoAlgorithmExecutionException(
                                self.tr('Geometry exception while computing '
                                        'intersection'))
                        else:
                            int_geom = QgsGeometry(int_geom)

                        if int_geom.wkbType() == QGis.WKBUnknown or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                            # Intersection produced different geomety types
                            temp_list = int_geom.asGeometryCollection()
                            for i in temp_list:
                                if i.type() == geom.type():
                                    int_geom = QgsGeometry(i)
                        try:
                            outFeat.setGeometry(int_geom)
                            attrs = []
                            attrs.extend(atMapA)
                            attrs.extend(atMapB)
                            outFeat.setAttributes(attrs)
                            writer.addFeature(outFeat)
                        except Exception as err:
                            raise GeoAlgorithmExecutionException(
                                self.tr('Feature exception while computing union'))

                try:
                    # the remaining bit of inFeatA's geometry
                    # if there is nothing left, this will just silently fail and we're good
                    diff_geom = QgsGeometry(geom)
                    if len(lstIntersectingB) != 0:
                        intB = QgsGeometry.unaryUnion(lstIntersectingB)
                        diff_geom = diff_geom.difference(intB)

                    if diff_geom.wkbType() == 0 or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                        temp_list = diff_geom.asGeometryCollection()
                        for i in temp_list:
                            if i.type() == geom.type():
                                diff_geom = QgsGeometry(i)
                    outFeat.setGeometry(diff_geom)
                    outFeat.setAttributes(atMapA)
                    writer.addFeature(outFeat)
                except Exception as err:
                    raise GeoAlgorithmExecutionException(
                        self.tr('Feature exception while computing union'))

        length = len(vproviderA.fields())

        featuresA = vector.features(vlayerB)
        nFeat = len(featuresA)
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
  def features(self, request=None, clipGeom=None):
    settings = self.writer.settings
    mapTo3d = settings.mapTo3d()
    baseExtent = settings.baseExtent
    baseExtentGeom = baseExtent.geometry()
    rotation = baseExtent.rotation()
    prop = self.prop

    useZ = prop.useZ()
    if useZ:
      srs_from = osr.SpatialReference()
      srs_from.ImportFromProj4(str(self.layer.crs().toProj4()))
      srs_to = osr.SpatialReference()
      srs_to.ImportFromProj4(str(self.writer.settings.crs.toProj4()))

      ogr_transform = osr.CreateCoordinateTransformation(srs_from, srs_to)
      clipGeomWkb = clipGeom.asWkb() if clipGeom else None
      ogr_clipGeom = ogr.CreateGeometryFromWkb(clipGeomWkb) if clipGeomWkb else None

    else:
      # z_func: function to get elevation at given point (x, y) on surface
      if prop.isHeightRelativeToDEM():
        if self.geomType == QGis.Polygon and prop.type_index == 1:  # Overlay
          z_func = lambda x, y: 0
        else:
          # get elevation from DEM
          z_func = lambda x, y: self.writer.demProvider.readValue(x, y)
      else:
        z_func = lambda x, y: 0

    request = request or QgsFeatureRequest()
    for f in self.layer.getFeatures(request):
      geometry = f.geometry()
      if geometry is None:
        logMessage("null geometry skipped")
        continue

      # coordinate transformation - layer crs to project crs
      geom = QgsGeometry(geometry)
      if geom.transform(self.transform) != 0:
        logMessage("Failed to transform geometry")
        continue

      # check if geometry intersects with the base extent (rotated rect)
      if rotation and not baseExtentGeom.intersects(geom):
        continue

      # create feature
      feat = Feature(self.writer, self, f)

      # transform_func: function to transform the map coordinates to 3d coordinates
      relativeHeight = prop.relativeHeight(f)

      def transform_func(x, y, z):
        return mapTo3d.transform(x, y, z + relativeHeight)

      if useZ:
        ogr_geom = ogr.CreateGeometryFromWkb(geometry.asWkb())

        # transform geometry from layer CRS to project CRS
        if ogr_geom.Transform(ogr_transform) != 0:
          logMessage("Failed to transform geometry")
          continue

        # clip geometry
        if ogr_clipGeom and self.geomType == QGis.Line:
          ogr_geom = ogr_geom.Intersection(ogr_clipGeom)
          if ogr_geom is None:
            continue

        # check if geometry is empty
        if ogr_geom.IsEmpty():
          logMessage("empty geometry skipped")
          continue

        feat.geom = self.geomClass.fromOgrGeometry25D(ogr_geom, transform_func)

      else:
        # clip geometry
        if clipGeom and self.geomType in [QGis.Line, QGis.Polygon]:
          geom = geom.intersection(clipGeom)
          if geom is None:
            continue

        # check if geometry is empty
        if geom.isGeosEmpty():
          logMessage("empty geometry skipped")
          continue

        if self.geomType == QGis.Polygon:
          feat.geom = self.geomClass.fromQgsGeometry(geom, z_func, transform_func, self.hasLabel())
          if prop.type_index == 1 and prop.isHeightRelativeToDEM():   # Overlay and relative to DEM
            feat.geom.splitPolygon(self.writer.triangleMesh())

        else:
          feat.geom = self.geomClass.fromQgsGeometry(geom, z_func, transform_func)

      if feat.geom is None:
        continue

#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT))
        vlayerB = dataobjects.getObjectFromUri(self.getParameterValue(Union.INPUT2))

        vproviderA = vlayerA.dataProvider()

        geomType = vproviderA.geometryType()
        fields = vector.combineVectorFields(vlayerA, vlayerB)
        writer = self.getOutputFromName(Union.OUTPUT).getVectorWriter(fields,
                                                                      geomType, vproviderA.crs())
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        indexA = vector.spatialindex(vlayerB)
        indexB = vector.spatialindex(vlayerA)

        count = 0
        nElement = 0
        featuresA = vector.features(vlayerA)
        nFeat = len(featuresA)
        for inFeatA in featuresA:
            progress.setPercentage(nElement / float(nFeat) * 50)
            nElement += 1
            lstIntersectingB = []
            geom = QgsGeometry(inFeatA.geometry())
            atMapA = inFeatA.attributes()
            intersects = indexA.intersects(geom.boundingBox())
            if len(intersects) < 1:
                try:
                    outFeat.setGeometry(geom)
                    outFeat.setAttributes(atMapA)
                    writer.addFeature(outFeat)
                except:
                    # This really shouldn't happen, as we haven't
                    # edited the input geom at all
                    ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                           self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
            else:
                for id in intersects:
                    count += 1
                    request = QgsFeatureRequest().setFilterFid(id)
                    inFeatB = vlayerB.getFeatures(request).next()
                    atMapB = inFeatB.attributes()
                    tmpGeom = QgsGeometry(inFeatB.geometry())

                    if geom.intersects(tmpGeom):
                        int_geom = geom.intersection(tmpGeom)
                        lstIntersectingB.append(tmpGeom)

                        if int_geom is None:
                            # There was a problem creating the intersection
                            ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                   self.tr('GEOS geoprocessing error: One or more input features have invalid geometry.'))
                            int_geom = QgsGeometry()
                        else:
                            int_geom = QgsGeometry(int_geom)

                        if int_geom.wkbType() == Qgis.WKBUnknown or QgsWKBTypes.flatType(int_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                            # Intersection produced different geomety types
                            temp_list = int_geom.asGeometryCollection()
                            for i in temp_list:
                                if i.type() == geom.type():
                                    int_geom = QgsGeometry(i)
                                    try:
                                        outFeat.setGeometry(int_geom)
                                        outFeat.setAttributes(atMapA + atMapB)
                                        writer.addFeature(outFeat)
                                    except:
                                        ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                               self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
                        else:
                            # Geometry list: prevents writing error
                            # in geometries of different types
                            # produced by the intersection
                            # fix #3549
                            if int_geom.wkbType() in wkbTypeGroups[wkbTypeGroups[int_geom.wkbType()]]:
                                try:
                                    outFeat.setGeometry(int_geom)
                                    outFeat.setAttributes(atMapA + atMapB)
                                    writer.addFeature(outFeat)
                                except:
                                    ProcessingLog.addToLog(ProcessingLog.LOG_INFO,
                                                           self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))

                # the remaining bit of inFeatA's geometry
                # if there is nothing left, this will just silently fail and we're good
                diff_geom = QgsGeometry(geom)
                if len(lstIntersectingB) != 0:
                    intB = QgsGeometry.unaryUnion(lstIntersectingB)
                    diff_geom = diff_geom.difference(intB)
                    if diff_geom.isGeosEmpty() or not diff_geom.isGeosValid():
                        ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                                               self.tr('GEOS geoprocessing error: One or more input features have invalid geometry.'))

                if diff_geom.wkbType() == 0 or QgsWKBTypes.flatType(diff_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
                    temp_list = diff_geom.asGeometryCollection()
                    for i in temp_list:
                        if i.type() == geom.type():
                            diff_geom = QgsGeometry(i)
                try:
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsGeometry [as 别名]
# 或者: from qgis.core.QgsGeometry import intersection [as 别名]
    def _preparePolygonLayer(self, theQgisLayer):
        """Create a new layer with no intersecting features to self.layer.

        A helper function to align the polygons to the postprocLayer
        polygons. If one input polygon is in two or more postprocLayer polygons
        then it is divided so that each part is within only one of the
        postprocLayer polygons. this allows to aggregate in postrocessing using
        centroid in polygon.

        The function assumes EPSG:4326 but no checks are enforced

        Args:
            theQgisLayer of the file to be processed
        Returns:
            QgisLayer of the processed file

        Raises:
            Any exceptions raised by the InaSAFE library will be propagated.
        """
#        import time
#        startTime = time.clock()

        myMessage = m.Message(
            m.Heading(self.tr('Preclipping input data...')),
            m.Paragraph(self.tr(
                'Modifying %1 to avoid intersections with the aggregation '
                'layer'
            ).arg(theQgisLayer.name())))
        self._sendMessage(myMessage)

        theLayerFilename = str(theQgisLayer.source())
        myPostprocPolygons = self.safeLayer.get_geometry()
        myPolygonsLayer = safe_read_layer(theLayerFilename)
        myRemainingPolygons = numpy.array(myPolygonsLayer.get_geometry())
#        myRemainingAttributes = numpy.array(myPolygonsLayer.get_data())
        myRemainingIndexes = numpy.array(range(len(myRemainingPolygons)))

        #used for unit tests only
        self.preprocessedFeatureCount = 0

        # FIXME (MB) the intersecting array is used only for debugging and
        # could be safely removed
        myIntersectingPolygons = []
        myInsidePolygons = []

        # FIXME (MB) maybe do raw geos without qgis
        #select all postproc polygons with no attributes
        aggregationProvider = self.layer.dataProvider()
        aggregationProvider.select([])

        # copy polygons to a memory layer
        myQgisMemoryLayer = create_memory_layer(theQgisLayer)

        polygonsProvider = myQgisMemoryLayer.dataProvider()
        allPolygonAttrs = polygonsProvider.attributeIndexes()
        polygonsProvider.select(allPolygonAttrs)
        myQgisPostprocPoly = QgsFeature()
        myQgisFeat = QgsFeature()
        myInsideFeat = QgsFeature()
        fields = polygonsProvider.fields()
        myTempdir = temp_dir(sub_dir='preprocess')
        myOutFilename = unique_filename(suffix='.shp',
                                        dir=myTempdir)

        self.keywordIO.copy_keywords(theQgisLayer, myOutFilename)
        mySHPWriter = QgsVectorFileWriter(myOutFilename,
                                          'UTF-8',
                                          fields,
                                          polygonsProvider.geometryType(),
                                          polygonsProvider.crs())
        if mySHPWriter.hasError():
            raise InvalidParameterError(mySHPWriter.errorMessage())
        # end FIXME

        for (myPostprocPolygonIndex,
             myPostprocPolygon) in enumerate(myPostprocPolygons):
            LOGGER.debug('PostprocPolygon %s' % myPostprocPolygonIndex)
            myPolygonsCount = len(myRemainingPolygons)
            aggregationProvider.featureAtId(
                myPostprocPolygonIndex, myQgisPostprocPoly, True, [])
            myQgisPostprocGeom = QgsGeometry(myQgisPostprocPoly.geometry())

            # myPostprocPolygon bounding box values
            A = numpy.array(myPostprocPolygon)
            minx = miny = sys.maxint
            maxx = maxy = -minx
            myPostprocPolygonMinx = min(minx, min(A[:, 0]))
            myPostprocPolygonMaxx = max(maxx, max(A[:, 0]))
            myPostprocPolygonMiny = min(miny, min(A[:, 1]))
            myPostprocPolygonMaxy = max(maxy, max(A[:, 1]))

            # create an array full of False to store if a BB vertex is inside
            # or outside the myPostprocPolygon
            myAreVerticesInside = numpy.zeros(myPolygonsCount * 4,
                                              dtype=numpy.bool)

            # Create Nx2 vector of vertices of bounding boxes
            myBBVertices = []
            # Compute bounding box for each geometry type
            for myPoly in myRemainingPolygons:
#.........这里部分代码省略.........