Retrieve pixel value with geographic coordinate as input with gdal

You can retrieve raster pixel values with any of the following scripts. List of X,Y coordinates (as tuples) are in the python list named 'points_list'.

With Python GDAL:

from osgeo import gdal

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/aleatorio.tif" #path to raster
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

transform = dataset.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

data = band.ReadAsArray(0, 0, cols, rows)

points_list = [ (355278.165927, 4473095.13829), (355978.319525, 4472871.11636) ] #list of X,Y coordinates

for point in points_list:
    col = int((point[0] - xOrigin) / pixelWidth)
    row = int((yOrigin - point[1] ) / pixelHeight)

    print row,col, data[row][col]

With PyQGIS:

filename = "/home/zeito/pyqgis_data/aleatorio.tif" #path to raster

layer = QgsRasterLayer(filename,
                       "my_raster")

provider = layer.dataProvider()

extent = layer.extent()

xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

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

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

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

points_list = [ (355278.165927, 4473095.13829), (355978.319525, 4472871.11636) ]#list of X,Y coordinates

for point in points_list:
    col = int((point[0] - xmin) / pixelWidth)
    row = int((ymax - point[1] ) / pixelHeight)

    print row,col, block.value(row, col)

I tried them with my particular raster and they worked. Result was, for both cases, the following:

4 4 36
7 13 42

The first and second valor (each line) are indices of row, column (for verification purposes). The third one is raster value.

Here is the function I came up with, using a function I found in another stack post (that I unfortunately cannot remember the title of). It was originally written to be used with a point vector file instead of manually inputting the points like I am doing. Below is the simplified function, using affine and gdal, where data_source is an opened gdal object of a GeoTIFF and coord is a tuple of a geo-coordinate. This tuple must be in the same coordinate system as the GeoTIFF.

def retrieve_pixel_value(geo_coord, data_source):
    """Return floating-point value that corresponds to given point."""
    x, y = geo_coord[0], geo_coord[1]
    forward_transform =  \
        affine.Affine.from_gdal(*data_source.GetGeoTransform())
    reverse_transform = ~forward_transform
    px, py = reverse_transform * (x, y)
    px, py = int(px + 0.5), int(py + 0.5)
    pixel_coord = px, py

    data_array = np.array(data_source.GetRasterBand(1).ReadAsArray())
    return data_array[pixel_coord[0]][pixel_coord[1]]