如何获取Geotif中单元格的坐标?
问题内容:
我有地理信息提示。使用gdal,我可以将栅格文件转换为数组(numpy)。
如何获取该数组中一个条目的坐标?
问题答案:
使用仿射变换矩阵,该矩阵将像素坐标映射到世界坐标。因此,例如,使用仿射软件包。(还有其他方法可以使用简单的数学方法完成此操作。)
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
。