4

It seems that gdal can caused in some crashes in QGIS, specific in used for raster layers (such as links : QGIS crashes on close if data was read through GDAL in Python Console or GDAL/QGIS broken Ubuntu 14.04 or Why close QGIS for extract pixel values by GDAL?. So Is it possible to get pixel values only with pyqgis and without gdal?

Improve this question
2

3 Answers 3

Reset to default
5

Yes, it is possible to get pixel values without gdal; only PyQGIS. Easier way is with a QgsRasterBlock object. This is the code:

layer = iface.activeLayer() #in my case, a 20x20 raster

provider = layer.dataProvider()

extent = provider.extent()

rows = layer.height()
cols = layer.width()

block = provider.block(1, extent, cols, rows)

for i in range(rows):
    for j in range(cols):
        print block.value(i,j),
    print

After running the code at Python Console of QGIS, I got:

27.0 17.0 48.0 60.0 32.0 53.0 73.0 5.0 94.0 35.0 18.0 31.0 35.0 63.0 54.0 26.0 79.0 17.0 3.0 98.0
32.0 44.0 7.0 9.0 54.0 52.0 12.0 19.0 75.0 26.0 27.0 28.0 36.0 58.0 25.0 83.0 83.0 67.0 70.0 6.0
27.0 17.0 44.0 63.0 82.0 76.0 64.0 69.0 42.0 52.0 81.0 4.0 46.0 37.0 97.0 86.0 88.0 44.0 17.0 16.0
36.0 70.0 53.0 91.0 34.0 9.0 90.0 47.0 39.0 23.0 30.0 27.0 36.0 89.0 99.0 48.0 23.0 55.0 56.0 77.0
75.0 96.0 75.0 64.0 36.0 72.0 55.0 87.0 35.0 56.0 23.0 6.0 97.0 58.0 97.0 54.0 14.0 93.0 9.0 54.0
41.0 58.0 52.0 3.0 96.0 45.0 32.0 85.0 32.0 98.0 9.0 26.0 63.0 20.0 83.0 10.0 67.0 24.0 58.0 19.0
61.0 5.0 15.0 17.0 53.0 57.0 57.0 59.0 79.0 43.0 7.0 65.0 83.0 65.0 32.0 59.0 22.0 24.0 25.0 22.0
66.0 60.0 3.0 94.0 53.0 93.0 10.0 8.0 24.0 16.0 92.0 95.0 79.0 42.0 97.0 45.0 82.0 87.0 31.0 8.0
2.0 40.0 34.0 5.0 6.0 37.0 9.0 83.0 28.0 11.0 54.0 84.0 35.0 30.0 46.0 24.0 49.0 15.0 30.0 36.0
30.0 45.0 65.0 64.0 18.0 49.0 78.0 70.0 83.0 64.0 48.0 73.0 69.0 2.0 69.0 1.0 30.0 20.0 45.0 74.0
7.0 93.0 76.0 96.0 60.0 84.0 7.0 81.0 79.0 6.0 47.0 80.0 9.0 61.0 43.0 19.0 53.0 100.0 28.0 49.0
77.0 4.0 79.0 84.0 65.0 14.0 91.0 62.0 41.0 1.0 78.0 26.0 100.0 66.0 37.0 8.0 96.0 51.0 83.0 94.0
100.0 28.0 77.0 7.0 3.0 72.0 82.0 97.0 3.0 25.0 64.0 66.0 25.0 75.0 54.0 29.0 15.0 36.0 25.0 57.0
7.0 63.0 64.0 46.0 41.0 13.0 6.0 34.0 55.0 32.0 99.0 87.0 12.0 100.0 9.0 31.0 48.0 22.0 7.0 52.0
45.0 43.0 9.0 17.0 39.0 98.0 20.0 51.0 80.0 91.0 52.0 18.0 40.0 39.0 5.0 97.0 74.0 5.0 9.0 99.0
67.0 53.0 24.0 67.0 32.0 17.0 35.0 98.0 37.0 68.0 18.0 31.0 71.0 95.0 59.0 100.0 97.0 60.0 22.0 6.0
97.0 73.0 69.0 93.0 13.0 87.0 29.0 69.0 10.0 62.0 35.0 26.0 62.0 51.0 30.0 15.0 25.0 62.0 1.0 61.0
59.0 53.0 22.0 74.0 5.0 52.0 13.0 29.0 14.0 18.0 9.0 99.0 70.0 3.0 53.0 100.0 39.0 20.0 73.0 30.0
66.0 72.0 92.0 26.0 77.0 67.0 74.0 60.0 55.0 68.0 71.0 61.0 32.0 72.0 100.0 49.0 80.0 74.0 8.0 1.0
53.0 4.0 82.0 72.0 86.0 82.0 92.0 22.0 52.0 8.0 31.0 30.0 6.0 2.0 14.0 73.0 34.0 60.0 84.0 96.0

I corroborated, with help of Value Tool plugin, that values extracted were correct.

Improve this answer
1
  • Thanks, It's simple but effective code.
    – HMadadi
    Jul 13, 2017 at 6:40
3

Based on @Michael Stimson in above and QGIS document (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/raster.html#query-values). In below example it is possible to get pixel values without gdal.

from qgis.core import *

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return QgsPoint(xp, yp)

# Open tif file
ds = QgsRasterLayer("/Data/Sample_Polygon.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

for x, y in pntRstList:
    ident = ds.dataProvider().identify(QgsPoint(x, y), 
    QgsRaster.IdentifyFormatValue)
    if ident.isValid():
        print ident.results()
Improve this answer
1
  • 1
    Yes, it is possible to get pixel values without gdal, only PyQGIS, but easier way is with a QgsRasterBlock object.
    – xunilk
    Jul 13, 2017 at 0:47
1

the first solution is a really good one! I just would like to propose instead of print or instead of put it in a list like in the second solution, put the values in numpy array. If you are planning to work with large maps is better to take this approach once the numpy array is much faster.

    from qgis.core import *
    import numpy as np

    def GetValuePixels(layer):
        provider = layer.dataProvider()

        extent = provider.extent()

        rows = layer.height()
        cols = layer.width()

        block = provider.block(1, extent, cols, rows)
        values = np.zeros(shape=(rows,cols)) #Creates space in the memory to avoid spend time...

        for i in range(rows):
            for j in range(cols):
                values[i,j] = block.value(i,j)

        return values

thanks!

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.