This article shows how to read the value of a raster cell with QGIS and GDAL.
QGIS and GDAL both have Python bindings, you can use both libraries to read a value from a raster cell, since QGIS uses GDAL libraries under the hood, we can expect to read the exact same value with both systems. Here is a short example about how to do it with the two different approaches, we assume that you are working inside the QGIS python console and the project has a raster file loaded, but with just a few modifications, the example can also be run from a standard python console. The example raster layer is a DTM with 1000 cells width and 2000 cells height, we want to read the value at the cell with coordinates x = 500 and y = 1000.
# First layer in QGIS project is a DTM 2 bands raster
from osgeo import gdal
# You need this to convert raw values readings from GDAL
import struct
# Read the cell with this raster coordinates
x = 500
y = 1000
# Get the map layer registry
reg = QgsMapLayerRegistry.instance()
# Get the first layer (the DTM raster)
qgis_layer = reg.mapLayers().values()[0]
# Open the raster with GDAL
gdal_layer = gdal.Open(rlayer.source())
"""
Fetches the coefficients for transforming between pixel/line (P,L) raster space,
and projection coordinates (Xp,Yp) space.
Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];
In a north up image, padfTransform[1] is the pixel width, and padfTransform[5]
is the pixel height. The upper left corner of the upper left pixel is
at position (padfTransform[0],padfTransform[3]).
"""
gt = gldal_layer.GetGeoTransform()
# o:origin, r:rotation, s:size
xo, xs, xr, yo, yr, ys = gt
# Read band 1 at the middle of the raster ( x = 500, y = 1000)
band = gdal_layer.GetRasterBand(1)
gdal_value = struct.unpack('f', band.ReadRaster(x, y, 1, 1, buf_type=band.DataType))[0]
xcoo = xo + xs * x + xr * y
ycoo = yo + yr * x + ys * y
# Read the value with QGIS, we must pass the map coordinates
# and the exact extent = 1 cell size
qgis_value = qgis_layer.dataProvider().identify(QgsPoint(xcoo, ycoo), \
QgsRaster.IdentifyFormatValue, \
theExtent=QgsRectangle( xcoo , ycoo, xcoo + xs, ycoo + ys) )\
.results()[1]
assert(gdal_value == qgis_value)