Java 調用gdal API(二)——柵格裁剪

gdal可以說是GIS數據處理比較好的工具之一，雖然也提供了Java API，但是官方文檔確實太過簡單，用起來確實太難受，每次都需要去參考對應的C++api，然後在對應使用。

因此小編決定從這篇文章開始，將自己以前以及今後用到的API，都記錄下，以幫助更多從事java GIS服務開發者。

今天主要來記錄柵格數據的裁剪功能，gdal中提供了多種裁剪方法，今天主要介紹以下兩種：

第一中是通過gdal自身的Warp方法，其實現有以下幾個介面；

其實這個幾個介面基本差不多

1、dstDs/dest 其實就是目標數據，兩種方式，要麼傳入目標路徑，要麼是目標Dataset

2、object_list_count 這個是源數據的Dataset 數據

3、warpAppOptions 這個參數很重要，是我們裁剪需要傳遞的重要參數，後面代碼會專門解釋

4、callback 處理進度

下麵我們結合最後一個方法，來處理一個柵格裁剪

private static void cutTif() throws FactoryException, IOException, ParseException {
        Dataset dataSet = gdal.Open("E:\\數據\\影像數據\\tiff\\xian18.tif");
        //通過最大最小經緯度生成一個面矢量
        createCutShapeRegion(108.583740,33.928528,108.883740,34.260242);
        Vector vector = new Vector();
        vector.add("-cutline");
        vector.add("D:\\data\\cut2.shp");
        vector.add("-crop_to_cutline");
        Dataset ds = gdal.Warp("D:\\data\\xian18_1.tif",new Dataset[]{dataSet},new WarpOptions(vector));
        if (ds != null){
            System.out.println("成功！");
        }
    }

createCutShapeRegion是用通過bbox來生成shape數據的，我這邊是用矩形來裁剪的，當然你也可以用多邊形，通過傳入多邊形的頂點來生層一個面矢量，記住一定是面矢量，

D:\\data\\cut2.shp是createCutShapeRegion生成的矢量文件，最主要的兩個參數一個是-cutline 一個是-crop_to_cutline，其中-cutline後面跟裁剪的矢量文件，-crop_to_cutline不需要具體值，這些具體參數可參考gdal的warp的參數介紹，其實和安裝版gdalwarp的參數是基本一致的。

但是這裡說明下，並不是gdalwarp的所有參數都在gdal的java api中進行了適配，就比如gdalwarp這個命令中除了上面使用shp作為裁剪之外，還有一個-te參數，這個參數後面跟一個經緯度的範圍
主要作用就是將原始數據裁剪為制定的範圍，就比如我們使用安裝的gdal中可以這樣調用裁剪
gdalwarp.exe -te 108.583740 33.928528 108.883740 34.260242  E:\\數據\\影像數據\\tiff\\xian18.tif D:\\data\\xian18_2.tif
然而這個-te參數如果使用代碼來調用，就會報沒有-te這個option，就比如下麵的這個例子其實是跑不通的

public static void main(String[] args) throws NoSuchFieldException, IllegalAccessException, FactoryException, IOException, ParseException {
        gdal.AllRegister();
        Dataset dataSet = gdal.Open("E:\\數據\\影像數據\\tiff\\xian18.tif");
        Vector vector = new Vector();
        vector.add("-te");
        vector.add("108.583740 33.928528 108.883740 34.260242"); //裁剪範圍
        vector.add("-crop_to_cutline");
        Dataset ds = gdal.Warp("D:\\data\\xian18_1.tif",new Dataset[]{dataSet},new WarpOptions(vector));
    }

上述代碼是跑不通的，雖然安裝板的gdalwarp可以使用-te來裁剪，而且效率比上面使用shp要更快，但是java代碼的api卻沒有集成這個-te參數，當然我不知道是我調用的問題，還是確實沒有適配

如果有瞭解的朋友，麻煩告知一聲。

除了上述使用gdalwarp的介面外，我們還可以通過使用driver來快速裁剪出一個tif,具體代碼如下

private static void clipRaster(Dataset dataset, Envelope boundingBox){
        //獲取當前數據範圍
        ReferencedEnvelope originalEnvelope = (ReferencedEnvelope) GdalDataAnalysis.readBboxFromFile(dataset);
        //判斷是否相交，不想交則直接返回
        if (!originalEnvelope.intersects(boundingBox)){
            return;
        }
        //求請求與原始數據的範圍
        ReferencedEnvelope interEnvelope = originalEnvelope.intersection(boundingBox);
        //原始數據上偏移位置
        int offset_x = (int) ((interEnvelope.getMinX() - originalEnvelope.getMinX()) / GdalDataAnalysis.getResoulution(dataset,true));
        int offset_y = (int) ((interEnvelope.getMaxY() - originalEnvelope.getMaxY()) / GdalDataAnalysis.getResoulution(dataset,false));
        //需要多少個原始像素
        int block_xsize = (int) ((interEnvelope.getMaxX() - interEnvelope.getMinX()) / GdalDataAnalysis.getResoulution(dataset,true));
        int block_ysize = (int)Math.abs((interEnvelope.getMaxY() - interEnvelope.getMinY()) / GdalDataAnalysis.getResoulution(dataset,false));
        //實際瓦片的長度，裁剪數據的解析度和原始圖像一致
        int image_Xbuf = block_xsize;
        int image_Ybuf =  block_ysize;
        //寫入瓦片的偏移位置
        int imageOffSetX = 0;
        int imageOffSetY = 0;
        //迴圈讀寫
        Driver driver = gdal.GetDriverByName("GTiff");
        int dataType = dataset.GetRasterBand(1).getDataType();
        int bands = dataset.getRasterCount();
        byte[] pixelArray = new byte[image_Xbuf*image_Ybuf];
        Dataset dsDataset = driver.Create("D:\\data\\xian18_1.tif",image_Xbuf,image_Ybuf,bands, dataType);
        for (int i = 1;i<=bands;i++){
            Band band = dataset.GetRasterBand(i);
            band.ReadRaster(offset_x,offset_y,block_xsize,block_ysize,image_Xbuf,image_Ybuf,dataType,pixelArray);
            //寫數據
            dsDataset.GetRasterBand(i).WriteRaster(imageOffSetX,imageOffSetY,block_xsize,block_ysize,pixelArray);
            //清理刪除緩存
            dsDataset.GetRasterBand(i).FlushCache();
            dsDataset.GetRasterBand(i).delete();
        }
        dsDataset.SetProjection(dataset.GetProjection());
        double[] rsGeoTransform = dataset.GetGeoTransform();
        rsGeoTransform[0] = interEnvelope.getMinX();
        rsGeoTransform[3] = interEnvelope.getMaxY();
        dsDataset.SetGeoTransform(rsGeoTransform);
        dsDataset.FlushCache();
        dsDataset.delete();
    }

這個函數中dataset是原始數據，boundingBox是裁剪矩形範圍，當然這個方法，無法裁剪不規則多邊形，也就是它只能裁剪矩形規則範圍。

相比使用shape文件裁剪效率，這個直接拷貝像素的效率更高，當然這些裁剪方法中，還是直接使用gdalwarp.exe的-te參數效率最高，但是

奈何gdal的java api中好像沒有適配這個參數，不知道為什麼？

好了今天的內容就講解到這裡，如有問題，歡迎評論區討論。