Retrieve pixel value with geographic coordinate as input with gdal

Question

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?

Answer 1

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.

Answer 2

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 =  \
        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]]

Answer 3

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]

Answer 4

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]

Answer 5

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.

<?php
date_default_timezone_set('America/Denver');
$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";  
        }
    fclose($handle);
}
?>