java+gdal實現影像重投影

java+gdal實現影像重投影

GDAL功能很強大，用來處理影像資料，今天我要做的是java程式碼寫的影像重投影，網上參考資料大都是c++和python寫的，也看了一些大牛寫的程式碼，最後寫出了java版的，eclipse寫的，直接引用一個gdal.jar包，不過要有一些dll檔案，網上有相關的java配置jdal庫的部落格，不配置jdal會報錯：本地庫錯誤。還有對於gdal讀取六引數geoTransform的理解，我在這方面沒理解好，耗費了一些時間，會在下文講到。我的源資料與目標資料具有相同的基準面，屬於嚴密轉換，實現起來較為簡單。實現重投影要需要以下幾個步驟：

• 源影像仿射變換系數及投影獲取
• 目標影像資訊寫入
• ReprojectImage重新投影
• 重取樣（本人需求）

參考博文： 
1.GDAL庫學習筆記(三):GDAL建立資料集 
2.GDAL之柵格重投影 
3.GDAL影像投影轉換

以上博文中，1和3是用C++寫的程式碼，2是用python寫的。關於背景知識以上博文講訴的較為詳細，這裡只做簡單描述：

影像要進行投影轉換，首先要自帶座標系，然後看源影像座標系統和目標座標系統。本文進行的是嚴密轉換，同一基準面下的地理座標系轉到投影座標系。

現將java程式碼中用的GDAl主要函式列出來，做簡要說明：

GetProjectionRef() 獲取源影像的座標參考 
GetGeoTransform() 讀取六引數，具體含義參考以上引用博文 
CloneGeogCS() 獲取一個投影系中的地理系，投影座標系是地理座標投影的結果， 
CoordinateTransformation ct = new CoordinateTransformation(src_Crs, oLatLong);兩種座標系統的變換關係 
ReprojectImage(Dataset src_ds, Dataset dst_ds, String src_wkt, String dst_wkt, int eResampleAlg)實現重投影的函式

先在eclipse中新建工程->包->類，匯入gdal.jar包，配置好gdal（實際執行在dll），以下是參考程式碼，註釋較為詳細，可以更改：

package readTifftest;

import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;
import org.gdal.osr.CoordinateTransformation;
import org.gdal.osr.SpatialReference;


/**
 * 
 * @author Administrator
 * @ClassName:Resample
 * @Version 版本
 * @ModifiedBy 修改人
 * @date 2018年1月16日 下午2:49:43
 */
public class GDAL_Resample {

    public static void main(String[] args) {
        // TODO Auto-generated method stub
        Resample("E:\\database\\overview\\J46_棄用\\J46D001010.png","E:\\database\\overview\\J46\\J46D001010.tif");

    }
    /**
     * 影像投影轉換並重取樣得到輸出影像，同一個橢球體基準面下的轉換就是嚴密的轉換
     * @param input  輸入影像路徑
     * @param output 輸出影像路徑
     */
    public static void Resample(String input,String output){
        //String agentfile=output.substring(0, output.lastIndexOf("."))+".tif";
        //註冊GDAL
        gdal.AllRegister();
        //只讀方式讀取資料
        Dataset pSrc = gdal.Open(input,gdalconstConstants.GA_ReadOnly);
        if (pSrc == null)
            return;
        //開啟的影像畫素、波段等資訊
        int numBands=pSrc.GetRasterCount(); //讀取影像波段數
        int xSize = pSrc.GetRasterXSize();//柵格尺寸
        int ySize = pSrc.GetRasterYSize();//
        int depth=pSrc.GetRasterBand(1).GetRasterDataType();//影象深度

        String src_wkt=pSrc.GetProjectionRef();//獲取源影象crs

        SpatialReference src_Crs=new SpatialReference(src_wkt);//構造投影座標系統的空間參考(wkt)

        double[] geoTransform = pSrc.GetGeoTransform();     
        //影象範圍
        double xmin = geoTransform[0];
        double ymax = geoTransform[3];
        double w_src = geoTransform[1];//畫素寬度
        double h_src = geoTransform[5];//畫素高度
        double xmax = geoTransform[0] + xSize * w_src;
        double ymin = geoTransform[3] + ySize * h_src;

        //設定輸出影象的座標 
        SpatialReference oLatLong;  
        oLatLong = src_Crs.CloneGeogCS(); //獲取該投影座標系統中的地理座標系   
        //構造一個從投影座標系到地理座標系的轉換關係  
        CoordinateTransformation ct = new CoordinateTransformation(src_Crs, oLatLong);

        // 此註釋為錯誤計算過程  
//        double dst1[]=ct.TransformPoint(xmin,ymax);
//        double dst2[]=ct.TransformPoint(xmax,ymin);
//      double dstxmin = dst1[0];
//      double dstymax = dst1[1];
//      double dstxmax = dst2[0];
//      double dstymin = dst2[1]; //39.67609572370327

        //計算目標影像的左上和右下座標,即目標影像的仿射變換引數,投影轉換為經緯度
        double a[]= ct.TransformPoint(xmin, ymax);
        double b[]= ct.TransformPoint(xmax, ymax);
        double c[]= ct.TransformPoint(xmax, ymin);
        double d[]= ct.TransformPoint(xmin, ymin);     
        double dbX[]={a[0],b[0],c[0],d[0]};
        double dbY[]={a[1],b[1],c[1],d[1]};

        GDAL_Resample test=new GDAL_Resample();
        test.res(dbX);//按從小到大排列
        double dstxmin=dbX[0];
        double dstxmax=dbX[3];
        test.res(dbY);
        double dstymin=dbY[0];
        double dstymax=a[1];
        double[] adfGeoTransform=new double[6];
        adfGeoTransform[0]=dstxmin;
        adfGeoTransform[3]=dstymin;
        adfGeoTransform[1]=(dbX[3]-dbX[0])/xSize;
        adfGeoTransform[5]=(dbY[3]-dbY[0])/ySize;
        //adfGeoTransform[2]=(dstxmax-adfGeoTransform[0]-xSize*adfGeoTransform[1])/ySize;
       // adfGeoTransform[4]=(dbY[3]-adfGeoTransform[3]-ySize*adfGeoTransform[5])/xSize;

        //建立輸出影象
        Driver driver = gdal.GetDriverByName("GTiff"); 
        driver.Register();  
        String[] options={"INTERLEAVE=PIXEL"};
        Dataset pDst = driver.Create(output,xSize,ySize,numBands,depth,options);

        //寫入仿射變換系數及投影 
        String dst_wkt=oLatLong.ExportToWkt();
        pDst.SetGeoTransform(adfGeoTransform);
        pDst.SetProjection(dst_wkt);

        /*eResampleAlg取樣模式
         * GRA_NearestNeighbour=0   最近鄰法，演算法簡單並能保持原光譜資訊不變；缺點是幾何精度差，灰度不連續，邊緣會出現鋸齒狀     
           GRA_Bilinear=1         雙線性法，計算簡單，影象灰度具有連續性且取樣精度比較精確；缺點是會喪失細節； 
           GRA_Cubic=2            三次卷積法，計算量大，影象灰度具有連續性且取樣精度高； 
           GRA_CubicSpline=3      三次樣條法，灰度連續性和取樣精度最佳； 
           GRA_Lanczos=4          分塊蘭索斯法，由匈牙利數學家、物理學家蘭索斯法創立，實驗發現效果和雙線性接近； 
         */
        //重新投影
        gdal.ReprojectImage(pSrc, pDst, src_wkt, dst_wkt, 3);

        pDst.delete();
        pSrc.delete();
    }
    /**
     * 比較大小，得到最大值與最小值
     * @param arr
     * @return
     */
      public double[] res(double[] arr){
            for(int i=0;i<arr.length;i++){
                for(int j=0;j<arr.length-i-1;j++){

                    if(arr[j]>arr[j+1]){
                        double temp=arr[j];
                        arr[j]=arr[j+1];
                        arr[j+1]=temp;
                    }
                }
            }                
            return arr;         
        }
}