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() # 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 + P*padfTransform + L*padfTransform; Yp = padfTransform + P*padfTransform + L*padfTransform; In a north up image, padfTransform is the pixel width, and padfTransform is the pixel height. The upper left corner of the upper left pixel is at position (padfTransform,padfTransform). """ 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)) 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() assert(gdal_value == qgis_value)