QGIS » QGIS Application

# Analyzing raster data

registry = QgsProviderRegistry.instance()

provider = registry.createProvider("gdal", "C:/Users/kimpr/Desktop/QGIS_TrainingTutorial_Data/raster/SRTM/srtm_41_19.tif")

extent = provider.extent()

width = provider.xSize()

height = provider.ySize()

no_data_value = provider.sourceHasNoDataValue(1)

histogram = {} # mapping the elevation to number of occurrences

# Read the raster by using the block method

block = provider.block(1, extent, width, height)

# Read all the data (1) from the raster file

for x in range(width):

    for y in range(height):

        elevation = block.value(x, y)

        if elevation == no_data_value: continue

        try:

            histogram[elevation] += 1

        except KeyError:

            histogram[elevation] = 1

for height in sorted(histogram.keys()):

    print (height, histogram[height])