1

In below code I want to print pixel values, but when run in QGIS python console QGIS software crashes and closed! what is the problem?

from osgeo import gdal
import numpy
import struct

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()

# extent of the Raster layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/Data/Sample_Polygon.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

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:
    structval=outBand.ReadRaster(x,y,1,1,buf_type=gdal.GDT_Float32)
    intval = struct.unpack('f' , structval)
    val=intval[0]
    print val

Also I test the @xunilk code in this link, and closed the QGIS software again. I use QGIS 2.18.9 in Linux Mint 18.

  • 'ReadRaster' GDAL method has a very different syntax that you used for unpacking data with struct. That's the issue. – xunilk Jul 12 '17 at 22:09
1

Issue is because 'ReadRaster' GDAL method has a very different syntax that you used for unpacking data with struct. Furthermore, you need 'BandType' to determine adequate 'fmttype' to be used in unpacking (see dictionary in last lines in next code). I change your paths to adequate to my example raster (20 x 20).

from osgeo import gdal
import numpy
import struct

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("/home/zeito/pyqgis_data/aleatorio1.tif","Study")

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

# extent of the Raster layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/home/zeito/pyqgis_data/aleatorio1.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

pntRstList = []

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

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
BandType = gdal.GetDataTypeName(outBand.DataType)

for y in range(outBand.YSize):
    scanline = outBand.ReadRaster(0, y, outBand.XSize, 1, outBand.YSize, 1, outBand.DataType)
    values = struct.unpack(fmttypes[BandType] * outBand.YSize, scanline)
    print values

After running it at Python Console of QGIS, I didn't have any error. This was result (20 tuples of 20 values):

(1, 6, 4, 9, 4, 7, 3, 2, 7, 9, 4, 8, 5, 5, 3, 6, 4, 7, 9, 6)
(10, 8, 1, 2, 4, 1, 1, 8, 8, 1, 4, 5, 7, 10, 8, 6, 3, 8, 8, 1)
(4, 5, 6, 9, 5, 4, 1, 10, 1, 4, 9, 8, 10, 6, 1, 9, 5, 7, 9, 7)
(3, 8, 5, 3, 2, 3, 10, 3, 7, 2, 3, 8, 3, 7, 7, 7, 10, 6, 5, 3)
(7, 9, 7, 7, 6, 3, 7, 7, 5, 10, 3, 1, 2, 2, 1, 3, 10, 2, 5, 4)
(7, 9, 1, 6, 3, 4, 3, 9, 8, 4, 5, 9, 8, 1, 4, 1, 5, 6, 6, 6)
(5, 4, 1, 8, 2, 4, 7, 4, 10, 9, 3, 8, 8, 3, 2, 8, 1, 2, 4, 6)
(7, 8, 8, 4, 9, 7, 5, 3, 4, 5, 7, 2, 8, 6, 5, 2, 4, 10, 6, 5)
(3, 9, 6, 9, 10, 1, 9, 5, 1, 3, 7, 4, 9, 5, 2, 4, 5, 4, 4, 9)
(1, 8, 9, 6, 2, 9, 8, 5, 10, 5, 7, 3, 5, 4, 6, 7, 1, 2, 8, 4)
(3, 6, 10, 3, 2, 9, 10, 6, 7, 7, 8, 3, 7, 7, 1, 5, 8, 1, 4, 9)
(10, 8, 6, 8, 6, 4, 8, 4, 3, 2, 10, 2, 10, 3, 1, 4, 9, 1, 2, 8)
(7, 4, 9, 2, 2, 10, 9, 10, 7, 4, 3, 3, 10, 4, 6, 6, 7, 3, 1, 7)
(4, 8, 9, 7, 5, 4, 10, 10, 1, 6, 9, 2, 2, 7, 9, 6, 4, 3, 8, 5)
(3, 5, 9, 4, 2, 4, 8, 8, 7, 1, 2, 4, 1, 4, 10, 7, 6, 3, 7, 3)
(2, 4, 8, 1, 10, 5, 2, 10, 7, 9, 7, 8, 3, 5, 5, 4, 4, 10, 2, 9)
(7, 5, 2, 9, 1, 1, 2, 8, 10, 10, 6, 7, 8, 6, 5, 7, 4, 3, 10, 9)
(7, 4, 10, 4, 2, 6, 1, 10, 7, 4, 9, 3, 4, 10, 1, 9, 3, 4, 9, 8)
(9, 10, 1, 10, 9, 10, 4, 6, 10, 2, 6, 7, 8, 8, 9, 4, 5, 4, 9, 8)
(3, 5, 4, 3, 1, 7, 10, 5, 2, 1, 7, 3, 10, 6, 6, 1, 1, 10, 6, 6)
  • Hi @xunilk, Thanks for your help. I used your code but unfortunately QGIS crashed again. Can I asked what is your version of QGIS and OS? – nickan Jul 12 '17 at 22:19
  • It could be your raster. Try with my 20x20 raster first. I will hang it in dropbox. Please, wait a moment. By the way, my QGIS is 2.18.10 and I use a Debian stable system. – xunilk Jul 12 '17 at 22:27
  • It's another raster but it is also 20 x20: dropbox.com/s/0el8ewzpf7rwjzd/aleatorio.tif?dl=0 – xunilk Jul 12 '17 at 22:32
  • Sorry @xunilk, I used your raster file but It crashed again. – nickan Jul 13 '17 at 6:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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