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

电玩女神 2022-08-20 15:23 26阅读 0赞

这里要使用的一个比较重要的函数是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 gdal
    
    gdal.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))