6

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?

9

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.

7

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]]
  • 2
    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 '17 at 19:34
0

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);
}
?>
0

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]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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