How to I get the coordinates of a cell in a geotif?

Question

I have a tif with geoinformation. With gdal I can transform the raster file to an array (numpy).

How can I get the coordinates for one entry in that array?

Answer 1

Use the affine transformation matrix, which maps pixel coordinates to world coordinates. So for example, use the affine package. (There are other ways to do the same, using simple math.)

from affine import Affine
fname = '/path/to/raster.tif'

Here are two ways to get the affine transformation matrix, T0. E.g., using GDAL/Python:

from osgeo import gdal
ds = gdal.Open(path, gdal.GA_ReadOnly)
T0 = Affine.from_gdal(*ds.GetGeoTransform())
ds = None  # close

E.g., using rasterio:

import rasterio
with rasterio.open(fname, 'r') as r:
    T0 = r.affine

The convention for the transform array as used by GDAL (T0) is to reference the pixel corner. You may want to instead reference the pixel centre, so it needs to be translated by 50%:

T1 = T0 * Affine.translation(0.5, 0.5)

And now to transform from pixel coordinates to world coordinates, multiply the coordinates with the matrix, which can be done with a simple function:

rc2xy = lambda r, c: T1 * (c, r)

Now, get the coordinates for a raster in the first row, second column (index [0, 1]):

print(rc2xy(0, 1))

Also, note that if you need to get the pixel coordinate from a world coordinate, you can use an inverted affine transformation matrix, ~T0.