Python FToolsUtils.createSpatialIndex方法代码示例

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        polyLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        pointLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POINTS))
        fieldName = self.getParameterValue(self.FIELD)
        fieldIdx = pointLayer.fieldNameIndex(self.getParameterValue(self.WEIGHT))

        polyProvider = polyLayer.dataProvider()

        idxCount, fieldList = utils.findOrCreateField(polyLayer, polyLayer.pendingFields(), fieldName)

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

        spatialIndex = utils.createSpatialIndex(pointLayer)

        ftPoint = QgsFeature()
        outFeat = QgsFeature()
        geom = QgsGeometry()

        current = 0
        hasIntersections = False

        features = QGisLayers.features(polyLayer)
        total = 100.0 / float(len(features))
        for ftPoly in features:
            geom = ftPoly.geometry()
            attrs = ftPoly.attributes()

            count = 0
            hasIntersections = False
            points = spatialIndex.intersects(geom.boundingBox())
            if len(points) > 0:
                hasIntersections = True

            if hasIntersections:
                progress.setText(str(len(points)))
                for i in points:
                    request = QgsFeatureRequest().setFilterFid(i)
                    ftPoint = pointLayer.getFeatures(request).next()
                    tmpGeom = QgsGeometry(ftPoint.geometry())
                    if geom.contains(tmpGeom):
                        weight = str(ftPoint.attributes()[fieldIdx])
                        try:
                            count += float(weight)
                        except:
                            pass  # ignore fields with non-numeric values

            outFeat.setGeometry(geom)
            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 sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        filename = self.getParameterValue(self.INPUT)
        inputLayer = QGisLayers.getObjectFromUri(filename)
        method = self.getParameterValue(self.METHOD)
        filename = self.getParameterValue(self.INTERSECT)
        selectLayer = QGisLayers.getObjectFromUri(filename)

        oldSelection = set(inputLayer.selectedFeaturesIds())
        index = spatialIndex = utils.createSpatialIndex(inputLayer)

        feat = QgsFeature()
        geom = QgsGeometry()
        selectedSet = []
        current = 0
        features = QGisLayers.features(selectLayer)
        total = 100.0 / float(len(features))
        for f in features:
            geom = QgsGeometry(f.geometry())
            intersects = index.intersects(geom.boundingBox())
            for i in intersects:
                request = QgsFeatureRequest().setFilterFid(i)
                feat = inputLayer.getFeatures(request).next()
                tmpGeom = QgsGeometry(feat.geometry())
                if geom.intersects(tmpGeom):
                    selectedSet.append(feat.id())
            current += 1
            progress.setPercentage(int(current * total))

        if method == 1:
            selectedSet = list(oldSelection.union(selectedSet))
        elif method == 2:
            selectedSet = list(oldSelection.difference(selectedSet))

        inputLayer.setSelectedFeatures(selectedSet)
        self.setOutputValue(self.OUTPUT, filename)

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        polyLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        pointLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POINTS))
        fieldName = self.getParameterValue(self.FIELD)

        polyProvider = polyLayer.dataProvider()
        pointProvider = pointLayer.dataProvider()
        if polyProvider.crs() != pointProvider.crs():
            SextanteLog.addToLog(SextanteLog.LOG_WARNING,
                                 "CRS warning: Input layers have non-matching CRS. This may cause unexpected results.")

        idxCount, fieldList = utils.findOrCreateField(polyLayer, polyLayer.pendingFields(), fieldName)

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

        spatialIndex = utils.createSpatialIndex(pointLayer)

        pointProvider.rewind()
        pointProvider.select()

        allAttrs = polyLayer.pendingAllAttributesList()
        polyLayer.select(allAttrs)

        ftPoly = QgsFeature()
        ftPoint = QgsFeature()
        outFeat = QgsFeature()
        geom = QgsGeometry()

        current = 0
        hasIntersections = False

        features = QGisLayers.features(polyLayer)
        total = 100.0 / float(len(features))
        for ftPoly in features:
            geom = ftPoly.geometry()
            atMap = ftPoly.attributeMap()

            count = 0
            hasIntersections = False
            points = spatialIndex.intersects(geom.boundingBox())
            if len(points) > 0:
                hasIntersections = True

            if hasIntersections:
                for i in points:
                    pointLayer.featureAtId(int(i), ftPoint, True, False)
                    tmpGeom = QgsGeometry(ftPoint.geometry())
                    if geom.contains(tmpGeom):
                        count += 1

            outFeat.setGeometry(geom)
            outFeat.setAttributeMap(atMap)
            outFeat.addAttribute(idxCount, QVariant(count))
            writer.addFeature(outFeat)

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

        del writer

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        polyLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        pointLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POINTS))
        fieldName = self.getParameterValue(self.FIELD)
        classFieldName = self.getParameterValue(self.CLASSFIELD)

        polyProvider = polyLayer.dataProvider()
        pointProvider = pointLayer.dataProvider()
        if polyProvider.crs() != pointProvider.crs():
            SextanteLog.addToLog(SextanteLog.LOG_WARNING,
                                 "CRS warning: Input layers have non-matching CRS. This may cause unexpected results.")

        classFieldIndex = pointProvider.fieldNameIndex(classFieldName)
        idxCount, fieldList = utils.findOrCreateField(polyLayer, polyLayer.pendingFields(), fieldName)

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

        spatialIndex = utils.createSpatialIndex(pointLayer)

        ftPoint = QgsFeature()
        outFeat = QgsFeature()
        geom = QgsGeometry()

        current = 0
        hasIntersections = False

        features = QGisLayers.features(polyLayer)
        total = 100.0 / float(len(features))
        for ftPoly in features:
            geom = ftPoly.geometry()
            atMap = ftPoly.attributes()

            classes = []
            hasIntersections = False
            points = spatialIndex.intersects(geom.boundingBox())
            if len(points) > 0:
                hasIntersections = True

            if hasIntersections:
                for i in points:
                    pointLayer.featureAtId(int(i), ftPoint, True, True)
                    tmpGeom = QgsGeometry(ftPoint.geometry())
                    if geom.contains(tmpGeom):
                        clazz = ftPoint.attributes()[classFieldIndex].toString()
                        if not clazz in classes:
                            classes.append(clazz)

            outFeat.setGeometry(geom)        
            if idxCount == len(atMap):
                atMap.append(QVariant(len(classes)))
            else:
                atMap[idxCount] =  QVariant(len(classes))
            outFeat.setAttributes(atMap)
            writer.addFeature(outFeat)

            current += 1
            progress.setPercentage(current / total)

        del writer

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        layerA = QGisLayers.getObjectFromUri(self.getParameterValue(Difference.INPUT))
        layerB = QGisLayers.getObjectFromUri(self.getParameterValue(Difference.OVERLAY))

        GEOS_EXCEPT = True

        FEATURE_EXCEPT = True

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

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

        index = utils.createSpatialIndex(layerB)

        selectionA = QGisLayers.features(layerA)

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

        for inFeatA in selectionA:
            add = True
            geom = QgsGeometry(inFeatA.geometry())
            diff_geom = QgsGeometry(geom)
            attrs = inFeatA.attributes()
            intersections = index.intersects(geom.boundingBox())
            for i in intersections:
                request = QgsFeatureRequest().setFilterFid(i)
                inFeatB = layerB.getFeatures(request).next()
                tmpGeom = QgsGeometry(inFeatB.geometry())
                try:
                    if diff_geom.intersects(tmpGeom):
                        diff_geom = QgsGeometry(diff_geom.difference(tmpGeom))
                except:
                    GEOS_EXCEPT = False
                    add = False
                    break

            if add:
                try:
                    outFeat.setGeometry(diff_geom)
                    outFeat.setAttributes(attrs)
                    writer.addFeature(outFeat)
                except:
                    FEATURE_EXCEPT = False
                    continue

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

        del writer

        if not GEOS_EXCEPT:
            SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Geometry exception while computing difference")
        if not FEATURE_EXCEPT:
            SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Feature exception while computing difference")

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = QGisLayers.getObjectFromUri(self.getParameterValue(Difference.INPUT))
        vlayerB = QGisLayers.getObjectFromUri(self.getParameterValue(Difference.INPUT2))
        GEOS_EXCEPT = True
        FEATURE_EXCEPT = True
        vproviderA = vlayerA.dataProvider()
        vproviderB = vlayerB.dataProvider()
        fields = vproviderA.fields()
        # check for crs compatibility
        crsA = vproviderA.crs()
        crsB = vproviderB.crs()
        if not crsA.isValid() or not crsB.isValid():
            SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Difference. Invalid CRS. Results might be unexpected")
        else:
            if crsA != crsB:
                SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Difference. Non-matching CRSs. Results might be unexpected")
        writer = self.getOutputFromName(Difference.OUTPUT).getVectorWriter(fields, vproviderA.geometryType(), vproviderA.crs() )
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        index = utils.createSpatialIndex(vlayerB)
        nElement = 0
        selectionA = QGisLayers.features(vlayerA)
        nFeat = len(selectionA)
        for inFeatA in selectionA:
            nElement += 1
            progress.setPercentage(nElement/float(nFeat) * 100)
            add = True
            geom = QgsGeometry( inFeatA.geometry() )
            diff_geom = QgsGeometry( geom )
            atMap = inFeatA.attributes()
            intersects = index.intersects( geom.boundingBox() )
            for id in intersects:
                vlayerB.featureAtId( int( id ), inFeatB , True)
                tmpGeom = QgsGeometry( inFeatB.geometry() )
                try:
                    if diff_geom.intersects( tmpGeom ):
                        diff_geom = QgsGeometry( diff_geom.difference( tmpGeom ) )
                except:
                    GEOS_EXCEPT = False
                    add = False
                    break
            if add:
                try:
                    outFeat.setGeometry( diff_geom )
                    outFeat.setAttributes( atMap )
                    writer.addFeature( outFeat )
                except:
                    FEATURE_EXCEPT = False
                    continue


        del writer
        if not GEOS_EXCEPT:
            SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Geometry exception while computing difference")
        if not FEATURE_EXCEPT:
            SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Feature exception while computing difference")

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        polyLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POLYGONS))
        pointLayer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POINTS))
        fieldName = self.getParameterValue(self.FIELD)
        classFieldName = self.getParameterValue(self.CLASSFIELD)

        polyProvider = polyLayer.dataProvider()

        classFieldIndex = pointLayer.fieldNameIndex(classFieldName)
        idxCount, fieldList = utils.findOrCreateField(polyLayer, polyLayer.pendingFields(), fieldName)

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

        spatialIndex = utils.createSpatialIndex(pointLayer)

        ftPoint = QgsFeature()
        outFeat = QgsFeature()
        geom = QgsGeometry()

        current = 0
        hasIntersections = False

        features = QGisLayers.features(polyLayer)
        total = 100.0 / float(len(features))
        for ftPoly in features:
            geom = ftPoly.geometry()
            attrs = ftPoly.attributes()

            classes = []
            hasIntersections = False
            points = spatialIndex.intersects(geom.boundingBox())
            if len(points) > 0:
                hasIntersections = True

            if hasIntersections:
                for i in points:
                    request = QgsFeatureRequest().setFilterFid(i)
                    ftPoint = pointLayer.getFeatures(request).next()
                    tmpGeom = QgsGeometry(ftPoint.geometry())
                    if geom.contains(tmpGeom):
                        clazz = ftPoint.attributes()[classFieldIndex]
                        if not clazz in classes:
                            classes.append(clazz)

            outFeat.setGeometry(geom)
            if idxCount == len(attrs):
                attrs.append(len(classes))
            else:
                attrs[idxCount] =  len(classes)
            outFeat.setAttributes(attrs)
            writer.addFeature(outFeat)

            current += 1
            progress.setPercentage(current / total)

        del writer

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        filename = self.getParameterValue(self.INPUT)
        inputLayer = QGisLayers.getObjectFromUri(filename)
        method = self.getParameterValue(self.METHOD)
        selection = self.getParameterValue(self.USE_SELECTED)

        filename = self.getParameterValue(self.INTERSECT)
        selectLayer = QGisLayers.getObjectFromUri(filename)

        inputProvider = inputLayer.dataProvider()
        selectProvider = selectLayer.dataProvider()

        index = utils.createSpatialIndex(inputLayer)

        inputProvider.select()
        selectProvider.select()

        feat = QgsFeature()
        infeat = QgsFeature()
        geom = QgsGeometry()
        selectedSet = []

        if selection:
            features = selectLayer.selectedFeatures()
            total = 100.0 / float(len(features))
            current = 0
            for feat in features:
                geom = QgsGeometry(feat.geometry())
                intersects = index.intersects(geom.boundingBox())
                for i in intersects:
                    inputProvider.featureAtId(i, infeat, True)
                    tmpGeom = QgsGeometry(infeat.geometry())
                    if geom.intersects(tmpGeom):
                        selectedSet.append(infeat.id())
                current += 1
                progress.setPercentage(int(current * total))
        else:
            total = 100.0 / float(selectProvider.featureCount())
            current = 0
            while selectProvider.nextFeature(feat):
                geom = QgsGeometry(feat.geometry())
                intersects = index.intersects(geom.boundingBox())
                for i in intersects:
                    inputProvider.featureAtId(i, infeat, True)
                    tmpGeom = QgsGeometry(infeat.geometry())
                    if geom.intersects(tmpGeom):
                        selectedSet.append(infeat.id())
                current += 1
                progress.setPercentage(int(current * total))

        if method == 1:
            selectedSet = list(set(inputLayer.selectedFeaturesIds()).union(selectedSet))
        elif method == 2:
            selectedSet = list(set(inputLayer.selectedFeaturesIds()).difference(selectedSet))

        inputLayer.setSelectedFeatures(selectedSet)
        self.setOutputValue(self.OUTPUT, filename)

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        layer = QGisLayers.getObjectFromUri(self.getParameterValue(self.POINTS))
        output = self.getOutputValue(self.OUTPUT)

        provider = layer.dataProvider()
        spatialIndex = utils.createSpatialIndex(provider)
        provider.rewind()
        provider.select()

        feat = QgsFeature()
        neighbour = QgsFeature()
        distance = QgsDistanceArea()

        sumDist = 0.00
        A = layer.extent()
        A = float(A.width() * A.height())

        current = 0
        total = 100.0 / float(provider.featureCount())

        while provider.nextFeature( feat ):
            neighbourID = spatialIndex.nearestNeighbor(feat.geometry().asPoint(), 2)[1]
            provider.featureAtId(neighbourID, neighbour, True)
            sumDist += distance.measureLine(neighbour.geometry().asPoint(), feat.geometry().asPoint())

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

        count = provider.featureCount()
        do = float(sumDist) / count
        de = float(0.5 / math.sqrt(count / A))
        d = float(do / de)
        SE = float(0.26136 / math.sqrt(( count ** 2) / A))
        zscore = float((do - de) / SE)

        data = []
        data.append("Observed mean distance: " + unicode(do))
        data.append("Expected mean distance: " + unicode(de))
        data.append("Nearest neighbour index: " + unicode(d))
        data.append("Number of points: " + unicode(count))
        data.append("Z-Score: " + unicode(zscore))

        self.createHTML(output, data)

        self.setOutputValue(self.OBSERVED_MD, float( data[ 0 ].split( ": " )[ 1 ] ) )
        self.setOutputValue(self.EXPECTED_MD, float( data[ 1 ].split( ": " )[ 1 ] ) )
        self.setOutputValue(self.NN_INDEX, float( data[ 2 ].split( ": " )[ 1 ] ) )
        self.setOutputValue(self.POINT_COUNT, float( data[ 3 ].split( ": " )[ 1 ] ) )
        self.setOutputValue(self.Z_SCORE, float( data[ 4 ].split( ": " )[ 1 ] ) )

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def linearMatrix(self, inLayer, inField, targetLayer, targetField, matType, nPoints, progress):
        if matType == 0:
            self.writer.writerow(["InputID", "TargetID", "Distance"])
        else:
            self.writer.writerow(["InputID", "MEAN", "STDDEV", "MIN", "MAX"])

        index = utils.createSpatialIndex(targetLayer);

        inIdx = inLayer.fieldNameIndex(inField)
        inLayer.select([inIdx])
        outIdx = targetLayer.fieldNameIndex(inField)

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

        features = QGisLayers.features(inLayer)
        current = 0
        total = 100.0 / float(len(features))
        for inFeat in features:
            inGeom = inFeat.geometry()
            inID = inFeat.attributes()[inIdx].toString()
            featList = index.nearestNeighbor(inGeom.asPoint(), nPoints)
            distList = []
            vari = 0.0
            for i in featList:
                request = QgsFeatureRequest().setFilterFid(i)
                outFeat = targetLayer.getFeatures(request).next()
                outID = outFeat.attributes()[outIdx].toString()
                outGeom = outFeat.geometry()
                dist = distArea.measureLine(inGeom.asPoint(), outGeom.asPoint())
                if matType == 0:
                    self.writer.writerow([unicode(inID), unicode(outID), unicode(dist)])
                else:
                    distList.append(float(dist))

            if matType == 2:
                mean = sum(distList) / len(distList)
                for i in distList:
                    vari += (i - mean) * (i - mean)
                vari = math.sqrt(vari / len(distList))
                self.writer.writerow([unicode(inID), unicode(mean), unicode(vari), unicode(min(distList)), unicode(max(distList))])

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

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = QGisLayers.getObjectFromUri(self.getParameterValue(Intersection.INPUT))
        vlayerB = QGisLayers.getObjectFromUri(self.getParameterValue(Intersection.INPUT2))
        vproviderA = vlayerA.dataProvider()

        fields = utils.combineVectorFields(vlayerA, vlayerB)
        writer = self.getOutputFromName(Intersection.OUTPUT).getVectorWriter(fields, vproviderA.geometryType(), vproviderA.crs() )
        inFeatA = QgsFeature()
        inFeatB = QgsFeature()
        outFeat = QgsFeature()
        index = utils.createSpatialIndex(vlayerB)
        nElement = 0
        selectionA = QGisLayers.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 id in intersects:
                vlayerB.featureAtId( int( id ), inFeatB , True)
                tmpGeom = QgsGeometry( inFeatB.geometry() )
                try:
                    if geom.intersects( tmpGeom ):
                        atMapB = inFeatB.attributes()
                        int_geom = QgsGeometry( geom.intersection( tmpGeom ) )
                        if int_geom.wkbType() == 7:
                            int_com = geom.combine( tmpGeom )
                            int_sym = geom.symDifference( tmpGeom )
                            int_geom = QgsGeometry( int_com.difference( int_sym ) )
                    try:
                        outFeat.setGeometry( int_geom )
                        attrs = []
                        attrs.extend(atMapA)
                        attrs.extend(atMapB)
                        outFeat.setAttributes(attrs)
                        writer.addFeature( outFeat )
                    except:
                        raise GeoAlgorithmExecutionException("Feature exception while computing intersection")
                except:
                    raise GeoAlgorithmExecutionException("Geometry exception while computing intersection")


        del writer

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def regularMatrix(self, inLayer, inField, targetLayer, targetField, nPoints, progress):
        index = utils.createSpatialIndex(targetLayer)

        inIdx = inLayer.fieldNameIndex(inField)
        outIdx = targetLayer.fieldNameIndex(inField)

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

        first = True
        current = 0
        features = QGisLayers.features(inLayer)
        total = 100.0 / float(len(features))

        for inFeat in features:
            inGeom = inFeat.geometry()
            inID = inFeat.attributes()[inIdx].toString()
            if first:
                featList = index.nearestNeighbor(inGeom.asPoint(), nPoints)
                first = False
                data = ["ID"]
                for i in featList:
                    request = QgsFeatureRequest().setFilterFid(i)
                    outFeat = targetLayer.getFeatures(request).next()
                    data.append(unicode(outFeat.attributes[outIdx].toString()))
                self.writer.writerow(data)

            data = [unicode(inID)]
            for i in featList:
                request = QgsFeatureRequest().setFilterFid(i)
                outFeat = targetLayer.getFeatures(request).next()
                outGeom = outFeat.geometry()
                dist = distArea.measureLine(inGeom.asPoint(), outGeom.asPoint())
                data.append(unicode(float(dist)))
            self.writer.writerow(data)

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

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        layerA = QGisLayers.getObjectFromUri(self.getParameterValue(self.INPUT_A))
        layerB = QGisLayers.getObjectFromUri(self.getParameterValue(self.INPUT_B))
        fieldA = self.getParameterValue(self.FIELD_A)
        fieldB = self.getParameterValue(self.FIELD_B)

        providerA = layerA.dataProvider()
        providerB = layerB.dataProvider()


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

        fieldList = {0 : layerA.pendingFields()[idxA],
                     1 : layerB.pendingFields()[idxB]
                    }

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

        spatialIndex = utils.createSpatialIndex(layerB)

        providerA.rewind()
        layerA.select([idxA])
        providerB.rewind()
        layerB.select([idxB])

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

        current = 0
        total = 100.0 / float(providerA.featureCount())
        hasIntersections = False

        while layerA.nextFeature(inFeatA):
            inGeom = inFeatA.geometry()
            hasIntersections = False
            lines = spatialIndex.intersects(inGeom.boundingBox())

            if len(lines) > 0:
                hasIntersections = True

            if hasIntersections:
                layerB.select([idxB])
                for i in lines:
                    layerB.featureAtId(int(i), inFeatB)
                    tmpGeom = QgsGeometry(inFeatB.geometry())

                    points = []
                    atMapA = inFeatA.attributeMap()
                    atMapB = inFeatB.attributeMap()

                    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.addAttribute(0, atMapA[idxA])
                                outFeat.addAttribute(1, atMapB[idxB])
                                writer.addFeature(outFeat)

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

        del writer

# 需要导入模块: from sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        layerA = QGisLayers.getObjectFromUri(self.getParameterValue(self.INPUT_A))
        layerB = QGisLayers.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 = utils.createSpatialIndex(layerB)

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

        features = QGisLayers.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 sextante.algs.ftools import FToolsUtils [as 别名]
# 或者: from sextante.algs.ftools.FToolsUtils import createSpatialIndex [as 别名]
    def processAlgorithm(self, progress):
        vlayerA = QGisLayers.getObjectFromUri(self.getParameterValue(Union.INPUT))
        vlayerB = QGisLayers.getObjectFromUri(self.getParameterValue(Union.INPUT2))
        GEOS_EXCEPT = True
        FEATURE_EXCEPT = True
        vproviderA = vlayerA.dataProvider()

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

        count = 0
        nElement = 0
        featuresA = QGisLayers.features(vlayerA)
        nFeat = len(featuresA)
        for inFeatA in featuresA:
          progress.setPercentage(nElement/float(nFeat) * 50)
          nElement += 1
          found = False
          geom = QgsGeometry( inFeatA.geometry() )
          diff_geom = QgsGeometry( geom )
          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("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 ):
                  found = True
                  int_geom = geom.intersection( tmpGeom )

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

                  if diff_geom.intersects( tmpGeom ):
                    diff_geom = diff_geom.difference( tmpGeom )
                    if diff_geom is None:
                      # It's possible there was an error here?
                      diff_geom = QgsGeometry()
                    else:
                      diff_geom = QgsGeometry(diff_geom)

                  if int_geom.wkbType() == 0:
                    # 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, err:
                    raise GeoAlgorithmExecutionException("Feature exception while computing union")
                else:
                  # this only happends if the bounding box
                  # intersects, but the geometry doesn't
                  try:
                    outFeat.setGeometry( geom )
                    outFeat.setAttributes( atMapA )
                    writer.addFeature( outFeat )
                  except:
                    # also shoudn't ever happen
                    raise GeoAlgorithmExecutionException("Feature exception while computing union")


            if found:
              try:
                if diff_geom.wkbType() == 0:
                  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 )
#.........这里部分代码省略.........