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