Как получить координаты ячейки в геотиме?

Question

У меня есть tif с геоинформацией. С gdal я могу преобразовать растровый файл в массив (numpy).

Как получить координаты для одной записи в этом массиве?

Answer 1

Используйте матрицу аффинного преобразования, которая отображает координаты пикселей в мировые координаты. Например, используйте пакет affine. (Есть и другие способы сделать то же самое, используя простую математику.)

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

Вот два способа получить матрицу аффинного преобразования, T0. Например, с использованием GDAL/Python:

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

например, с использованием rasterio:

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

В конвенции для преобразования массива, которые используются GDAL (T0), чтобы ссылаться на пиксельный угол. Вы можете вместо ссылки на пиксель центра, поэтому он должен быть переведен на 50%:

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

И теперь, чтобы преобразовать от пиксельных координат на мировые координаты, умножать координаты с матрицей, которая может быть сделана с простая функция:

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

Теперь, получить координаты растра в первой строке, второй столбец (индекс [0, 1]):

print(rc2xy(0, 1)) 

Кроме того, обратите внимание: если вам нужно получить координату пикселя из мировой координаты, вы можете использовать матрицу инвертированного аффинного преобразования, ~T0.