1. 程式人生 > >java+gdal實現影像重投影

java+gdal實現影像重投影

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;         
        }
}