I have a list of X,Y coordinates in UTM called coords. I also have a .tif of a digital terrain model (DTM) also referenced in UTM. I would like to use the Python wrapper for gdal to return the pixel values (i.e. elevation values) for each coordinate in coords.

In my search for an answer so far I have only found this answer, which is related, but does not just use a simple list of coordinates and is seemingly a bit more complex than the simple task I am trying to achieve.

My concern is not so much the particular script that needs to be written, but the methods or functions within gdal that can perform this type of function. Is there a simple function within gdal that takes a geographic coordinate and returns a pixel value?

5 Answers 5


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,

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.

  • 1
    how can I optimize when I have large data array? I do not have any point, I need to traverse each cell in nested for loop to check the raster value. In my case, it is taking several hours because of huge data array Feb 29, 2020 at 3:13

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.

from osgeo import gdal 
import affine
import numpy as np

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 =  \
    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]]
  • 3
    Not only does this work well, it also handles a geotransformation where the skew is non-zero. The above, accepted answer will not work properly if elements 2 and 4 (xskew, yskew in the GDAL documentation) are non-zero.
    – Jay Laura
    May 9, 2017 at 19:34
  • In some coordinate systems, notably EPSG:3857 WGS84, the coordinate order is lat,lon or y,x, which confuses this routine with a transposed data_array. See gis.stackexchange.com/questions/124077/…
    – Dave X
    Aug 9, 2021 at 19:50
  • To get consistent data out of EPSG:4326 datasets such as the Geoid models from usna.edu/Users/oceano/pguth/md_help/html/egm96.htm , I needed to transpose the numpy array, data_array = np.array(data_source.GetRasterBand(1).ReadAsArray().transpose())
    – Dave X
    Aug 9, 2021 at 21:29
  • 1
    If you use coordinates outside the bounds of the data to the left or up on the raster, this can deliver nonsense results because of numpy negative indexing. For example, coordinates mapping to pixel (-1,-1) ends up delivering the value for the lower right corner.
    – Dave X
    Aug 9, 2021 at 22:21
  • Be awere that the function doesn't work well with negative coords. It returns (0,0) pixel_coord when it's for example (-0.99, -0.99) because int() doesn't round, uses first digits till dot. You should use round() instead! also adding this 0.5 may change the result to incorrect, for example when pixel_coords without rounding would be (-0.90 , -0.90), it should end up in (-1,-1) not (0, 0). This is important when you want to extract for example misalignment of top left corner between 2 rasters. Apr 20, 2023 at 8:53

The basic steps are as follows:

  1. Find out coordinate of the top left pixel in whichever CRS the raster is projected in, which is stored in the raster metadata
  2. Find out the size of each pixel in the x and y dimensions
  3. Use these pieces of information to convert the coordinate to an index of the raster grid (i.e. the row and column of the desired pixel), and sample this pixel.

In python:

from osgeo import gdal

raster_file = "raster.tif"
ds = gdal.Open(raster_file)

# GetGeoTransform gives the (x, y) origin of the top left pixel,
# the x and y resolution of the pixels, and the rotation of the
# raster. If the raster is rotated (i.e. the rotation values are
# anything other than 0) this method will not work.

# The information is returned as tuple:
# (TL x, X resolution, X rotation, TL y, Y rotation, y resolution)
TL_x, x_res, _, TL_y, _, y_res = gds.GetGeoTransform()

# The point where you wish to sample the raster
coordinate = (x, y)

# Divide the difference between the x value of the point and origin,
# and divide this by the resolution to get the raster index
x_index = (x - TL_x) / x_res

# and the same as the y
y_index = (y - TL_y) / y_res

# Read the raster as an array
array = ds.ReadAsArray()

# Sample with the indexs, not that y_index should be first as the index is
# [rows, columns] in a 2d grid in python
pixel_val = array[y_index, x_index]

Here you can find solution which doesn't use any external library other than GDAL and reads pixel value without loading the whole raster into memory.

from osgeo import gdal 
import math

def get_raster_value(geo_x: float, geo_y: float, ds: gdal.Dataset, band_index: int = 1):
    """Return raster value that corresponds to given coordinates."""
    forward_transform = ds.GetGeoTransform()
    reverse_transform = gdal.InvGeoTransform(forward_transform)
    pixel_coord = gdal.ApplyGeoTransform(reverse_transform, geo_x, geo_y)
    pixel_x = math.floor(pixel_coord[0])
    pixel_y = math.floor(pixel_coord[1])
    band: gdal.Band = ds.GetRasterBand(band_index)
    val_arr = band.ReadAsArray(pixel_x, pixel_y, 1, 1) # Avoid reading the whole raster into memory - read 1x1 array
    return val_arr[0][0]
  • what are pixel_width and pixel_height used for? Feb 9, 2023 at 1:01
  • 1
    Good point @Mike'Pomax'Kamermans, pixel_width and pixel_height are not used in this example. Code updated. Feb 10, 2023 at 10:25
  • cheers. I've (suggested an) update to the code to also make x and y "be the right value" immediately rather than needing separate floors. Feb 10, 2023 at 16:53

This script is not Python but could easily be adapted. It shows the requested gdal command which will return a pixel value from a GeoTIFF based on X,Y query. The projection would be irrelevant, as long as they match (CVS input txt and GeoTIFF projection).

The input CSV txt file would have 3 columns with one point per line: POINT NAME, X, Y

Edit: For direct GDAL & Python interaction, check out Shapely and its built in GDAL functions.

$incomingfile = "./DroneMapper-GCPCheckXY.txt";

if (($handle = fopen($incomingfile, "r")) !== FALSE) {
    while (($data = fgetcsv($handle, 1000, ",")) !== FALSE) {
        $value = shell_exec("/usr/local/bin/gdallocationinfo -geoloc -valonly folder/myGeoTIFF.tif ".$data[1]." ".$data[2]);
        $delta = $value - $data[3];
        print "GCP: " . $data[0] . " Z Delta: " . $delta . "\n";  

Your Answer

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

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