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

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 

13759 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。