GDAL获取栅格数据各个像素对应的经纬度（Python版）

本文主要是介绍GDAL获取栅格数据各个像素对应的经纬度（Python版），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

这里要使用的一个比较重要的函数是GetGeoTransform函数，GDAL官网的API中对GetGeoTransform的解释如下：
CPLErr GDALDataset::GetGeoTransform ( double * padfTransform )
virtual
Fetch the affine transformation coefficients.

Fetches the coefficients for transforming between pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.

Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];

In a north up image, padfTransform[1] is the pixel width, and padfTransform[5] is the pixel height. The upper left corner of the upper left pixel is at position (padfTransform[0],padfTransform[3]).

The default transform is (0,1,0,0,0,1) and should be returned even when a CE_Failure error is returned, such as for formats that don’t support transformation to projection coordinates.

This method does the same thing as the C GDALGetGeoTransform() function.

Parameters
padfTransform an existing six double buffer into which the transformation will be placed.
Returns
CE_None on success, or CE_Failure if no transform can be fetched.
Reimplemented in GDALProxyPoolDataset, VRTDataset, GDALPamDataset, GDALProxyDataset, and GDALGeorefPamDataset.

我们将使用如下公式进行坐标的计算：

 Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

from osgeo import gdalgdal.AllRegister()filePath = '/home/theone/Data/GreatKhingan/DEM/Slope_GreatKhingan_500m.tif'
dataset = gdal.Open(filePath)adfGeoTransform = dataset.GetGeoTransform()# 左上角地理坐标
print(adfGeoTransform[0])
print(adfGeoTransform[3])nXSize = dataset.RasterXSize #列数
nYSize = dataset.RasterYSize #行数arrSlope = [] # 用于存储每个像素的（X，Y）坐标
for i in range(nYSize):row = []for j in range(nXSize):px = adfGeoTransform[0] + i * adfGeoTransform[1] + j * adfGeoTransform[2]py = adfGeoTransform[3] + i * adfGeoTransform[4] + j * adfGeoTransform[5]col = [px, py]row.append(col)arrSlope.append(row)print(len(arrSlope))