關於GDAL 仿射變化參數和圖像投影

轉自： http://blog.csdn.net/ivan_ljf/article/details/9226463

 

關於GDAL計算圖像坐標的幾個問題

使用GDAL處理地理圖像時，不可避免的會遇到一個問題，圖像的地理坐標問題，因為有了這個地理坐標，地理圖像才和普通圖像有了最本質的區別，那么在使用GDAL時，如何處理與地理坐標相關的信息呢？下面進行簡單的說明。

1：如何使用行列號計算圖像的地理坐標？或者如何通過地理坐標來定位在圖像的某個位置？

2：如何獲取圖像的四至范圍？或者如果通過指定的地理范圍計算圖像的所在區域？

要解決上面三個問題，首先需要知道和了解GDAL的數據模型，其中里面有個非常重要的就是投影和六參數。這兩個可以使用GDALDataset類中的GeoTransform()函數和GetProjectionRef()函數來進行獲取。

第一個參數獲取的是圖像的六參數（我自己起的名字，是一個仿射變化的參數），第二個是圖像的投影（也就是空間參考系統）。

下面先說說第一個六參數，六參數其實是圖像行列號坐標和地理坐標轉換的一組轉換系數。下面是用GT來表示六參數，圖像行列號與圖像的地理坐標之間的數學關系式如下：

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

上式中，Xgeo和Ygeo表示的圖像的地理坐標，Xpixel表示圖像的列號，Yline表示圖像的行號，GT(i)就是上面所說的六參數，一共是六個值。

這六個值大致可以分為三組:

　　GT(0)和GT(3)是第一組，表示圖像左上角的地理坐標；

　　GT(1)和GT(5)是第二組，表示圖像橫向和縱向的分辨率（一般這兩者的值相等，符號相反，橫向分辨率為正數，縱向分辨率為負數）；

　　GT(2)和GT(4)是第三組，表示圖像旋轉系數，對於一般圖像來說，這兩個值都為0。

 

為什么說圖像的GT(0)和GT(3)表示圖像左上角的坐標，對於圖像行列號坐標系統來說，坐標的原點在左上角，所以左上角的行列號是(0,0)，將坐標帶入上式，可以得到：

    Xgeo = GT(0)
    Ygeo = GT(3)

所以說GT(0)和GT(3)表示圖像左上角的坐標。

 

GT(1)和GT(5)表示圖像橫向和縱向的分辨率。圖像的分辨率就是圖像每個像素所能表示的面積，一般都是正方形的格網，所以也就是沒兩個相鄰像元坐標的差值。

基於這個原理，使用兩個坐標進行驗證。假設當前點行列號坐標為A（i，j），相鄰的右側點坐標為B（i+1，j）。分別計算A和B的橫向地理坐標，並計算差值，即：

dX  = XgeoB - XgeoA

　　= [GT(0) + (i+1)*GT(1) + j*GT(2)] - [GT(0) + i*GT(1) + j*GT(2)]

　　=  GT(0) -GT(0) + (i+1)*GT(1) - i*GT(1) + j*GT(2) - j*GT(2) 

　　=  (i+1)*GT(1) - i*GT(1)

　　= GT(1)

同理可以得到 dY= GT(5)。

對於一個普通的標准圖像來說（這里的標准圖像是指GT(2)和GT(4)都為0），如圖1所示，圖像的行列號坐標為XOY，每個網格代表一個圖像像素區域，i表示列號，j表示行號，淡藍色右下角的行列坐標為（i，j），

圖中紅色方塊縱向長度為dy，橫向長度為dx，分別為圖像的分辨率；圖中O點的地理坐標就是（GT(0),GT(3)）。                   

   

圖1 一個標准的圖像行列號坐標及其地理坐標說明

有了上面的說明，那么就可以很簡單的來進行圖像的行列號與地理坐標進行相互轉換，具體的代碼如下，共有兩個，一個正算，一個反算。

 

bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow)
{
    try
    {
        double dTemp = adfGeoTransform[1]*adfGeoTransform[5] - adfGeoTransform[2]*adfGeoTransform[4];
        double dCol = 0.0, dRow = 0.0;
        dCol = (adfGeoTransform[5]*(dProjX - adfGeoTransform[0]) - 
            adfGeoTransform[2]*(dProjY - adfGeoTransform[3])) / dTemp + 0.5;
        dRow = (adfGeoTransform[1]*(dProjY - adfGeoTransform[3]) - 
            adfGeoTransform[4]*(dProjX - adfGeoTransform[0])) / dTemp + 0.5;

        iCol = static_cast(dCol);
        iRow = static_cast(dRow);
        return true;
    }
    catch(...)
    {
        return false;
    }
}

bool ImageRowCol2Projection(double *adfGeoTransform, int iCol, int iRow, double &dProjX, double &dProjY)
{
    //adfGeoTransform[6]  數組adfGeoTransform保存的是仿射變換中的一些參數，分別含義見下
    //adfGeoTransform[0]  左上角x坐標 
    //adfGeoTransform[1]  東西方向分辨率
    //adfGeoTransform[2]  旋轉角度, 0表示圖像 "北方朝上"
    //adfGeoTransform[3]  左上角y坐標 
    //adfGeoTransform[4]  旋轉角度, 0表示圖像 "北方朝上"
    //adfGeoTransform[5]  南北方向分辨率

    try
    {
        dProjX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;
        dProjY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;
        return true;
    }
    catch(...)
    {
        return false;
    }
}

 

現在我們再回到之前開頭的兩個問題。

對於第一個問題，實際就是圖像行列號坐標與地理坐標的相互轉換，上面的代碼就可以用來解決。對於第二個問題，可以轉換為第一個問題，第二個問題其實就是兩個點的坐標轉換，分別是左上角點和右下角點。

比如第一個，如何計算圖像的四至范圍，圖像的四至范圍從圖1中可以看出，圖像的四至其實就是圖像左上角坐標和右下角的坐標為起來的矩形區域，那么就分別將左上角和右下角的行列號按照上面的公式進行轉換即可得到四至范圍；對與第二問就是將這個坐標進行反算得到。

 

關於地理圖像的坐標問題就說到這里，關於上面的投影信息（空間參考信息），需要說明一下，這個六參數里面的坐標范圍是和空間參考是一一對應的，

比如空間參考是一個WGS84橢球體，那么這個六參數一般的單位就是度，如果是北京54等分帶投影，那么六參數的單位就是米。