Python analysis.QgsZonalStatistics类代码示例

 def processAlgorithm(self, parameters, context, feedback):
     zs = QgsZonalStatistics(self.vectorLayer,
                             self.rasterLayer,
                             self.columnPrefix,
                             self.bandNumber,
                             QgsZonalStatistics.Statistics(self.selectedStats))
     zs.calculateStatistics(feedback)
     return {self.INPUT_VECTOR: self.vectorLayer}

    def processAlgorithm(self, feedback):
        rasterPath = self.getParameterValue(self.INPUT_RASTER)
        vectorPath = self.getParameterValue(self.INPUT_VECTOR)
        bandNumber = self.getParameterValue(self.RASTER_BAND)
        columnPrefix = self.getParameterValue(self.COLUMN_PREFIX)
        st = self.getParameterValue(self.STATISTICS)

        vectorLayer = dataobjects.getObjectFromUri(vectorPath)

        keys = list(self.STATS.keys())
        selectedStats = 0
        for i in st:
            selectedStats |= self.STATS[keys[i]]

        zs = QgsZonalStatistics(vectorLayer, rasterPath, columnPrefix, bandNumber, selectedStats)
        zs.calculateStatistics(None)

        self.setOutputValue(self.OUTPUT_LAYER, vectorPath)

    def testStatistics(self):
        """Test zonal stats"""
        sep = os.sep
        TEST_DATA_DIR = unitTestDataPath() + sep + "zonalstatistics" + sep
        myTempPath = QDir.tempPath() + sep
        testDir = QDir(TEST_DATA_DIR)
        for f in testDir.entryList(QDir.Files):
            QFile.remove(myTempPath + f)
            QFile.copy(TEST_DATA_DIR + f, myTempPath + f)

        myVector = QgsVectorLayer(myTempPath + "polys.shp", "poly", "ogr")
        myRasterPath = myTempPath + "edge_problem.asc"
        zs = QgsZonalStatistics(myVector, myRasterPath, "", 1)
        zs.calculateStatistics(None)

        feat = QgsFeature()
        # validate statistics for each feature
        request = QgsFeatureRequest().setFilterFid(0)
        feat = myVector.getFeatures(request).next()
        myMessage = ('Expected: %f\nGot: %f\n' % (12.0, feat[1]))
        assert feat[1] == 12.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (8.0, feat[2]))
        assert feat[2] == 8.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.666666666666667, feat[3]))
        assert abs(feat[3] - 0.666666666666667) < 0.00001, myMessage

        request.setFilterFid(1)
        feat = myVector.getFeatures(request).next()
        myMessage = ('Expected: %f\nGot: %f\n' % (9.0, feat[1]))
        assert feat[1] == 9.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (5.0, feat[2]))
        assert feat[2] == 5.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.555555555555556, feat[3]))
        assert abs(feat[3] - 0.555555555555556) < 0.00001, myMessage

        request.setFilterFid(2)
        feat = myVector.getFeatures(request).next()
        myMessage = ('Expected: %f\nGot: %f\n' % (6.0, feat[1]))
        assert feat[1] == 6.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (5.0, feat[2]))
        assert feat[2] == 5.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.833333333333333, feat[3]))
        assert abs(feat[3] - 0.833333333333333) < 0.00001, myMessage

    def processAlgorithm(self, parameters, context, feedback):
        bandNumber = self.parameterAsInt(parameters, self.RASTER_BAND, context)
        columnPrefix = self.parameterAsString(parameters, self.COLUMN_PREFIX, context)
        st = self.parameterAsEnums(parameters, self.STATISTICS, context)

        vectorLayer = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR, context)
        rasterLayer = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER, context)

        keys = list(self.STATS.keys())
        selectedStats = 0
        for i in st:
            selectedStats |= self.STATS[keys[i]]

        zs = QgsZonalStatistics(vectorLayer,
                                rasterLayer,
                                columnPrefix,
                                bandNumber,
                                selectedStats)
        zs.calculateStatistics(feedback)

        return {self.INPUT_VECTOR: vectorLayer}

    def testStatistics(self):
        """Test zonal stats"""
        TEST_DATA_DIR = unitTestDataPath() + "/zonalstatistics/"
        myTempPath = QDir.tempPath() + "/"
        testDir = QDir(TEST_DATA_DIR)
        for f in testDir.entryList(QDir.Files):
            QFile.remove(myTempPath + f)
            QFile.copy(TEST_DATA_DIR + f, myTempPath + f)

        myVector = QgsVectorLayer(myTempPath + "polys.shp", "poly", "ogr")
        myRaster = QgsRasterLayer(myTempPath + "edge_problem.asc", "raster", "gdal")
        zs = QgsZonalStatistics(myVector, myRaster, "", 1, QgsZonalStatistics.All)
        zs.calculateStatistics(None)

        feat = QgsFeature()
        # validate statistics for each feature
        request = QgsFeatureRequest().setFilterFid(0)
        feat = next(myVector.getFeatures(request))
        myMessage = ('Expected: %f\nGot: %f\n' % (12.0, feat[1]))
        assert feat[1] == 12.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (8.0, feat[2]))
        assert feat[2] == 8.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.666666666666667, feat[3]))
        assert abs(feat[3] - 0.666666666666667) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[4]))
        assert feat[4] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.47140452079103201, feat[5]))
        assert abs(feat[5] - 0.47140452079103201) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[6]))
        assert feat[6] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[7]))
        assert feat[7] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[8]))
        assert feat[8] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[9]))
        assert feat[9] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[10]))
        assert feat[10] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (2.0, feat[11]))
        assert feat[11] == 2.0, myMessage

        request.setFilterFid(1)
        feat = next(myVector.getFeatures(request))
        myMessage = ('Expected: %f\nGot: %f\n' % (9.0, feat[1]))
        assert feat[1] == 9.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (5.0, feat[2]))
        assert feat[2] == 5.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.555555555555556, feat[3]))
        assert abs(feat[3] - 0.555555555555556) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[4]))
        assert feat[4] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.49690399499995302, feat[5]))
        assert abs(feat[5] - 0.49690399499995302) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[6]))
        assert feat[6] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[7]))
        assert feat[7] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[8]))
        assert feat[8] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[9]))
        assert feat[9] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[10]))
        assert feat[10] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (2.0, feat[11]))
        assert feat[11] == 2.0, myMessage

        request.setFilterFid(2)
        feat = next(myVector.getFeatures(request))
        myMessage = ('Expected: %f\nGot: %f\n' % (6.0, feat[1]))
        assert feat[1] == 6.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (5.0, feat[2]))
        assert feat[2] == 5.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.833333333333333, feat[3]))
        assert abs(feat[3] - 0.833333333333333) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[4]))
        assert feat[4] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.372677996249965, feat[5]))
        assert abs(feat[5] - 0.372677996249965) < 0.00001, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[6]))
        assert feat[6] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[7]))
        assert feat[7] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[8]))
        assert feat[8] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (0.0, feat[9]))
        assert feat[9] == 0.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (1.0, feat[10]))
        assert feat[10] == 1.0, myMessage
        myMessage = ('Expected: %f\nGot: %f\n' % (2.0, feat[11]))
        assert feat[11] == 2.0, myMessage

def zonal_stats(raster, vector, callback=None):
    """Reclassify a continuous raster layer.

    Issue https://github.com/inasafe/inasafe/issues/3190


    :param raster: The raster layer.
    :type raster: QgsRasterLayer

    :param vector: The vector layer.
    :type vector: QgsVectorLayer

    :param callback: A function to all to indicate progress. The function
        should accept params 'current' (int), 'maximum' (int) and 'step' (str).
        Defaults to None.
    :type callback: function

    :return: The output of the zonal stats.
    :rtype: QgsVectorLayer

    .. versionadded:: 4.0
    """
    output_layer_name = zonal_stats_steps['output_layer_name']
    processing_step = zonal_stats_steps['step_name']

    layer = create_memory_layer(
        output_layer_name,
        vector.geometryType(),
        vector.crs(),
        vector.fields()
    )

    copy_layer(vector, layer)

    analysis = QgsZonalStatistics(
        layer, raster.source(), 'exposure_', 1, QgsZonalStatistics.Sum)
    result = analysis.calculateStatistics(None)
    LOGGER.debug(tr('Zonal stats on %s : %s' % (raster.source(), result)))

    layer.startEditing()
    exposure = raster.keywords['exposure']
    output_field = exposure_count_field['field_name'] % exposure
    fields_to_rename = {
        'exposure_sum': output_field
    }
    copy_fields(layer, fields_to_rename)
    remove_fields(layer, fields_to_rename.keys())

    layer.commitChanges()

    # The zonal stats is producing some None values. We need to fill these
    # with 0. See issue : #3778
    # We should start a new editing session as previous fields need to be
    # commited first.
    layer.startEditing()
    request = QgsFeatureRequest()
    expression = '\"%s\" is None' % output_field
    request.setFilterExpression(expression)
    request.setFlags(QgsFeatureRequest.NoGeometry)
    index = layer.fieldNameIndex(output_field)
    for feature in layer.getFeatures():
        if feature[output_field] is None:
            layer.changeAttributeValue(feature.id(), index, 0)
    layer.commitChanges()

    layer.keywords = raster.keywords.copy()
    layer.keywords['inasafe_fields'] = vector.keywords['inasafe_fields'].copy()
    layer.keywords['inasafe_default_values'] = (
        raster.keywords['inasafe_default_values'].copy())

    key = exposure_count_field['key'] % raster.keywords['exposure']

    # Special case here, one field is the exposure count and the total.
    layer.keywords['inasafe_fields'][key] = output_field
    layer.keywords['inasafe_fields'][total_field['key']] = output_field

    layer.keywords['exposure_keywords'] = raster.keywords.copy()
    layer.keywords['hazard_keywords'] = vector.keywords[
        'hazard_keywords'].copy()
    layer.keywords['aggregation_keywords'] = (
        vector.keywords['aggregation_keywords'])
    layer.keywords['layer_purpose'] = (
        layer_purpose_aggregate_hazard_impacted['key'])

    layer.keywords['title'] = output_layer_name

    check_layer(layer)
    return layer

def zonal_stats(raster, vector):
    """Reclassify a continuous raster layer.

    Issue https://github.com/inasafe/inasafe/issues/3190

    The algorithm will take care about projections.
    We don't want to reproject the raster layer.
    So if CRS are different, we reproject the vector layer and then we do a
    lookup from the reprojected layer to the original vector layer.

    :param raster: The raster layer.
    :type raster: QgsRasterLayer

    :param vector: The vector layer.
    :type vector: QgsVectorLayer

    :return: The output of the zonal stats.
    :rtype: QgsVectorLayer

    .. versionadded:: 4.0
    """
    output_layer_name = zonal_stats_steps['output_layer_name']

    exposure = raster.keywords['exposure']
    if raster.crs().authid() != vector.crs().authid():
        layer = reproject(vector, raster.crs())

        # We prepare the copy
        output_layer = create_memory_layer(
            output_layer_name,
            vector.geometryType(),
            vector.crs(),
            vector.fields()
        )
        copy_layer(vector, output_layer)
    else:
        layer = create_memory_layer(
            output_layer_name,
            vector.geometryType(),
            vector.crs(),
            vector.fields()
        )
        copy_layer(vector, layer)

    input_band = layer.keywords.get('active_band', 1)
    analysis = QgsZonalStatistics(
        layer,
        raster,
        'exposure_',
        input_band,
        QgsZonalStatistics.Sum)
    result = analysis.calculateStatistics(None)
    LOGGER.debug(tr('Zonal stats on %s : %s' % (raster.source(), result)))

    output_field = exposure_count_field['field_name'] % exposure
    if raster.crs().authid() != vector.crs().authid():
        output_layer.startEditing()
        field = create_field_from_definition(
            exposure_count_field, exposure)

        output_layer.addAttribute(field)
        new_index = output_layer.fields().lookupField(field.name())
        old_index = layer.fields().lookupField('exposure_sum')
        for feature_input, feature_output in zip(
                layer.getFeatures(), output_layer.getFeatures()):
            output_layer.changeAttributeValue(
                feature_input.id(), new_index, feature_input[old_index])
        output_layer.commitChanges()
        layer = output_layer
    else:
        fields_to_rename = {
            'exposure_sum': output_field
        }
        if qgis_version() >= 21600:
            rename_fields(layer, fields_to_rename)
        else:
            copy_fields(layer, fields_to_rename)
            remove_fields(layer, list(fields_to_rename.keys()))
        layer.commitChanges()

    # The zonal stats is producing some None values. We need to fill these
    # with 0. See issue : #3778
    # We should start a new editing session as previous fields need to be
    # committed first.
    layer.startEditing()
    request = QgsFeatureRequest()
    expression = '\"%s\" is None' % output_field
    request.setFilterExpression(expression)
    request.setFlags(QgsFeatureRequest.NoGeometry)
    index = layer.fields().lookupField(output_field)
    for feature in layer.getFeatures():
        if feature[output_field] is None:
            layer.changeAttributeValue(feature.id(), index, 0)
    layer.commitChanges()

    layer.keywords = raster.keywords.copy()
    layer.keywords['inasafe_fields'] = vector.keywords['inasafe_fields'].copy()
    layer.keywords['inasafe_default_values'] = (
        raster.keywords['inasafe_default_values'].copy())

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

####################################################

import os, sys
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry,  QgsZonalStatistics

# Add polygon layer 

#specify polygon shapefile vector
polygonLayer = QgsVectorLayer('/Users/amirsiraj/GoFlex/DataNow/Admin_units/south_America.shp', 'zonepolygons', "ogr") 


vtypes = ["cumI_SA_repl_", "cumPI_SA_repl_", "cumBI_SA_repl_"]
sufs = ["I","PI","BI"]

rplst = range(1,8,1) + range(9,101,1) 
rplst = range(1,8,1)+ range(9,59,1) + range(60,92,1) + range(93,101,1)
rplst = range (501,1001,1)
rng2 = range (1,2,1) + range(3,4,1)

for vtype in rng2 :
	for repl in rplst :  
		# Add raster 1
		rasterFilePath = '/Users/amirsiraj/GoFlex/zika/outputbk/'+ vtypes[vtype-1] + str(repl) +'_stat.bil'
		zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, sufs[vtype-1]+str(repl), 1, QgsZonalStatistics.Sum)
		zoneStat.calculateStatistics(None)
		print(vtype)
		print(repl)
		

    def _aggregateRasterImpact(self, theQGISImpactLayer):
        """
        Performs Aggregation postprocessing step on raster impact layers by
        calling QgsZonalStatistics
        Args:
            QgsMapLayer: theQGISImpactLayer a valid QgsVectorLayer

        Returns: None
        """
        myZonalStatistics = QgsZonalStatistics(
            self.layer,
            theQGISImpactLayer.dataProvider().dataSourceUri(),
            self.prefix)
        myProgressDialog = QtGui.QProgressDialog(
            self.tr('Calculating zonal statistics'),
            self.tr('Abort...'),
            0,
            0)
        startTime = time.clock()
        myZonalStatistics.calculateStatistics(myProgressDialog)
        if myProgressDialog.wasCanceled():
            QtGui.QMessageBox.error(
                self, self.tr('ZonalStats: Error'),
                self.tr('You aborted aggregation, '
                        'so there are no data for analysis. Exiting...'))
        cppDuration = time.clock() - startTime
        print 'CPP duration: %ss' % (cppDuration)

        startTime = time.clock()
        # new way
        # myZonalStatistics = {
        # 0L: {'count': 50539,
        #      'sum': 12015061.876953125,
        #      'mean': 237.73841739949594},
        # 1L: {
        #   'count': 19492,
        #   'sum': 2945658.1220703125,
        #   'mean': 151.12138939412642},
        # 2L: {
        #   'count': 57372,
        #   'sum': 1643522.3984985352, 'mean': 28.6467684323108},
        # 3L: {
        #   'count': 0.00013265823369700314,
        #   'sum': 0.24983273179242008,
        #   'mean': 1883.2810058593748},
        # 4L: {
        #   'count': 1.8158245316933218e-05,
        #   'sum': 0.034197078505115275,
        #   'mean': 1883.281005859375},
        # 5L: {
        #   'count': 73941,
        #   'sum': 10945062.435424805,
        #   'mean': 148.024268476553},
        # 6L: {
        #   'count': 54998,
        #   'sum': 11330910.488220215,
        #   'mean': 206.02404611477172}}

        myZonalStatistics = calculateZonalStats(theQGISImpactLayer, self.layer)
        pyDuration = time.clock() - startTime
        print 'CPP duration: %ss' % (pyDuration)

        try:
            ratio = pyDuration / cppDuration
        except ZeroDivisionError:
            ratio = 1

        print 'py to CPP: %s%%' % (ratio * 100)
        # FIXME (MB) remove this once fully implemented
        oldPrefix = self.prefix

        self.prefix = 'newAggr'
        myProvider = self.layer.dataProvider()
        self.layer.startEditing()

        # add fields for stats to aggregation layer
        # { 1: {'sum': 10, 'count': 20, 'min': 1, 'max': 4, 'mean': 2},
        #             QgsField(self._minFieldName(), QtCore.QVariant.Double),
        #             QgsField(self._maxFieldName(), QtCore.QVariant.Double)]
        myFields = [QgsField(self._countFieldName(), QtCore.QVariant.Double),
                    QgsField(self._sumFieldName(), QtCore.QVariant.Double),
                    QgsField(self._meanFieldName(), QtCore.QVariant.Double)
                    ]
        myProvider.addAttributes(myFields)
        self.layer.commitChanges()

        sumIndex = myProvider.fieldNameIndex(self._sumFieldName())
        countIndex = myProvider.fieldNameIndex(self._countFieldName())
        meanIndex = myProvider.fieldNameIndex(self._meanFieldName())
        # minIndex = myProvider.fieldNameIndex(self._minFieldName())
        # maxIndex = myProvider.fieldNameIndex(self._maxFieldName())

        self.layer.startEditing()
        allPolygonAttrs = myProvider.attributeIndexes()
        myProvider.select(allPolygonAttrs)
        myFeature = QgsFeature()

        while myProvider.nextFeature(myFeature):
            myFid = myFeature.id()
            myStats = myZonalStatistics[myFid]
#.........这里部分代码省略.........