如何获取geotif中单元格的坐标？

Question

我有一个地理信息的tif.使用gdal,我可以将光栅文件转换为数组(numpy).

如何获取该数组中一个条目的坐标？

Answer 1

使用仿射变换矩阵,将像素坐标映射到世界坐标.例如,使用仿射包.(还有其他方法可以使用简单的数学方法.)

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.